我们从Python开源项目中，提取了以下20个代码示例，用于说明如何使用scipy.signal.medfilt()。
def calculate_features_for_VAD(sound_frames, frequencies_axis, spectrogram): features = numpy.empty((spectrogram.shape[0], 3)) # smooted_spectrogram, smoothed_frequencies_axis = smooth_spectrogram(spectrogram, frequencies_axis, 24) for time_ind in range(spectrogram.shape[0]): mean_spectrum = spectrogram[time_ind].mean() if mean_spectrum > 0.0: sfm = -10.0 * math.log10(stats.gmean(spectrogram[time_ind]) / mean_spectrum) else: sfm = 0.0 # max_freq = smoothed_frequencies_axis[smooted_spectrogram[time_ind].argmax()] max_freq = frequencies_axis[spectrogram[time_ind].argmax()] features[time_ind][0] = numpy.square(sound_frames[time_ind]).mean() features[time_ind][1] = sfm features[time_ind][2] = max_freq """medfilt_order = 3 for feature_ind in range(features.shape[0]): features[feature_ind] = signal.medfilt(features[feature_ind], medfilt_order)""" return features
def plot_numvisc(diagfile): plt.figure() nc = Dataset(diagfile) t=nc.variables['t'][:] ke=nc.variables['ke'][:] dkdt=np.diff(ke)/np.diff(t) ens=nc.variables['enstrophy'][:] ensm=0.5*(ens[1:]+ens[:-1]) # deltake[visc,res]=-(ke[-1]-ke[0]) # deltaens[visc,res]=max(medfilt(ens,21))-ens[5] visc_tseries = -dkdt/ensm*4.4*np.pi visc_num = max(visc_tseries[t[1:]>0.02]) #print('N=%4i / visc = %4.1e / num = %4.2e'%(N[res],Kdiff[visc],visc_num[res])) plt.semilogy(t[1:],visc_tseries) plt.xlabel('time') plt.ylabel('viscosity (-(1/2V)dE/dt)') plt.grid('on') plt.show()
def plot_singlewave(file): sensor = import_sensorfile(file) sensor_processed = process_input(sensor) timestamps = [] [timestamps.append(str(item[0])) for item in sensor] sensor_processed_calibrated = mf.calibrate_median(sensor_processed) sensor_filtered = mf.medfilt(sensor_processed_calibrated, 3) plt.plot(sensor_filtered, linewidth='0.8') plt.xlim([0, 12000]) plt.ylim([-5, 5]) plt.ylabel('Acceleration (g)') plt.xlabel('Time (ms)') # plt.xticks(sensor_filtered, timestamps, rotation='vertical') plt.show()
def get_master_sky(wavelength, spectrum, fiber_to_fiber, path, exptime): masterwave = [] masterspec = [] for wave, spec, ftf in zip(wavelength,spectrum,fiber_to_fiber): masterwave.append(wave) masterspec.append(np.where(ftf>1e-8, spec/ftf, 0.0)) masterwave = np.hstack(masterwave) ind = np.argsort(masterwave) masterwave[:] = masterwave[ind] masterspec = np.hstack(masterspec) masterspec[:] = masterspec[ind] mastersky = medfilt(masterspec, 281) wv = np.arange(masterwave.min(),masterwave.max()+0.05,0.05) s = np.zeros((len(wv),2)) s[:,0] = wv s[:,1] = np.interp(wv, masterwave, mastersky / exptime) np.savetxt(op.join(path, 'sky_model.txt'), s)
def process(self, X): onset_func = zscore(self.func(X)) onset_func = signal.filtfilt(self.moving_avg_filter, 1, onset_func) onset_func = onset_func - signal.medfilt( onset_func[:,np.newaxis], self.median_kernel )[:,0] peaks = signal.argrelmax(onset_func) onsets = peaks[0][np.where(onset_func[peaks[0]] > self.threshold)] return onsets
def process(self, S): mag = S.magnitude() H = np.empty_like(mag) H[:] = signal.medfilt(mag, self.h_kernel) P = np.empty_like(mag) P[:] = signal.medfilt(mag, self.p_kernel) h_mask = self.mask_class(H, P, self.mask_exp) p_mask = self.mask_class(P, H, self.mask_exp) return (S * h_mask, S * p_mask)
def load_data(indir, smooth, bin_size): datas = [] infiles = glob.glob(os.path.join(indir, '*.monitor.csv')) for inf in infiles: with open(inf, 'r') as f: f.readline() f.readline() for line in f: tmp = line.split(',') t_time = float(tmp[2]) tmp = [t_time, int(tmp[1]), float(tmp[0])] datas.append(tmp) datas = sorted(datas, key=lambda d_entry: d_entry[0]) result = [] timesteps = 0 for i in range(len(datas)): result.append([timesteps, datas[i][-1]]) timesteps += datas[i][1] if len(result) < bin_size: return [None, None] x, y = np.array(result)[:, 0], np.array(result)[:, 1] if smooth == 1: x, y = smooth_reward_curve(x, y) if smooth == 2: y = medfilt(y, kernel_size=9) x, y = fix_point(x, y, bin_size) return [x, y]
def median_baseline(intensities, window_size=501): '''Perform median filtering baseline removal. Window should be wider than FWHM of the peaks. "A Model-free Algorithm for the Removal of Baseline Artifacts" Friedrichs 1995 ''' # Ensure the window size is odd if window_size % 2 == 0: window_size += 1 # Enable batch mode if intensities.ndim == 2: window_size = (1, window_size) return medfilt(intensities, window_size)
def ContourSet(img): med_img = medfilt(img,kernel_size=11) cs = plt.contour(med_img, np.arange(0.1,2,0.1), origin='lower') return cs
def medfilt_along_axis(x, N, axis=-1): """Applies median filter smoothing on one axis of an N-dimensional array. """ kernel_size = np.array(x.shape) kernel_size[:] = 1 kernel_size[axis] = N return medfilt(x, kernel_size) # TO UTILS
def shift_to_f0(v_shift, v_voi, fs, out='f0', b_smooth=True): v_f0 = v_voi * fs / v_shift.astype('float64') if b_smooth: v_f0 = v_voi * signal.medfilt(v_f0) if out == 'lf0': v_f0 = la.f0_to_lf0(v_f0) return v_f0 #==============================================================================
def format_for_modelling(m_mag, m_real, m_imag, v_f0, fs, nbins_mel=60, nbins_phase=45): # alpha: alpha = define_alpha(fs) # f0 to Smoothed lf0: v_voi = (v_f0>0).astype('float') v_f0_smth = v_voi * signal.medfilt(v_f0) v_lf0_smth = la.f0_to_lf0(v_f0_smth) # Mag to Log-Mag-Mel (compression): m_mag_mel = la.sp_mel_warp(m_mag, nbins_mel, alpha=alpha, in_type=3) m_mag_mel_log = la.log(m_mag_mel) # Phase feats to Mel-phase (compression): m_imag_mel = la.sp_mel_warp(m_imag, nbins_mel, alpha=alpha, in_type=2) m_real_mel = la.sp_mel_warp(m_real, nbins_mel, alpha=alpha, in_type=2) # Cutting phase vectors: m_real_mel = m_real_mel[:,:nbins_phase] m_imag_mel = m_imag_mel[:,:nbins_phase] # Removing phase in unvoiced frames ("random" values): m_real_mel = m_real_mel * v_voi[:,None] m_imag_mel = m_imag_mel * v_voi[:,None] # Clipping between -1 and 1: m_real_mel = np.clip(m_real_mel, -1, 1) m_imag_mel = np.clip(m_imag_mel, -1, 1) return m_mag_mel_log, m_real_mel, m_imag_mel, v_lf0_smth
def median_filter(data): return medfilt(data, kernel_size=3)
def calibrate_eye(self,eye_channel,realign_mark,realign_timebin,eye_medfilt_win=21,eye_gausfilt_sigma=3): ''' Args eye_channel (list): the first element is the channel name for the horizontal eye position the second element is the channel name for the vertial eye position realign_mark (string): event marker used to align eye positions realign_timebin (list): a period of time relative to the realign_mark. Example: [0,100] eye_medfilt_win (int): parameter for the median filter to smooth the eye movement trajectory eye_gausfilt_sigma (int): sigma of the gaussian kernel to smooth the eye movement trajectory Return: - ''' samp_time = 1000.0/self.sampling_rate[eye_channel[0]] # medfilt eye x, y position lamb_medfilt = lambda ite:signal.medfilt(ite,eye_medfilt_win) self.data_df[eye_channel[0]] = self.data_df[eye_channel[0]].apply(lamb_medfilt) self.data_df[eye_channel[1]] = self.data_df[eye_channel[1]].apply(lamb_medfilt) # gaussian filt eye x,y position lamb_gausfilt = lambda ite:ndimage.filters.gaussian_filter1d(ite,eye_gausfilt_sigma) self.data_df[eye_channel[0]] = self.data_df[eye_channel[0]].apply(lamb_gausfilt) self.data_df[eye_channel[1]] = self.data_df[eye_channel[1]].apply(lamb_gausfilt) # align eye to realign_mark, realign_timebin uses realign_mark as reference realign_poinnum = (self.data_df[realign_mark]/samp_time).values start_points = realign_poinnum + realign_timebin[0]/samp_time points_num = int((realign_timebin[1]-realign_timebin[0])/samp_time) for channel in eye_channel: align_points = list() for idx in self.data_df.index: start_point = start_points[idx] if ~np.isnan(start_point): start_point = int(start_point) end_point = start_point + points_num align_point = self.data_df[channel].loc[idx][start_point:end_point] align_point = align_point.mean() else: align_point = np.nan align_points.append(align_point) self.data_df[channel] = self.data_df[channel] - align_points # find all saccades for all trials
def classify_active(self): """Creates a classification for the active song using classifier """ self.logger.info('Classifying {0}'.format(str(self.active_song))) batch_size = self.params.get('batch_size', 100) indices = np.arange(self.active_song.time.size) input = self.get_data_sample(indices) prbs = self.classifier.predict_proba( input, batch_size=batch_size, verbose=1).T # for i in range(prbs.shape[1]): # print self.active_song.time[i], ':', prbs[:, i] #new_prbs = self.probs_to_classes(prbs) # print new_prbs.shape # for i in range(new_prbs.shape[1]): # print self.active_song.time[i], ':', new_prbs[:, i] unfiltered_classes = self.probs_to_classes(prbs) try: power_threshold = self.params['power_threshold'] except KeyError: thresholded_classes = unfiltered_classes else: self.logger.debug('Thresholding at {0} dB'.format(power_threshold)) below_threshold = np.flatnonzero( 10 * np.log10(self.active_song.power) < power_threshold) self.logger.debug( '{0} indices found with low power'.format(below_threshold.size)) thresholded_classes = unfiltered_classes thresholded_classes[below_threshold] = 0 # no need to be wasteful, filter if there is a filter try: medfilt_time = self.params['medfilt_time'] except KeyError: filtered_classes = thresholded_classes else: dt = self.active_song.time[1] - self.active_song.time[0] windowsize = int(np.round(medfilt_time / dt)) windowsize = windowsize + (windowsize + 1) % 2 filtered_classes = signal.medfilt(thresholded_classes, windowsize) self.active_song.classification = filtered_classes
def plot_example(missed, acknowledged): sensor_miss = import_sensorfile(missed) sensor_ack = import_sensorfile(acknowledged) # Window data mag_miss = window_data(process_input(sensor_miss)) mag_ack = window_data(process_input(sensor_ack)) # Window data mag_miss = window_data(process_input(sensor_miss)) mag_ack = window_data(process_input(sensor_ack)) # Filter setup kernel = 15 # apply filter mag_miss_filter = sci.medfilt(mag_miss, kernel) mag_ack_filter = sci.medfilt(mag_ack, kernel) # calibrate data mag_miss_cal = mf.calibrate_median(mag_miss) mag_miss_cal_filter = mf.calibrate_median(mag_miss_filter) mag_ack_cal = mf.calibrate_median(mag_ack) mag_ack_cal_filter = mf.calibrate_median(mag_ack_filter) # PLOT sns.set_style("white") current_palette = sns.color_palette('muted') sns.set_palette(current_palette) plt.figure(0) # Plot RAW missed and acknowledged reminders ax1 = plt.subplot2grid((2, 1), (0, 0)) plt.ylim([-1.5, 1.5]) plt.ylabel('Acceleration (g)') plt.plot(mag_miss_cal, label='Recording 1') plt.legend(loc='lower left') ax2 = plt.subplot2grid((2, 1), (1, 0)) # Plot Missed Reminder RAW plt.ylim([-1.5, 1.5]) plt.ylabel('Acceleration (g)') plt.xlabel('t (ms)') plt.plot(mag_ack_cal, linestyle='-', label='Recording 2') plt.legend(loc='lower left') # CALC AND SAVE STATS stats_one = sp.calc_stats_for_data_stream_as_dictionary(mag_miss_cal) stats_two = sp.calc_stats_for_data_stream_as_dictionary(mag_ack_cal) data = [stats_one, stats_two] write_to_csv(data, 'example_waves') plt.show()
def f_psd(x, Fs, method, Hzmed=0, welch_params={'window':'hanning','nperseg':1000,'noverlap':None}): ''' Calculate the power spectrum of a signal Parameters ---------- x : array temporal signal Fs : integer sampling rate method : str in ['fftmed','welch'] Method for calculating PSD Hzmed : float relevant if method == 'fftmed' Frequency width of the median filter welch_params : dict relevant if method == 'welch' Parameters to sp.signal.welch Returns ------- f : array frequencies corresponding to the PSD output psd : array power spectrum ''' if method == 'fftmed': # Calculate frequencies N = len(x) f = np.arange(0,Fs/2,Fs/N) # Calculate PSD rawfft = np.fft.fft(x) psd = np.abs(rawfft[:len(f)])**2 # Median filter if Hzmed > 0: sampmed = np.argmin(np.abs(f-Hzmed/2.0)) psd = signal.medfilt(psd,sampmed*2+1) elif method == 'welch': f, psd = sp.signal.welch(x, fs=Fs, **welch_params) else: raise ValueError('input for PSD method not recognized') return f, psd
def lidarCB(self, msg): if not self.received_data: rospy.loginfo("success! first message received") # populate cached constants if self.direction == RIGHT: center_angle = -math.pi / 2 self.direction_muliplier = -1 else: center_angle = math.pi / 2 self.direction_muliplier = 1 self.min_angle = center_angle - FAN_ANGLE self.max_angle = center_angle + FAN_ANGLE self.laser_angles = (np.arange(len(msg.ranges)) * msg.angle_increment) + msg.angle_min self.data = msg.ranges values = np.array(msg.ranges) # remove out of range values ranges = values[(values > msg.range_min) & (values < msg.range_max)] angles = self.laser_angles[(values > msg.range_min) & (values < msg.range_max)] # apply median filter to clean outliers filtered_ranges = signal.medfilt(ranges, MEDIAN_FILTER_SIZE) # apply a window function to isolate values to the side of the car window = (angles > self.min_angle) & (angles < self.max_angle) filtered_ranges = filtered_ranges[window] filtered_angles = angles[window] # convert from polar to euclidean coordinate space self.ys = filtered_ranges * np.cos(filtered_angles) self.xs = filtered_ranges * np.sin(filtered_angles) self.fit_line() self.compute_pd_control() # filter lidar data to clean it up and remove outlisers self.received_data = True if PUBLISH_LINE: self.publish_line() if SHOW_VIS: # cache data for development visualization self.filtered_ranges = filtered_ranges self.filtered_angles = filtered_angles
def get_wavelength_offsets(fibers, binsize, wvstep=0.2, colsize=300): bins = np.hstack([np.arange(0,fibers[0].D,binsize),fibers[0].D]) mn = np.min([fiber.wavelength.min() for fiber in fibers]) mx = np.min([fiber.wavelength.max() for fiber in fibers]) yf = [] wv = np.arange(mn,mx,wvstep) pd = [] nbins = len(bins)-1 for fiber in fibers: y = fiber.spectrum w = fiber.wavelength x = np.arange(fiber.D) s = -999*np.ones((nbins,fiber.D)) for i in np.arange(nbins): x0 = int(np.max([0,(bins[i]+bins[i+1])/2 - colsize/2])) x1 = int(np.min([1032,(bins[i]+bins[i+1])/2 + colsize/2])) s[i,x0:x1], flag = polynomial_normalization(x, y, x0, x1) s = np.ma.array(s,mask=s==-999.) ym = s.mean(axis=0) ym = biweight_filter(ym, 5, ignore_central=1) pr, ph = find_maxima(w, -1.*(y/ym-1)+1, interp_window=3) pd.append(pr[~np.isnan(pr)]) yf.append(np.interp(wv,w,-1.*(y/ym-1)+1,left=0.0,right=0.0)) master = biweight_location(yf,axis=(0,)) pr,ph = find_maxima(wv, master, y_window=10, interp_window=2.5, repeat_length=3, first_order_iter=5, max_to_min_dist=15) sel = (ph>1.03) * (~np.isnan(pr)) pr = pr[sel] ph = ph[sel] wk = [] for p,fiber in zip(pd,fibers): c = p[:,None] - pr a = np.where(np.abs(c)<4.) x = pr[a[1]] y = p[a[0]] - pr[a[1]] sel = np.argsort(x) x = x[sel] y = y[sel] m = medfilt(y,15) good = ~is_outlier(y-m) p0 = np.polyfit(x[good],y[good],3) wk.append(np.polyval(p0,fiber.wavelength)) return wk
