python类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 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