Python scipy.signal.medfilt() Examples

def find_angle_graph(velocity, vertical_velocity, interp=False):
    angle = []

    for i in range(len(velocity)):
        if velocity[i] == 0:
            angle.append(angle[-1])
        else:
            ratio = max(-1, min(vertical_velocity[i] / velocity[i], 1))
            angle.append(asin(ratio))

    angle = savgol_filter(angle, 5, 1)

    if interp:
        angle = savgol_filter(angle, 11, 1)
        return ss.medfilt(angle, kernel_size=7)

    return angle 

def unmeasure_np(self, hparams, x_measured_val, theta_val):
        if hparams.unmeasure_type == 'medfilt':
            unmeasure_func = lambda image, mask: signal.medfilt(image)
        elif hparams.unmeasure_type == 'inpaint-telea':
            inpaint_type = cv2.INPAINT_TELEA
            unmeasure_func = measure_utils.get_inpaint_func_opencv(hparams, inpaint_type)
        elif hparams.unmeasure_type == 'inpaint-ns':
            inpaint_type = cv2.INPAINT_NS
            unmeasure_func = measure_utils.get_inpaint_func_opencv(hparams, inpaint_type)
        elif hparams.unmeasure_type == 'inpaint-tv':
            unmeasure_func = measure_utils.get_inpaint_func_tv()
        elif hparams.unmeasure_type == 'blur':
            unmeasure_func = measure_utils.get_blur_func()
        else:
            raise NotImplementedError

        x_unmeasured_val = np.zeros_like(x_measured_val)
        for i in range(x_measured_val.shape[0]):
            x_unmeasured_val[i] = unmeasure_func(x_measured_val[i], theta_val[i])

        return x_unmeasured_val 

def circular_filter_1d(signal, window_size, kernel='gaussian'):

    """ This function filters circularly the signal inputted with a median filter of inputted size, in this context
    circularly means that the signal is wrapped around and then filtered
    inputs :
        - signal : 1D numpy array
        - window_size : size of the kernel, an int
    outputs :
        - signal_smoothed : 1D numpy array, same size as signal"""

    signal_extended = np.concatenate((signal, signal, signal))  # replicate signal at both ends
    if kernel == 'gaussian':
        signal_extended_smooth = ndimage.gaussian_filter(signal_extended, window_size)  # gaussian
    elif kernel == 'median':
        signal_extended_smooth = medfilt(signal_extended, window_size)  # median filtering
    else:
        raise Exception("Unknow type of kernel")

    signal_smoothed = signal_extended_smooth[len(signal):2*len(signal)]  # truncate back the signal

    return signal_smoothed 

def test_basic(self):
        f = [[50, 50, 50, 50, 50, 92, 18, 27, 65, 46],
             [50, 50, 50, 50, 50, 0, 72, 77, 68, 66],
             [50, 50, 50, 50, 50, 46, 47, 19, 64, 77],
             [50, 50, 50, 50, 50, 42, 15, 29, 95, 35],
             [50, 50, 50, 50, 50, 46, 34, 9, 21, 66],
             [70, 97, 28, 68, 78, 77, 61, 58, 71, 42],
             [64, 53, 44, 29, 68, 32, 19, 68, 24, 84],
             [3, 33, 53, 67, 1, 78, 74, 55, 12, 83],
             [7, 11, 46, 70, 60, 47, 24, 43, 61, 26],
             [32, 61, 88, 7, 39, 4, 92, 64, 45, 61]]

        d = signal.medfilt(f, [7, 3])
        e = signal.medfilt2d(np.array(f, np.float), [7, 3])
        assert_array_equal(d, [[0, 50, 50, 50, 42, 15, 15, 18, 27, 0],
                               [0, 50, 50, 50, 50, 42, 19, 21, 29, 0],
                               [50, 50, 50, 50, 50, 47, 34, 34, 46, 35],
                               [50, 50, 50, 50, 50, 50, 42, 47, 64, 42],
                               [50, 50, 50, 50, 50, 50, 46, 55, 64, 35],
                               [33, 50, 50, 50, 50, 47, 46, 43, 55, 26],
                               [32, 50, 50, 50, 50, 47, 46, 45, 55, 26],
                               [7, 46, 50, 50, 47, 46, 46, 43, 45, 21],
                               [0, 32, 33, 39, 32, 32, 43, 43, 43, 0],
                               [0, 7, 11, 7, 4, 4, 19, 19, 24, 0]])
        assert_array_equal(d, e) 

def compute_medfilt(arr):
    """Median filtered signal as features.

    Parameters
    ----------
    arr : ndarray, shape (n_channels, n_times)

    Returns
    -------
    output : (n_channels * n_times,)
    """
    return medfilt(arr, kernel_size=(1, 5)).ravel()


###############################################################################
# Prepare for the classification task
# -----------------------------------
#
# In addition to the new feature function, we also propose to extract the
# mean of the data: 

def transform(self,X):
        """Detrend each column of X

        Parameters
        ----------
        X : DataFrame in `traces` structure [n_samples, n_traces]

        Returns
        -------
        Xt : DataFrame in `traces` structure [n_samples, n_traces]
            The detrended data.
        """
        self.fit_params = {}
        X_new = X.copy()
        for col in X.columns:
            tmp_data = X[col].values.astype(np.double)
            mf = medfilt(tmp_data, self.window)
            mf = np.minimum(mf, self.peak_std_threshold * self._robust_std(mf))
            self.fit_params[col] = dict(mf=mf)
            X_new[col] = tmp_data - mf

        return X_new 

def loadData(log_folder, smooth, bin_size, is_es=False):
    """
    :param log_folder: (str)
    :param smooth: (int) Smoothing method
    :param bin_size: (int)
    :param is_es: (bool)
    :return:
    """
    result, timesteps = loadCsv(log_folder, is_es=is_es)

    if len(result) < bin_size:
        return [None, None]

    x, y = np.array(result)[:, 0], np.array(result)[:, 1]

    if smooth == 1:
        x, y = smoothRewardCurve(x, y)

    if smooth == 2:
        y = medfilt(y, kernel_size=9)

    x, y = fixPoint(x, y, bin_size)
    return [x, y] 

def test_basic(self):
        f = [[50, 50, 50, 50, 50, 92, 18, 27, 65, 46],
             [50, 50, 50, 50, 50, 0, 72, 77, 68, 66],
             [50, 50, 50, 50, 50, 46, 47, 19, 64, 77],
             [50, 50, 50, 50, 50, 42, 15, 29, 95, 35],
             [50, 50, 50, 50, 50, 46, 34, 9, 21, 66],
             [70, 97, 28, 68, 78, 77, 61, 58, 71, 42],
             [64, 53, 44, 29, 68, 32, 19, 68, 24, 84],
             [3, 33, 53, 67, 1, 78, 74, 55, 12, 83],
             [7, 11, 46, 70, 60, 47, 24, 43, 61, 26],
             [32, 61, 88, 7, 39, 4, 92, 64, 45, 61]]

        d = signal.medfilt(f, [7, 3])
        e = signal.medfilt2d(np.array(f, float), [7, 3])
        assert_array_equal(d, [[0, 50, 50, 50, 42, 15, 15, 18, 27, 0],
                               [0, 50, 50, 50, 50, 42, 19, 21, 29, 0],
                               [50, 50, 50, 50, 50, 47, 34, 34, 46, 35],
                               [50, 50, 50, 50, 50, 50, 42, 47, 64, 42],
                               [50, 50, 50, 50, 50, 50, 46, 55, 64, 35],
                               [33, 50, 50, 50, 50, 47, 46, 43, 55, 26],
                               [32, 50, 50, 50, 50, 47, 46, 45, 55, 26],
                               [7, 46, 50, 50, 47, 46, 46, 43, 45, 21],
                               [0, 32, 33, 39, 32, 32, 43, 43, 43, 0],
                               [0, 7, 11, 7, 4, 4, 19, 19, 24, 0]])
        assert_array_equal(d, e) 

def smooth_bbox_params(bbox_params, kernel_size=11, sigma=8):
    """
    Applies median filtering and then gaussian filtering to bounding box
    parameters.

    Args:
        bbox_params (Nx3): [cx, cy, scale].
        kernel_size (int): Kernel size for median filtering (must be odd).
        sigma (float): Sigma for gaussian smoothing.

    Returns:
        Smoothed bounding box parameters (Nx3).
    """
    smoothed = np.array([signal.medfilt(param, kernel_size)
                         for param in bbox_params.T]).T
    return np.array([gaussian_filter1d(traj, sigma) for traj in smoothed.T]).T 

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 signal.medfilt(x, kernel_size) 

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 MedianFilter(x, kernel_size=5):
    '''
    A silly wrapper around :py:func:`scipy.signal.medfilt`.

    '''

    if kernel_size % 2 == 0:
        kernel_size += 1
    return medfilt(x, kernel_size=kernel_size) 

def apply_median_filter(self, window_size=3):
        if self.mean_signal is None:
            self.x = signal.medfilt(self.x, window_size)
            self.y = signal.medfilt(self.y, window_size)
            self.z = signal.medfilt(self.z, window_size)
        else:
            self.mean_signal = signal.medfilt(self.mean_signal, window_size) 

def interpolate(self):
        pitch = np.zeros((self.nframes))
        pitch[:] = self.samp_values
        pitch2 = medfilt(self.samp_values, self.SMOOTH_FACTOR)

        # This part in the original code is kind of confused and caused
        # some problems with the extrapolated points before the first
        # voiced frame and after the last voiced frame. So, I made some
        # small modifications in order to make it work better.
        edges = self.edges_finder(pitch)
        first_sample = pitch[0]
        last_sample = pitch[-1]

        if len(np.nonzero(pitch2)[0]) < 2:
            pitch[pitch == 0] = self.PTCH_TYP
        else:
            nz_pitch = pitch2[pitch2 > 0]
            pitch2 = scipy_interp.pchip(np.nonzero(pitch2)[0],
                                        nz_pitch)(range(self.nframes))
            pitch[pitch == 0] = pitch2[pitch == 0]
        if self.SMOOTH > 0:
            pitch = medfilt(pitch, self.SMOOTH_FACTOR)
        try:
            if first_sample == 0:
                pitch[:edges[0]-1] = pitch[edges[0]]
            if last_sample == 0:
                pitch[edges[-1]+1:] = pitch[edges[-1]]
        except:
            pass
        self.samp_interp = pitch 

def denoise(volume, kernel_size=3):
    return medfilt(volume, kernel_size) 

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 on_epoch_end(self, callback_data, model, epoch):
        # convert to numpy arrays
        data_batch = model.data_batch.get()
        noise_batch = model.noise_batch.get()
        # value transform
        data_batch = self._value_transform(data_batch)
        noise_batch = self._value_transform(noise_batch)
        # shape transform
        data_canvas = self._shape_transform(data_batch)
        noise_canvas = self._shape_transform(noise_batch)
        # plotting options
        im_args = dict(interpolation="nearest", vmin=0., vmax=1.)
        if self.nchan == 1:
            im_args['cmap'] = plt.get_cmap("gray")
        fname = self.filename+'_data_'+'{:03d}'.format(epoch)+'.png'
        Image.fromarray(np.uint8(data_canvas*255)).convert('RGB').save(fname)
        fname = self.filename+'_noise_'+'{:03d}'.format(epoch)+'.png'
        Image.fromarray(np.uint8(noise_canvas*255)).convert('RGB').save(fname)

        # plot logged WGAN costs if logged
        if model.cost.costfunc.func == 'wasserstein':
            giter = callback_data['gan/gen_iter'][:]
            nonzeros = np.where(giter)
            giter = giter[nonzeros]
            cost_dis = callback_data['gan/cost_dis'][:][nonzeros]
            w_dist = medfilt(np.array(-cost_dis, dtype='float64'), kernel_size=101)
            plt.figure(figsize=(400/self.dpi, 300/self.dpi), dpi=self.dpi)
            plt.plot(giter, -cost_dis, 'k-', lw=0.25)
            plt.plot(giter, w_dist, 'r-', lw=2.)
            plt.title(self.filename, fontsize=self.font_size)
            plt.xlabel("Generator Iterations", fontsize=self.font_size)
            plt.ylabel("Wasserstein estimate", fontsize=self.font_size)
            plt.margins(0, 0, tight=True)
            plt.savefig(self.filename+'_training.png', bbox_inches='tight')
            plt.close() 

def medianFilter(self, kernel_size=5):
        """
        **SUMMARY**

        Apply median filter on the data

        **PARAMETERS**

        * *kernel_size* - Size of the filter (should be odd int) - int

        **RETURNS**

        A LineScan object with the median filter applied to the values.

        **EXAMPLE**

        >>> ls = img.getLineScan(x=10)
        >>> mf = ls.medianFilter()
        >>> plt.plot(ls)
        >>> plt.plot(mf)
        """
        try:
            from scipy.signal import medfilt
        except ImportError:
            warnings.warn("Scipy vesion >= 0.11 requierd.")
            return None
        if kernel_size % 2 == 0:
            kernel_size-=1
            print "Kernel Size should be odd. New kernel size =" , (kernel_size)
        
        medfilt_array = medfilt(np.asarray(self[:]), kernel_size)
        retVal = LineScan(medfilt_array.astype("uint8").tolist(), image=self.image,pointLoc=self.pointLoc,pt1=self.pt1,pt2=self.pt2, x=self.col, y=self.row)
        retVal._update(self)
        return retVal 

def detect_with_thresholding(metric,
                             thrh_type,
                             thrh_value,
                             proc_type,
                             proc_value,
                             debug_file=None):

    assert (thrh_type in ['max', 'mean'])
    assert (proc_type in ['dilation', 'median'])

    out_detections = []

    if thrh_type == 'max':
        mask = metric > thrh_value

    elif thrh_type == 'mean':
        mask = metric > (thrh_value * metric.mean())

    if proc_type == 'dilation':
        mask = dilation(mask, np.array([[1] for _ in range(proc_value)]))
    elif proc_type == 'median':
        mask = medfilt(mask[:, 0], kernel_size=proc_value)
        # kernel_size should be odd
        mask = np.expand_dims(mask, axis=1)

    return mask


################ Output Detection To Files ################ 

def median_filter(self) -> None:
        """
        Applies a median filter on the internal 1d signal. Can be useful for cosmic ray correction
        after temporal de-noising

        Returns
        -------
        NoneType
            None
        """

        self._m_data = medfilt(self._m_data, 19) 

def load_exp_sequence(path, use_medfilt=False, ksize=3):
    exp = np.load(path).astype(np.float32)
    if use_medfilt:
        exp = medfilt(exp, kernel_size=(ksize,1)).astype(np.float32)
    return exp


# estimate audio sequence -> exp 

def unmeasure_np(self, hparams, x_measured_val, theta_val):
        if hparams.unmeasure_type == 'medfilt':
            unmeasure_func = lambda image, mask: signal.medfilt(image)
        # elif hparams.unmeasure_type == 'inpaint-telea':
        #     inpaint_type = cv2.INPAINT_TELEA
        #     unmeasure_func = amb_measure_utils.get_inpaint_func_opencv(inpaint_type)
        # elif hparams.unmeasure_type == 'inpaint-ns':
        #     inpaint_type = cv2.INPAINT_NS
        #     unmeasure_func = amb_measure_utils.get_inpaint_func_opencv(inpaint_type)
        # elif hparams.unmeasure_type == 'inpaint-tv':
        #     assert hparams.dataset == 'mnist'  # Single channel support only
        #     unmeasure_func = amb_measure_utils.get_inpaint_func_tv()
        elif hparams.unmeasure_type == 'blur':
            # TODO(abora): Move radius and size to hparams
            gaussian_filter = amb_measure_utils.get_gaussian_filter(radius=1, size=5)
            def unmeasure_func(image, mask):
                blurred = np.zeros_like(image)
                for c in range(image.shape[2]):
                    blurred[:, :, c] = signal.convolve2d(image[:, :, c], gaussian_filter, mode='same')
                return blurred
        else:
            raise NotImplementedError

        x_unmeasured_val = np.zeros_like(x_measured_val)
        for i in range(x_measured_val.shape[0]):
            x_unmeasured_val[i] = unmeasure_func(x_measured_val[i], theta_val[i])
        return x_unmeasured_val 

def smooth_magseries_signal_medfilt(mags, windowsize):
    '''This smooths the magseries with a simple median filter.

    This function pads with zeros near the boundary, see:

    https://stackoverflow.com/questions/24585706/scipy-medfilt-wrong-result

    Typically this is bad.

    Parameters
    ----------

    mags : np.array
        The input mags/flux time-series to smooth.

    windowsize : int
        This is a odd integer containing the smoothing window size.

    Returns
    -------

    np.array
        The smoothed mag/flux time-series array.

    '''

    return medfilt(mags, windowsize) 

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 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 smooth_bbox_params(bbox_params, kernel_size=11, sigma=8):
    """
    Applies median filtering and then gaussian filtering to bounding box
    parameters.
    Args:
        bbox_params (Nx3): [cx, cy, scale].
        kernel_size (int): Kernel size for median filtering (must be odd).
        sigma (float): Sigma for gaussian smoothing.
    Returns:
        Smoothed bounding box parameters (Nx3).
    """
    smoothed = np.array([signal.medfilt(param, kernel_size)
                         for param in bbox_params.T]).T
    return np.array([gaussian_filter1d(traj, sigma) for traj in smoothed.T]).T 

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 test_none(self):
        # Ticket #1124. Ensure this does not segfault.
        try:
            signal.medfilt(None)
        except:
            pass
        # Expand on this test to avoid a regression with possible contiguous
        # numpy arrays that have odd strides. The stride value below gets
        # us into wrong memory if used (but it does not need to be used)
        dummy = np.arange(10, dtype=np.float64)
        a = dummy[5:6]
        a.strides = 16
        assert_(signal.medfilt(a, 1) == 5.) 

def test_refcounting(self):
        # Check a refcounting-related crash
        a = Decimal(123)
        x = np.array([a, a], dtype=object)
        if hasattr(sys, 'getrefcount'):
            n = 2 * sys.getrefcount(a)
        else:
            n = 10
        # Shouldn't segfault:
        for j in range(n):
            signal.medfilt(x)
        if hasattr(sys, 'getrefcount'):
            assert_(sys.getrefcount(a) < n)
        assert_equal(x, [a, a]) 

def test_none(self):
        # Ticket #1124. Ensure this does not segfault.
        try:
            signal.medfilt(None)
        except:
            pass
        # Expand on this test to avoid a regression with possible contiguous
        # numpy arrays that have odd strides. The stride value below gets
        # us into wrong memory if used (but it does not need to be used)
        dummy = np.arange(10, dtype=np.float64)
        a = dummy[5:6]
        a.strides = 16
        assert_(signal.medfilt(a, 1) == 5.)