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_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 bandPassFilterConcurrent(observations):
    """
    500Hz alpha 채널(8~13Hz)만 분리
    :param observations:
    :return:
    """
    return band_pass_filter(np.float64(observations),500,8,13)

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 filter_and_make_analytic_signal(data, sfreq, l_phase_freq, h_phase_freq,
                                    l_amp_freq, h_amp_freq, method='fft',
                                    n_jobs=1):
    """ Filter data to required range and compute analytic signal from it.

    Parameters
    ----------
    data : ndarray
        The signal to be analysed.
    l_phase_freq, h_phase_freq : float
        Low and high phase modulating frequencies.
    l_amp_freq, h_amp_freq : float
        Low and high amplitude modulated frequencies.
    method : 'fft' or 'iir'
        Filter method to be used. (mne.filter.band_pass_filter)
    n_jobs : int
        Number of parallel jobs to run.

    Returns
    -------
    theta : ndarray
        Low frequency filtered signal (modulating)
    gamma : ndarray
        High frequency filtered signal (modulated)
    phase : ndarray
        Phase of low frequency signal above.
    amp : ndarray
        Amplitude envelope of the high freq. signal above.
    """
    # filter theta and gamma signals
    n_jobs = 4
    method = 'fft'
    l_phase_freq, h_phase_freq, l_amp_freq, h_amp_freq = 6, 10, 60, 150

    theta = band_pass_filter(data, sfreq, l_phase_freq,
                             h_phase_freq, method=method, n_jobs=n_jobs)
    gamma = band_pass_filter(data, sfreq, l_amp_freq, h_amp_freq,
                             method=method, n_jobs=n_jobs)

    # phase of the low freq modulating signal
    phase = np.angle(hilbert(theta))
    # amplitude envelope of the high freq modulated signal
    amp = np.abs(hilbert(gamma))

    return theta, gamma, phase, amp

def modulation_index2d(data, sfreq):
    """ Compute the two dimensional modulation index.

    Parameters
    ----------
    data : ndarray
        The signal data
    sfreq: float
        Sampling frequency

    Returns
    -------
    mod2d : ndarray
        2 dimensional modulation index
    """

    from mne.filter import band_pass_filter
    from scipy.signal import hilbert

    flow = np.arange(2, 40, 1)
    flow_step = 1.0
    fhigh = np.arange(5, 205, 5)
    fhigh_step = 5.0

    mod2d = np.zeros((flow.size, fhigh.size))
    method = 'fft'
    n_jobs = 2

    for i in range(0, flow.size):
        theta = band_pass_filter(data, sfreq, flow[i], flow[i] + flow_step,
                                 method=method, n_jobs=n_jobs)
        theta = theta[sfreq: data.size - sfreq]
        phase = np.angle(hilbert(theta))

        for j in range(0, fhigh.size):
            gamma = band_pass_filter(data, sfreq, fhigh[j], fhigh[j] + fhigh_step,
                                     method=method, n_jobs=n_jobs)
            gamma = gamma[sfreq: data.size - sfreq]
            amp = np.abs(hilbert(gamma))

            # compute the modulation index
            m_norm_length = modulation_index1d(amp, phase, sfreq)
            mod2d[i, j] = m_norm_length

    return mod2d

 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 filterABDT(arrExp):
    frDelta = 1
    frTheta = 4
    frAlpha = 8
    frBeta = 13
    expnum = arrExp.shape[1]
    arrFiltered = np.ndarray((4,30,expnum,750), dtype=np.float64)
    for i in range(4): #alpha, beta, tetha, delta
        if i == 0: frFreq = frDelta; toFreq = frTheta;
        elif i == 1: frFreq = frTheta; toFreq = frAlpha;
        elif i == 2: frFreq = frAlpha; toFreq = frBeta;
        elif i == 3: frFreq = frBeta; toFreq = 20;
        for j in range(expnum): # number of experiments
            for k in range(arrExp.shape[0]): #channel
                arrFiltered[i,k,j,:] = band_pass_filter(np.float64(arrExp[k,j,:]),500,frFreq, toFreq)

    return arrFiltered

    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 filterByFrequency(arrExperiment,freqs):
    """
    실험단위별로 분리된 3차원 배열을 프리퀀시별로 분리하여 4차원 배열로 저장
    :param arrExperiment: 실험단위 별로 분리된 3차원 배열
    :return:
    """
    frDelta = 1
    frTheta = 4
    frAlpha = 8
    frBeta = 13
    expnum = arrExperiment.shape[1]
    obsnum = arrExperiment.shape[2]
    arrExperimentFrequency = np.ndarray((len(freqs),30,expnum,obsnum), dtype=np.float64)
    for i in range(len(freqs)): #alpha, beta, tetha, delta
        freq_val = freqs[i]
        if freq_val == 0: frFreq = frDelta; toFreq = frTheta;
        elif freq_val == 1: frFreq = frTheta; toFreq = frAlpha;
        elif freq_val == 2: frFreq = frAlpha; toFreq = frBeta;
        elif freq_val == 3: frFreq = frBeta; toFreq = 20;
        for j in range(expnum): # number of experiments
            for k in range(arrExperiment.shape[0]): #channel
                arrExperimentFrequency[i,k,j,:] = band_pass_filter(np.float64(arrExperiment[k,j,:]),500,frFreq, toFreq)

    return arrExperimentFrequency

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])

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

def filterAlphaPhase(_chnnDsetLst):
    return band_pass_filter(np.float64(_chnnDsetLst),500,8,13)

def mne_bandpass_filter(self):
        filter.band_pass_filter(self.data, self.fs, self.lowcut, self.highcut)
        return mne_filter

#         def __init__(self, a, b, signal):
#                 self.a = []
#                 self.b = []
#                 self.signal = signal
   
#         def apply(self, signal):
#                  return lfilter(self.b, self.a, signal)

# from scipy.io import loadmat
# mat = loadmat('bp_filter_coeff.mat')
# filters = [TimeDomainFilter(b=mat['bp_filter_coeff']['b'][0, 0].squeeze(),
#                                 a=mat['bp_filter_coeff']['a'][0, 0].squeeze()),
#         TimeDomainFilter(b=mat['bp_filter_coeff']['b_notch'][0, 0].squeeze(),
#                                  a=mat['bp_filter_coeff']['a_notch'][0, 0].squeeze()), ]
# for filter in filters:
#         y = TimeDomainFilter.apply(y)


hz_per_bin = float(250) / 256
data = np.genfromtxt('dan_pythondata.csv', delimiter=' ')
filt_data = filter.band_pass_filter(data[:, 1], 250, 1, 60)
psd, fft_data = windowed_fft(filt_data, 250)
psd_per_bin = psd / hz_per_bin


plt.figure(1)
plt.plot(fft_data, np.sqrt(psd_per_bin))

plt.figure(2)
plt.plot(filt_data)
plt.show() 

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

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

    # let's test our catchers
    for fl in ['blah', [0, 1], 1000.5, '10ss', '10']:
        assert_raises(ValueError, band_pass_filter, a, Fs, 4, 8,
                      filter_length=fl)
    for nj in ['blah', 0.5, 0]:
        assert_raises(ValueError, band_pass_filter, a, Fs, 4, 8, n_jobs=nj)
    assert_raises(ValueError, band_pass_filter, a, Fs, 4, Fs / 2.)  # > Nyq/2
    assert_raises(ValueError, low_pass_filter, a, Fs, Fs / 2.)  # > Nyq/2
    # check our short-filter warning:
    with warnings.catch_warnings(record=True) as w:
        # Warning for low attenuation
        band_pass_filter(a, Fs, 1, 8, filter_length=1024)
        # Warning for too short a filter
        band_pass_filter(a, Fs, 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, Fs, 4, 8, filter_length=fl)
        bs = band_stop_filter(a, Fs, 4 - 0.5, 8 + 0.5, filter_length=fl)
        lp = low_pass_filter(a, Fs, 8, filter_length=fl, n_jobs=2)
        hp = high_pass_filter(lp, Fs, 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, 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],
                              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(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)

    # 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)

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]))