def test_cuda():
    """Test CUDA-based filtering
    """
    Fs = 500
    sig_len_secs = 20
    a = np.random.randn(sig_len_secs * Fs)

    set_log_file(log_file, overwrite=True)
    for fl in [None, 2048]:
        bp = band_pass_filter(a, Fs, 4, 8, n_jobs=1, filter_length=fl)
        bs = band_stop_filter(a, Fs, 4 - 0.5, 8 + 0.5, n_jobs=1,
                              filter_length=fl)
        lp = low_pass_filter(a, Fs, 8, n_jobs=1, filter_length=fl)
        hp = high_pass_filter(lp, Fs, 4, n_jobs=1, filter_length=fl)

        bp_c = band_pass_filter(a, Fs, 4, 8, n_jobs='cuda', filter_length=fl,
                                verbose='INFO')
        bs_c = band_stop_filter(a, Fs, 4 - 0.5, 8 + 0.5, n_jobs='cuda',
                                filter_length=fl, verbose='INFO')
        lp_c = low_pass_filter(a, Fs, 8, n_jobs='cuda', filter_length=fl,
                               verbose='INFO')
        hp_c = high_pass_filter(lp, Fs, 4, n_jobs='cuda', filter_length=fl,
                                verbose='INFO')

        assert_array_almost_equal(bp, bp_c, 12)
        assert_array_almost_equal(bs, bs_c, 12)
        assert_array_almost_equal(lp, lp_c, 12)
        assert_array_almost_equal(hp, hp_c, 12)

    # check to make sure we actually used CUDA
    set_log_file()
    out = open(log_file).readlines()
    assert_true(sum(['Using CUDA for FFT FIR filtering' in o
                     for o in out]) == 8)

def test_filters():
    Fs = 500
    sig_len_secs = 60

    # Filtering of short signals (filter length = len(a))
    a = np.random.randn(sig_len_secs * Fs)
    bp = band_pass_filter(a, Fs, 4, 8)
    lp = low_pass_filter(a, Fs, 8)
    hp = high_pass_filter(lp, Fs, 4)
    assert_array_almost_equal(hp, bp, 2)

    # Overlap-add filtering with a fixed filter length
    filter_length = 8192
    bp_oa = band_pass_filter(a, Fs, 4, 8, filter_length)
    lp_oa = low_pass_filter(a, Fs, 8, filter_length)
    hp_oa = high_pass_filter(lp_oa, Fs, 4, filter_length)
    assert_array_almost_equal(hp_oa, bp_oa, 2)

    # The two methods should give the same result
    # As filtering for short signals uses a circular convolution (FFT) and
    # the overlap-add filter implements a linear convolution, the signal
    # boundary will be slightly different and we ignore it
    n_edge_ignore = 1000
    assert_array_almost_equal(hp[n_edge_ignore:-n_edge_ignore],
                              hp_oa[n_edge_ignore:-n_edge_ignore], 2)

def test_cuda():
    """Test CUDA-based filtering"""
    # NOTE: don't make test_cuda() the last test, or pycuda might spew
    # some warnings about clean-up failing
    # Also, using `n_jobs='cuda'` on a non-CUDA system should be fine,
    # as it should fall back to using n_jobs=1.
    sfreq = 500
    sig_len_secs = 20
    a = rng.randn(sig_len_secs * sfreq)

    with catch_logging() as log_file:
        for fl in ['auto', '10s', 2048]:
            bp = band_pass_filter(a, sfreq, 4, 8, fl, 1.0, 1.0, n_jobs=1,
                                  phase='zero')
            bs = band_stop_filter(a, sfreq, 4 - 1.0, 8 + 1.0, fl, 1.0, 1.0,
                                  n_jobs=1, phase='zero')
            lp = low_pass_filter(a, sfreq, 8, fl, 1.0, n_jobs=1, phase='zero')
            hp = high_pass_filter(lp, sfreq, 4, fl, 1.0, n_jobs=1,
                                  phase='zero')

            bp_c = band_pass_filter(a, sfreq, 4, 8, fl, 1.0, 1.0,
                                    n_jobs='cuda', verbose='INFO',
                                    phase='zero')
            bs_c = band_stop_filter(a, sfreq, 4 - 1.0, 8 + 1.0, fl, 1.0, 1.0,
                                    n_jobs='cuda', verbose='INFO',
                                    phase='zero')
            lp_c = low_pass_filter(a, sfreq, 8, fl, 1.0,
                                   n_jobs='cuda', verbose='INFO',
                                   phase='zero')
            hp_c = high_pass_filter(lp, sfreq, 4, fl, 1.0,
                                    n_jobs='cuda', verbose='INFO',
                                    phase='zero')

            assert_array_almost_equal(bp, bp_c, 12)
            assert_array_almost_equal(bs, bs_c, 12)
            assert_array_almost_equal(lp, lp_c, 12)
            assert_array_almost_equal(hp, hp_c, 12)

    # check to make sure we actually used CUDA
    out = log_file.getvalue().split('\n')[:-1]
    # triage based on whether or not we actually expected to use CUDA
    from mne.cuda import _cuda_capable  # allow above funs to set it
    tot = 12 if _cuda_capable else 0
    assert_true(sum(['Using CUDA for FFT FIR filtering' in o
                     for o in out]) == tot)

    # check resampling
    for window in ('boxcar', 'triang'):
        for N in (997, 1000):  # one prime, one even
            a = rng.randn(2, N)
            for fro, to in ((1, 2), (2, 1), (1, 3), (3, 1)):
                a1 = resample(a, fro, to, n_jobs=1, npad='auto',
                              window=window)
                a2 = resample(a, fro, to, n_jobs='cuda', npad='auto',
                              window=window)
                assert_allclose(a1, a2, rtol=1e-7, atol=1e-14)
    assert_array_almost_equal(a1, a2, 14)
    assert_array_equal(resample([0, 0], 2, 1, n_jobs='cuda'), [0., 0., 0., 0.])
    assert_array_equal(resample(np.zeros(2, np.float32), 2, 1, n_jobs='cuda'),
                       [0., 0., 0., 0.])

def test_cuda():
    """Test CUDA-based filtering
    """
    # NOTE: don't make test_cuda() the last test, or pycuda might spew
    # some warnings about clean-up failing
    # Also, using `n_jobs='cuda'` on a non-CUDA system should be fine,
    # as it should fall back to using n_jobs=1.
    tempdir = _TempDir()
    log_file = op.join(tempdir, 'temp_log.txt')
    sfreq = 500
    sig_len_secs = 20
    a = np.random.randn(sig_len_secs * sfreq)

    set_log_file(log_file, overwrite=True)
    for fl in ['10s', None, 2048]:
        bp = band_pass_filter(a, sfreq, 4, 8, n_jobs=1, filter_length=fl)
        bs = band_stop_filter(a, sfreq, 4 - 0.5, 8 + 0.5, n_jobs=1,
                              filter_length=fl)
        lp = low_pass_filter(a, sfreq, 8, n_jobs=1, filter_length=fl)
        hp = high_pass_filter(lp, sfreq, 4, n_jobs=1, filter_length=fl)

        bp_c = band_pass_filter(a, sfreq, 4, 8, n_jobs='cuda',
                                filter_length=fl, verbose='INFO')
        bs_c = band_stop_filter(a, sfreq, 4 - 0.5, 8 + 0.5, n_jobs='cuda',
                                filter_length=fl, verbose='INFO')
        lp_c = low_pass_filter(a, sfreq, 8, n_jobs='cuda', filter_length=fl,
                               verbose='INFO')
        hp_c = high_pass_filter(lp, sfreq, 4, n_jobs='cuda', filter_length=fl,
                                verbose='INFO')

        assert_array_almost_equal(bp, bp_c, 12)
        assert_array_almost_equal(bs, bs_c, 12)
        assert_array_almost_equal(lp, lp_c, 12)
        assert_array_almost_equal(hp, hp_c, 12)

    # check to make sure we actually used CUDA
    set_log_file()
    with open(log_file) as fid:
        out = fid.readlines()
    # triage based on whether or not we actually expected to use CUDA
    from mne.cuda import _cuda_capable  # allow above funs to set it
    tot = 12 if _cuda_capable else 0
    assert_true(sum(['Using CUDA for FFT FIR filtering' in o
                     for o in out]) == tot)

    # check resampling
    a = np.random.RandomState(0).randn(3, sig_len_secs * sfreq)
    a1 = resample(a, 1, 2, n_jobs=2, npad=0)
    a2 = resample(a, 1, 2, n_jobs='cuda', npad=0)
    a3 = resample(a, 2, 1, n_jobs=2, npad=0)
    a4 = resample(a, 2, 1, n_jobs='cuda', npad=0)
    assert_array_almost_equal(a3, a4, 14)
    assert_array_almost_equal(a1, a2, 14)
    assert_array_equal(resample([0, 0], 2, 1, n_jobs='cuda'), [0., 0., 0., 0.])
    assert_array_equal(resample(np.zeros(2, np.float32), 2, 1, n_jobs='cuda'),
                       [0., 0., 0., 0.])

def test_iir_stability():
    """Test IIR filter stability check."""
    sig = np.empty(1000)
    sfreq = 1000
    # This will make an unstable filter, should throw RuntimeError
    assert_raises(RuntimeError, high_pass_filter, sig, sfreq, 0.6,
                  method='iir', iir_params=dict(ftype='butter', order=8,
                                                output='ba'))
    # This one should work just fine
    high_pass_filter(sig, sfreq, 0.6, method='iir',
                     iir_params=dict(ftype='butter', order=8, output='sos'))
    # bad system type
    assert_raises(ValueError, high_pass_filter, sig, sfreq, 0.6, method='iir',
                  iir_params=dict(ftype='butter', order=8, output='foo'))
    # missing ftype
    assert_raises(RuntimeError, high_pass_filter, sig, sfreq, 0.6,
                  method='iir', iir_params=dict(order=8, output='sos'))
    # bad ftype
    assert_raises(RuntimeError, high_pass_filter, sig, sfreq, 0.6,
                  method='iir',
                  iir_params=dict(order=8, ftype='foo', output='sos'))
    # missing gstop
    assert_raises(RuntimeError, high_pass_filter, sig, sfreq, 0.6,
                  method='iir', iir_params=dict(gpass=0.5, output='sos'))
    # can't pass iir_params if method='fft'
    assert_raises(ValueError, high_pass_filter, sig, sfreq, 0.1,
                  method='fft', iir_params=dict(ftype='butter', order=2,
                                                output='sos'))
    # method must be string
    assert_raises(TypeError, high_pass_filter, sig, sfreq, 0.1,
                  method=1)
    # unknown method
    assert_raises(ValueError, high_pass_filter, sig, sfreq, 0.1,
                  method='blah')
    # bad iir_params
    assert_raises(TypeError, high_pass_filter, sig, sfreq, 0.1,
                  method='iir', iir_params='blah')
    assert_raises(ValueError, high_pass_filter, sig, sfreq, 0.1,
                  method='fft', iir_params=dict())

    # should pass because dafault trans_bandwidth is not relevant
    iir_params = dict(ftype='butter', order=2, output='sos')
    x_sos = high_pass_filter(sig, 250, 0.5, method='iir',
                             iir_params=iir_params)
    iir_params_sos = construct_iir_filter(iir_params, f_pass=0.5, sfreq=250,
                                          btype='highpass')
    x_sos_2 = high_pass_filter(sig, 250, 0.5, method='iir',
                               iir_params=iir_params_sos)
    assert_allclose(x_sos[100:-100], x_sos_2[100:-100])
    x_ba = high_pass_filter(sig, 250, 0.5, method='iir',
                            iir_params=dict(ftype='butter', order=2,
                                            output='ba'))
    # Note that this will fail for higher orders (e.g., 6) showing the
    # hopefully decreased numerical error of SOS
    assert_allclose(x_sos[100:-100], x_ba[100:-100])

def test_cuda():
    """Test CUDA-based filtering
    """
    # NOTE: don't make test_cuda() the last test, or pycuda might spew
    # some warnings about clean-up failing
    # Also, using `n_jobs='cuda'` on a non-CUDA system should be fine,
    # as it should fall back to using n_jobs=1.
    sfreq = 500
    sig_len_secs = 20
    a = rng.randn(sig_len_secs * sfreq)

    with catch_logging() as log_file:
        for fl in ['10s', None, 2048]:
            bp = band_pass_filter(a, sfreq, 4, 8, n_jobs=1, filter_length=fl)
            bs = band_stop_filter(a, sfreq, 4 - 0.5, 8 + 0.5, n_jobs=1,
                                  filter_length=fl)
            lp = low_pass_filter(a, sfreq, 8, n_jobs=1, filter_length=fl)
            hp = high_pass_filter(lp, sfreq, 4, n_jobs=1, filter_length=fl)

            bp_c = band_pass_filter(a, sfreq, 4, 8, n_jobs='cuda',
                                    filter_length=fl, verbose='INFO')
            bs_c = band_stop_filter(a, sfreq, 4 - 0.5, 8 + 0.5, n_jobs='cuda',
                                    filter_length=fl, verbose='INFO')
            lp_c = low_pass_filter(a, sfreq, 8, n_jobs='cuda',
                                   filter_length=fl, verbose='INFO')
            hp_c = high_pass_filter(lp, sfreq, 4, n_jobs='cuda',
                                    filter_length=fl, verbose='INFO')

            assert_array_almost_equal(bp, bp_c, 12)
            assert_array_almost_equal(bs, bs_c, 12)
            assert_array_almost_equal(lp, lp_c, 12)
            assert_array_almost_equal(hp, hp_c, 12)

    # check to make sure we actually used CUDA
    out = log_file.getvalue().split('\n')[:-1]
    # triage based on whether or not we actually expected to use CUDA
    from mne.cuda import _cuda_capable  # allow above funs to set it
    tot = 12 if _cuda_capable else 0
    assert_true(sum(['Using CUDA for FFT FIR filtering' in o
                     for o in out]) == tot)

    # check resampling
    a = rng.randn(3, sig_len_secs * sfreq)
    a1 = resample(a, 1, 2, n_jobs=2, npad=0)
    a2 = resample(a, 1, 2, n_jobs='cuda', npad=0)
    a3 = resample(a, 2, 1, n_jobs=2, npad=0)
    a4 = resample(a, 2, 1, n_jobs='cuda', npad=0)
    assert_array_almost_equal(a3, a4, 14)
    assert_array_almost_equal(a1, a2, 14)
    assert_array_equal(resample([0, 0], 2, 1, n_jobs='cuda'), [0., 0., 0., 0.])
    assert_array_equal(resample(np.zeros(2, np.float32), 2, 1, n_jobs='cuda'),
                       [0., 0., 0., 0.])

def test_cuda():
    """Test CUDA-based filtering
    """
    # NOTE: don't make test_cuda() the last test, or pycuda might spew
    # some warnings about clean-up failing
    Fs = 500
    sig_len_secs = 20
    a = np.random.randn(sig_len_secs * Fs)

    set_log_file(log_file, overwrite=True)
    for fl in ['10s', None, 2048]:
        bp = band_pass_filter(a, Fs, 4, 8, n_jobs=1, filter_length=fl)
        bs = band_stop_filter(a, Fs, 4 - 0.5, 8 + 0.5, n_jobs=1,
                              filter_length=fl)
        lp = low_pass_filter(a, Fs, 8, n_jobs=1, filter_length=fl)
        hp = high_pass_filter(lp, Fs, 4, n_jobs=1, filter_length=fl)

        bp_c = band_pass_filter(a, Fs, 4, 8, n_jobs='cuda', filter_length=fl,
                                verbose='INFO')
        bs_c = band_stop_filter(a, Fs, 4 - 0.5, 8 + 0.5, n_jobs='cuda',
                                filter_length=fl, verbose='INFO')
        lp_c = low_pass_filter(a, Fs, 8, n_jobs='cuda', filter_length=fl,
                               verbose='INFO')
        hp_c = high_pass_filter(lp, Fs, 4, n_jobs='cuda', filter_length=fl,
                                verbose='INFO')

        assert_array_almost_equal(bp, bp_c, 12)
        assert_array_almost_equal(bs, bs_c, 12)
        assert_array_almost_equal(lp, lp_c, 12)
        assert_array_almost_equal(hp, hp_c, 12)

    # check to make sure we actually used CUDA
    set_log_file()
    with open(log_file) as fid:
        out = fid.readlines()
    assert_true(sum(['Using CUDA for FFT FIR filtering' in o
                     for o in out]) == 12)

    # check resampling
    a = np.random.RandomState(0).randn(3, sig_len_secs * Fs)
    a1 = resample(a, 1, 2, n_jobs=2, npad=0)
    a2 = resample(a, 1, 2, n_jobs='cuda', npad=0)
    a3 = resample(a, 2, 1, n_jobs=2, npad=0)
    a4 = resample(a, 2, 1, n_jobs='cuda', npad=0)
    assert_array_almost_equal(a3, a4, 14)
    assert_array_almost_equal(a1, a2, 14)

def test_iir_stability():
    """Test IIR filter stability check
    """
    sig = np.empty(1000)
    sfreq = 1000
    # This will make an unstable filter, should throw RuntimeError
    assert_raises(
        RuntimeError, high_pass_filter, sig, sfreq, 0.6, method="iir", iir_params=dict(ftype="butter", order=8)
    )
    # can't pass iir_params if method='fir'
    assert_raises(ValueError, high_pass_filter, sig, sfreq, 0.1, method="fir", iir_params=dict(ftype="butter", order=2))
    # method must be string
    assert_raises(TypeError, high_pass_filter, sig, sfreq, 0.1, method=1)
    # unknown method
    assert_raises(ValueError, high_pass_filter, sig, sfreq, 0.1, method="blah")
    # bad iir_params
    assert_raises(ValueError, high_pass_filter, sig, sfreq, 0.1, method="fir", iir_params="blah")

    # should pass because dafault trans_bandwidth is not relevant
    high_pass_filter(sig, 250, 0.5, method="iir", iir_params=dict(ftype="butter", order=6))

 def process(self, data):
     if self.type == 'low-pass':
         return low_pass_filter(data, **self.params)
     elif self.type == 'high-pass':
         return high_pass_filter(data, **self.params)
     elif self.type == 'band-pass':
         return band_pass_filter(data, **self.params)
     elif self.type == 'band-stop':
         return band_stop_filter(data, **self.params)
     elif self.type == 'notch':
         return notch_filter(data, **self.params)
     else:
         raise ValueError('Unsupported filter type: {}'.format(self.type))

def test_filters():
    """Test low-, band-, and high-pass filters"""
    Fs = 500
    sig_len_secs = 60

    # Filtering of short signals (filter length = len(a))
    a = np.random.randn(sig_len_secs * Fs)
    bp = band_pass_filter(a, Fs, 4, 8)
    lp = low_pass_filter(a, Fs, 8)
    hp = high_pass_filter(lp, Fs, 4)
    assert_array_almost_equal(hp, bp, 2)

    # Overlap-add filtering with a fixed filter length
    filter_length = 8192
    bp_oa = band_pass_filter(a, Fs, 4, 8, filter_length)
    lp_oa = low_pass_filter(a, Fs, 8, filter_length)
    hp_oa = high_pass_filter(lp_oa, Fs, 4, filter_length)
    assert_array_almost_equal(hp_oa, bp_oa, 2)

    # The two methods should give the same result
    # As filtering for short signals uses a circular convolution (FFT) and
    # the overlap-add filter implements a linear convolution, the signal
    # boundary will be slightly different and we ignore it
    n_edge_ignore = 1000
    assert_array_almost_equal(hp[n_edge_ignore:-n_edge_ignore],
                              hp_oa[n_edge_ignore:-n_edge_ignore], 2)

    # and since these are low-passed, downsampling/upsampling should be close
    n_resamp_ignore = 10
    bp_up_dn = resample(resample(bp_oa, 2, 1), 1, 2)
    assert_array_almost_equal(bp_oa[n_resamp_ignore:-n_resamp_ignore],
                              bp_up_dn[n_resamp_ignore:-n_resamp_ignore], 2)
    # make sure we don't alias
    t = np.array(range(Fs*sig_len_secs))/float(Fs)
    # make sinusoid close to the Nyquist frequency
    sig = np.sin(2*np.pi*Fs/2.2*t)
    # signal should disappear with 2x downsampling
    sig_gone = resample(sig,1,2)[n_resamp_ignore:-n_resamp_ignore]
    assert_array_almost_equal(np.zeros_like(sig_gone), sig_gone, 2)

    def process(self, data):
        # fix for new MNE requirements
        import numpy as np
        data = np.asarray(data, dtype=np.float64)

        if self.type == 'low-pass':
            return low_pass_filter(data, **self.params)
        elif self.type == 'high-pass':
            return high_pass_filter(data, **self.params)
        elif self.type == 'band-pass':
            return band_pass_filter(data, **self.params)
        elif self.type == 'band-stop':
            return band_stop_filter(data, **self.params)
        elif self.type == 'notch':
            return notch_filter(data, **self.params)
        else:
            raise ValueError('Unsupported filter type: {}'.format(self.type))

 def apply_filters(self):
   Fs = int(self.frequency)
   Nf = int(Fs/2)
   lr = np.array([(x*60)-2 for x in range(1, Nf/60+1)], dtype=np.float64)
   hr = np.array([(x*60)+2 for x in range(1, Nf/60+1)], dtype=np.float64)
   for e in xrange(self.data.shape[0]):
     if np.all(self.data[e][:1000] == 0):
       continue
     self.data[e] = high_pass_filter(x = self.data[e], Fs = Fs, Fp = 0.4, 
                                     trans_bandwidth = 0.05, copy = False,
                                     filter_length=int(np.floor(Fs*32)), verbose=False)
     '''
     efft = self.fft['fft_e' + str(e)]
     line_ratio = (np.mean(efft[:,60])/np.mean(efft[:,50:58]) + 
                   np.mean(efft[:,60])/np.mean(efft[:,63:70])) / 2
     if line_ratio > 1.5:
     '''
     self.data[e] = band_stop_filter(x = self.data[e], Fs = Fs, Fp1 = lr, Fp2 = hr, 
                                       copy = False, verbose=False)

#.........这里部分代码省略.........
        units of "s" or "ms" (e.g., "10s" or "5500ms") will be converted
        to the shortest power-of-two length at least that duration.
        Not used for 'iir' filters.
    l_trans_bandwidth : float
        Width of the transition band at the low cut-off frequency in Hz
        (high pass or cutoff 1 in bandpass). Not used if 'order' is
        specified in iir_params.
    h_trans_bandwidth : float
        Width of the transition band at the high cut-off frequency in Hz
        (low pass or cutoff 2 in bandpass). Not used if 'order' is
        specified in iir_params.
    n_jobs : int | str
        Number of jobs to run in parallel. Can be 'cuda' if scikits.cuda
        is installed properly, CUDA is initialized, and method='fft'.
    method : str
        'fft' will use overlap-add FIR filtering, 'iir' will use IIR
        forward-backward filtering (via filtfilt).
    iir_params : dict | None
        Dictionary of parameters to use for IIR filtering.
        See mne.filter.construct_iir_filter for details. If iir_params
        is None and method="iir", 4th order Butterworth will be used.
    verbose : bool, str, int, or None
        If not None, override default verbose level (see mne.verbose).
        Defaults to raw.verbose.
    See Also
    --------
    mne.Epochs.savgol_filter
    """
    fname='ec_rest_before_tsss_mc_rsl.fif'
    raw = Raw(fname, preload=False)
    raw.preload_data() #  data becomes numpy.float64

    if verbose is None:
        verbose = raw.verbose
    fs = float(raw.info['sfreq'])
    if l_freq == 0:
        l_freq = None
    if h_freq is not None and h_freq > (fs / 2.):
        h_freq = None
    if l_freq is not None and not isinstance(l_freq, float):
        l_freq = float(l_freq)
    if h_freq is not None and not isinstance(h_freq, float):
        h_freq = float(h_freq)

    if not raw.preload:
        raise RuntimeError('Raw data needs to be preloaded to filter. Use '
                           'preload=True (or string) in the constructor.')
    if picks is None:
        if 'ICA ' in ','.join(raw.ch_names):
            pick_parameters = dict(misc=True, ref_meg=False)
        else:
            pick_parameters = dict(meg=True, eeg=True, ref_meg=False)
        picks = pick_types(raw.info, exclude=[], **pick_parameters)
        # let's be safe.
        if len(picks) < 1:
            raise RuntimeError('Could not find any valid channels for '
                               'your Raw object. Please contact the '
                               'MNE-Python developers.')

        # update info if filter is applied to all data channels,
        # and it's not a band-stop filter
        if h_freq is not None:
            if (l_freq is None or l_freq < h_freq) and \
               (raw.info["lowpass"] is None or
               h_freq < raw.info['lowpass']):
                    raw.info['lowpass'] = h_freq
        if l_freq is not None:
            if (h_freq is None or l_freq < h_freq) and \
               (raw.info["highpass"] is None or
               l_freq > raw.info['highpass']):
                    raw.info['highpass'] = l_freq
    if l_freq is None and h_freq is not None:
        low_pass_filter(raw._data, fs, h_freq,
                        filter_length=filter_length,
                        trans_bandwidth=h_trans_bandwidth, method=method,
                        iir_params=iir_params, picks=picks, n_jobs=n_jobs,
                        copy=False)
    if l_freq is not None and h_freq is None:
        high_pass_filter(raw._data, fs, l_freq,
                         filter_length=filter_length,
                         trans_bandwidth=l_trans_bandwidth, method=method,
                         iir_params=iir_params, picks=picks, n_jobs=n_jobs,
                         copy=False)
    if l_freq is not None and h_freq is not None:
        if l_freq < h_freq:
            raw._data = band_pass_filter(
                raw._data, fs, l_freq, h_freq,
                filter_length=filter_length,
                l_trans_bandwidth=l_trans_bandwidth,
                h_trans_bandwidth=h_trans_bandwidth,
                method=method, iir_params=iir_params, picks=picks,
                n_jobs=n_jobs, copy=False)
        else:
            raw._data = band_stop_filter(
                raw._data, fs, h_freq, l_freq,
                filter_length=filter_length,
                l_trans_bandwidth=h_trans_bandwidth,
                h_trans_bandwidth=l_trans_bandwidth, method=method,
                iir_params=iir_params, picks=picks, n_jobs=n_jobs,
                copy=False)

def test_filters():
    """Test low-, band-, high-pass, and band-stop filters plus resampling
    """
    sfreq = 500
    sig_len_secs = 30

    a = np.random.randn(2, sig_len_secs * sfreq)

    # let's test our catchers
    for fl in ['blah', [0, 1], 1000.5, '10ss', '10']:
        assert_raises(ValueError, band_pass_filter, a, sfreq, 4, 8,
                      filter_length=fl)
    for nj in ['blah', 0.5]:
        assert_raises(ValueError, band_pass_filter, a, sfreq, 4, 8, n_jobs=nj)
    # > Nyq/2
    assert_raises(ValueError, band_pass_filter, a, sfreq, 4, sfreq / 2.)
    assert_raises(ValueError, low_pass_filter, a, sfreq, sfreq / 2.)
    # check our short-filter warning:
    with warnings.catch_warnings(record=True) as w:
        # Warning for low attenuation
        band_pass_filter(a, sfreq, 1, 8, filter_length=1024)
        # Warning for too short a filter
        band_pass_filter(a, sfreq, 1, 8, filter_length='0.5s')
    assert_true(len(w) >= 2)

    # try new default and old default
    for fl in ['10s', '5000ms', None]:
        bp = band_pass_filter(a, sfreq, 4, 8, filter_length=fl)
        bs = band_stop_filter(a, sfreq, 4 - 0.5, 8 + 0.5, filter_length=fl)
        lp = low_pass_filter(a, sfreq, 8, filter_length=fl, n_jobs=2)
        hp = high_pass_filter(lp, sfreq, 4, filter_length=fl)
        assert_array_almost_equal(hp, bp, 2)
        assert_array_almost_equal(bp + bs, a, 1)

    # Overlap-add filtering with a fixed filter length
    filter_length = 8192
    bp_oa = band_pass_filter(a, sfreq, 4, 8, filter_length)
    bs_oa = band_stop_filter(a, sfreq, 4 - 0.5, 8 + 0.5, filter_length)
    lp_oa = low_pass_filter(a, sfreq, 8, filter_length)
    hp_oa = high_pass_filter(lp_oa, sfreq, 4, filter_length)
    assert_array_almost_equal(hp_oa, bp_oa, 2)
    # Our filters are no longer quite complementary with linear rolloffs :(
    # this is the tradeoff for stability of the filtering
    # obtained by directly using the result of firwin2 instead of
    # modifying it...
    assert_array_almost_equal(bp_oa + bs_oa, a, 1)

    # The two methods should give the same result
    # As filtering for short signals uses a circular convolution (FFT) and
    # the overlap-add filter implements a linear convolution, the signal
    # boundary will be slightly different and we ignore it
    n_edge_ignore = 0
    assert_array_almost_equal(hp[n_edge_ignore:-n_edge_ignore],
                              hp_oa[n_edge_ignore:-n_edge_ignore], 2)

    # and since these are low-passed, downsampling/upsampling should be close
    n_resamp_ignore = 10
    bp_up_dn = resample(resample(bp_oa, 2, 1, n_jobs=2), 1, 2, n_jobs=2)
    assert_array_almost_equal(bp_oa[n_resamp_ignore:-n_resamp_ignore],
                              bp_up_dn[n_resamp_ignore:-n_resamp_ignore], 2)
    # note that on systems without CUDA, this line serves as a test for a
    # graceful fallback to n_jobs=1
    bp_up_dn = resample(resample(bp_oa, 2, 1, n_jobs='cuda'), 1, 2,
                        n_jobs='cuda')
    assert_array_almost_equal(bp_oa[n_resamp_ignore:-n_resamp_ignore],
                              bp_up_dn[n_resamp_ignore:-n_resamp_ignore], 2)
    # test to make sure our resamling matches scipy's
    bp_up_dn = sp_resample(sp_resample(bp_oa, 2 * bp_oa.shape[-1], axis=-1,
                                       window='boxcar'),
                           bp_oa.shape[-1], window='boxcar', axis=-1)
    assert_array_almost_equal(bp_oa[n_resamp_ignore:-n_resamp_ignore],
                              bp_up_dn[n_resamp_ignore:-n_resamp_ignore], 2)

    # make sure we don't alias
    t = np.array(list(range(sfreq * sig_len_secs))) / float(sfreq)
    # make sinusoid close to the Nyquist frequency
    sig = np.sin(2 * np.pi * sfreq / 2.2 * t)
    # signal should disappear with 2x downsampling
    sig_gone = resample(sig, 1, 2)[n_resamp_ignore:-n_resamp_ignore]
    assert_array_almost_equal(np.zeros_like(sig_gone), sig_gone, 2)

    # let's construct some filters
    iir_params = dict(ftype='cheby1', gpass=1, gstop=20)
    iir_params = construct_iir_filter(iir_params, 40, 80, 1000, 'low')
    # this should be a third order filter
    assert_true(iir_params['a'].size - 1 == 3)
    assert_true(iir_params['b'].size - 1 == 3)
    iir_params = dict(ftype='butter', order=4)
    iir_params = construct_iir_filter(iir_params, 40, None, 1000, 'low')
    assert_true(iir_params['a'].size - 1 == 4)
    assert_true(iir_params['b'].size - 1 == 4)

    # check that picks work for 3d array with one channel and picks=[0]
    a = np.random.randn(5 * sfreq, 5 * sfreq)
    b = a[:, None, :]

    with warnings.catch_warnings(record=True) as w:
        a_filt = band_pass_filter(a, sfreq, 4, 8)
        b_filt = band_pass_filter(b, sfreq, 4, 8, picks=[0])

#.........这里部分代码省略.........

# Testing of the idea of "cleaning" the data before t-sneeing. Does not make a difference
#  Clean the data set to be fed into T-sne
# High pass over 2000Hz
num_of_spikes = ivm_data_filtered.shape[2]
ivm_data_double_filtered = np.zeros(shape=ivm_data_filtered.shape)
for spike in np.arange(0, num_of_spikes):
    ivm_data_double_filtered[:, :, spike] = filters.high_pass_filter(ivm_data_filtered[:, :, spike], Fs=sampling_freq,
                                                                     Fp=2000)
# Set to zero the channels that have to spike features
num_of_spikes = ivm_data_filtered.shape[2]
ivm_data_double_filtered_zeroed = np.zeros(shape=ivm_data_filtered.shape)
for spike in np.arange(0, num_of_spikes):
    for channel in np.arange(0, ivm_data_filtered.shape[0]):
        if np.std(ivm_data_filtered[channel, :, spike]) > 150:
            ivm_data_double_filtered_zeroed[channel, :, spike] = ivm_data_filtered[channel, :, spike]


spikes_to_include = 1500  # (about 3.8 minutes)
fin_time_point = all_cells_spike_triggers['9'][spikes_to_include] + 500

raw_data_file_ivm = os.path.join(data_folder, 'amplifier'+date+'T'+cell_capture_times['9']+'.bin')
raw_data_ivm = ephys.load_raw_data(raw_data_file_ivm, numchannels=num_ivm_channels, dtype=amp_dtype)
raw_data_ivm = raw_data_ivm.dataMatrix[:, :fin_time_point]

ivm_data_filtered = np.memmap(os.path.join(analysis_folder, types_of_data_to_load[data_to_load].format(good_cells)),
                                dtype=filtered_data_type,
                                mode='w+',
                                shape=shape_of_filt_spike_trig_ivm)
for spike in np.arange(0, num_of_spikes):
    trigger_point = spike_triggers[spike]
    start_point = int(trigger_point - (num_of_points_in_spike_trig_ivm + num_of_points_for_padding)/2)
    if start_point < 0:
        break
    end_point = int(trigger_point + (num_of_points_in_spike_trig_ivm + num_of_points_for_padding)/2)
    if end_point > raw_data_ivm.shape()[1]:
        break
    temp_unfiltered = raw_data_ivm.dataMatrix[:, start_point:end_point]
    temp_unfiltered = temp_unfiltered.astype(filtered_data_type)
    iir_params = {'order': 4, 'ftype': 'butter', 'padlen': 0}
    temp_filtered = filters.high_pass_filter(temp_unfiltered, sampling_freq, high_pass_freq, method='iir',
                                             iir_params=iir_params)  # 4th order Butter with no padding
    temp_filtered = temp_filtered[:, int(num_of_points_for_padding / 2):-int(num_of_points_for_padding / 2)]
    ivm_data_filtered[:, :, spike] = temp_filtered


# Load the already saved ivm_data_filtered
shape_of_filt_spike_trig_ivm = (num_ivm_channels,
                                num_of_points_in_spike_trig_ivm,
                                num_of_spikes)
time_axis = np.arange(-num_of_points_in_spike_trig_ivm/(2*sampling_freq),
              num_of_points_in_spike_trig_ivm/(2*sampling_freq),
              1/sampling_freq)
ivm_data_filtered = np.memmap(os.path.join(analysis_folder, types_of_data_to_load[data_to_load].format(good_cells)),
                                dtype=filtered_data_type,
                                mode='r',
                                shape=shape_of_filt_spike_trig_ivm)

def test_filters():
    """Test low-, band-, high-pass, and band-stop filters plus resampling."""
    sfreq = 100
    sig_len_secs = 15

    a = rng.randn(2, sig_len_secs * sfreq)

    # let's test our catchers
    for fl in ['blah', [0, 1], 1000.5, '10ss', '10']:
        assert_raises(ValueError, band_pass_filter, a, sfreq, 4, 8, fl, 1., 1.)
    for nj in ['blah', 0.5]:
        assert_raises(ValueError, band_pass_filter, a, sfreq, 4, 8, 100,
                      1., 1., n_jobs=nj)
    assert_raises(ValueError, band_pass_filter, a, sfreq, 4, 8, 100,
                  1., 1., fir_window='foo')
    # > Nyq/2
    assert_raises(ValueError, band_pass_filter, a, sfreq, 4, sfreq / 2.,
                  100, 1.0, 1.0)
    assert_raises(ValueError, low_pass_filter, a, sfreq, sfreq / 2.,
                  100, 1.0)
    # check our short-filter warning:
    with warnings.catch_warnings(record=True) as w:
        # Warning for low attenuation
        band_pass_filter(a, sfreq, 1, 8, filter_length=256)
    assert_true(any('attenuation' in str(ww.message) for ww in w))
    with warnings.catch_warnings(record=True) as w:
        # Warning for too short a filter
        band_pass_filter(a, sfreq, 1, 8, filter_length='0.5s')
    assert_true(any('Increase filter_length' in str(ww.message) for ww in w))

    # try new default and old default
    for fl in ['auto', '10s', '5000ms', 1024]:
        bp = band_pass_filter(a, sfreq, 4, 8, fl, 1.0, 1.0)
        bs = band_stop_filter(a, sfreq, 4 - 1.0, 8 + 1.0, fl, 1.0, 1.0)
        lp = low_pass_filter(a, sfreq, 8, fl, 1.0, n_jobs=2)
        hp = high_pass_filter(lp, sfreq, 4, fl, 1.0)
        assert_array_almost_equal(hp, bp, 4)
        assert_array_almost_equal(bp + bs, a, 4)

    # and since these are low-passed, downsampling/upsampling should be close
    n_resamp_ignore = 10
    bp_up_dn = resample(resample(bp, 2, 1, n_jobs=2), 1, 2, n_jobs=2)
    assert_array_almost_equal(bp[n_resamp_ignore:-n_resamp_ignore],
                              bp_up_dn[n_resamp_ignore:-n_resamp_ignore], 2)
    # note that on systems without CUDA, this line serves as a test for a
    # graceful fallback to n_jobs=1
    bp_up_dn = resample(resample(bp, 2, 1, n_jobs='cuda'), 1, 2, n_jobs='cuda')
    assert_array_almost_equal(bp[n_resamp_ignore:-n_resamp_ignore],
                              bp_up_dn[n_resamp_ignore:-n_resamp_ignore], 2)
    # test to make sure our resamling matches scipy's
    bp_up_dn = sp_resample(sp_resample(bp, 2 * bp.shape[-1], axis=-1,
                                       window='boxcar'),
                           bp.shape[-1], window='boxcar', axis=-1)
    assert_array_almost_equal(bp[n_resamp_ignore:-n_resamp_ignore],
                              bp_up_dn[n_resamp_ignore:-n_resamp_ignore], 2)

    # make sure we don't alias
    t = np.array(list(range(sfreq * sig_len_secs))) / float(sfreq)
    # make sinusoid close to the Nyquist frequency
    sig = np.sin(2 * np.pi * sfreq / 2.2 * t)
    # signal should disappear with 2x downsampling
    sig_gone = resample(sig, 1, 2)[n_resamp_ignore:-n_resamp_ignore]
    assert_array_almost_equal(np.zeros_like(sig_gone), sig_gone, 2)

    # let's construct some filters
    iir_params = dict(ftype='cheby1', gpass=1, gstop=20, output='ba')
    iir_params = construct_iir_filter(iir_params, 40, 80, 1000, 'low')
    # this should be a third order filter
    assert_equal(iir_params['a'].size - 1, 3)
    assert_equal(iir_params['b'].size - 1, 3)
    iir_params = dict(ftype='butter', order=4, output='ba')
    iir_params = construct_iir_filter(iir_params, 40, None, 1000, 'low')
    assert_equal(iir_params['a'].size - 1, 4)
    assert_equal(iir_params['b'].size - 1, 4)
    iir_params = dict(ftype='cheby1', gpass=1, gstop=20, output='sos')
    iir_params = construct_iir_filter(iir_params, 40, 80, 1000, 'low')
    # this should be a third order filter, which requires 2 SOS ((2, 6))
    assert_equal(iir_params['sos'].shape, (2, 6))
    iir_params = dict(ftype='butter', order=4, output='sos')
    iir_params = construct_iir_filter(iir_params, 40, None, 1000, 'low')
    assert_equal(iir_params['sos'].shape, (2, 6))

    # check that picks work for 3d array with one channel and picks=[0]
    a = rng.randn(5 * sfreq, 5 * sfreq)
    b = a[:, None, :]

    a_filt = band_pass_filter(a, sfreq, 4, 8, 400, 2.0, 2.0)
    b_filt = band_pass_filter(b, sfreq, 4, 8, 400, 2.0, 2.0, picks=[0])

    assert_array_equal(a_filt[:, None, :], b_filt)

    # check for n-dimensional case
    a = rng.randn(2, 2, 2, 2)
    with warnings.catch_warnings(record=True):  # filter too long
        assert_raises(ValueError, band_pass_filter, a, sfreq, 4, 8, 100,
                      1.0, 1.0, picks=np.array([0, 1]))

def test_filters():
    """Test low-, band-, high-pass, and band-stop filters plus resampling
    """
    Fs = 500
    sig_len_secs = 30

    # Filtering of short signals (filter length = len(a))
    a = np.random.randn(2, sig_len_secs * Fs)
    bp = band_pass_filter(a, Fs, 4, 8)
    bs = band_stop_filter(a, Fs, 4 - 0.5, 8 + 0.5)
    lp = low_pass_filter(a, Fs, 8)
    hp = high_pass_filter(lp, Fs, 4)
    assert_array_almost_equal(hp, bp, 2)
    assert_array_almost_equal(bp + bs, a, 1)

    # Overlap-add filtering with a fixed filter length
    filter_length = 8192
    bp_oa = band_pass_filter(a, Fs, 4, 8, filter_length)
    bs_oa = band_stop_filter(a, Fs, 4 - 0.5, 8 + 0.5, filter_length)
    lp_oa = low_pass_filter(a, Fs, 8, filter_length)
    hp_oa = high_pass_filter(lp_oa, Fs, 4, filter_length)
    assert_array_almost_equal(hp_oa, bp_oa, 2)
    assert_array_almost_equal(bp_oa + bs_oa, a, 2)

    # The two methods should give the same result
    # As filtering for short signals uses a circular convolution (FFT) and
    # the overlap-add filter implements a linear convolution, the signal
    # boundary will be slightly different and we ignore it
    n_edge_ignore = 0
    assert_array_almost_equal(hp[n_edge_ignore:-n_edge_ignore],