Python filter.resample函数代码示例

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)

def test_resample():
    "Test resampling"
    x = np.random.normal(0, 1, (10, 10, 10))
    x_rs = resample(x, 1, 2, 10)
    assert_equal(x.shape, (10, 10, 10))
    assert_equal(x_rs.shape, (10, 10, 5))

    x_2 = x.swapaxes(0, 1)
    x_2_rs = resample(x_2, 1, 2, 10)
    assert_array_equal(x_2_rs.swapaxes(0, 1), x_rs)

    x_3 = x.swapaxes(0, 2)
    x_3_rs = resample(x_3, 1, 2, 10, 0)
    assert_array_equal(x_3_rs.swapaxes(0, 2), x_rs)

def test_resample_scipy():
    """Test resampling against SciPy."""
    n_jobs_test = (1, 'cuda')
    for window in ('boxcar', 'hann'):
        for N in (100, 101, 102, 103):
            x = np.arange(N).astype(float)
            err_msg = '%s: %s' % (N, window)
            x_2_sp = sp_resample(x, 2 * N, window=window)
            for n_jobs in n_jobs_test:
                x_2 = resample(x, 2, 1, 0, window=window, n_jobs=n_jobs)
                assert_allclose(x_2, x_2_sp, atol=1e-12, err_msg=err_msg)
            new_len = int(round(len(x) * (1. / 2.)))
            x_p5_sp = sp_resample(x, new_len, window=window)
            for n_jobs in n_jobs_test:
                x_p5 = resample(x, 1, 2, 0, window=window, n_jobs=n_jobs)
                assert_allclose(x_p5, x_p5_sp, atol=1e-12, err_msg=err_msg)

def test_resample():
    """Test resampling."""
    x = rng.normal(0, 1, (10, 10, 10))
    x_rs = resample(x, 1, 2, 10)
    assert_equal(x.shape, (10, 10, 10))
    assert_equal(x_rs.shape, (10, 10, 5))

    x_2 = x.swapaxes(0, 1)
    x_2_rs = resample(x_2, 1, 2, 10)
    assert_array_equal(x_2_rs.swapaxes(0, 1), x_rs)

    x_3 = x.swapaxes(0, 2)
    x_3_rs = resample(x_3, 1, 2, 10, 0)
    assert_array_equal(x_3_rs.swapaxes(0, 2), x_rs)

    # make sure we cast to array if necessary
    assert_array_equal(resample([0, 0], 2, 1), [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_cuda_resampling():
    """Test CUDA 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.])
    from mne.cuda import _cuda_capable  # allow above funs to set it
    if not _cuda_capable:
        pytest.skip('CUDA not enabled')

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 _prepare_stim(zfn, path_out, sex, tal, cal, col, num, fs_out, dtype,
                  ref_rms, n_jobs):
    """Read in a binary CRM file and write out a scaled resampled wav.
    """
    from mne.filter import resample
    zip_file = ZipFile(zfn)
    x = _read_binary(zip_file, cal, col, num, 0)
    fn = '%i%i%i%i%i.wav' % (sex, tal, cal, col, num)
    if int(np.round(fs_out)) != int(np.round(_fs_binary)):
        x = resample(x, fs_out, _fs_binary, n_jobs=n_jobs, verbose=0)
    x *= ref_rms / _rms_binary
    write_wav(join(path_out, fn), x, fs_out, overwrite=True, dtype=dtype,
              verbose=False)

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

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 convolve_hrtf(data, fs, angle, source='cipic', interp=False):
    """Convolve a signal with a head-related transfer function

    Technically we will be convolving with binaural room impulse
    responses (BRIRs), but HRTFs (freq-domain equiv. representations)
    are the common terminology.

    Parameters
    ----------
    data : 1-dimensional or 1xN array-like
        Data to operate on.
    fs : float
        The sample rate of the data. (HRTFs will be resampled if necessary.)
    angle : float
        The azimuthal angle of the HRTF.
    source : str
        Source to use for HRTFs. Currently `'barb'` and `'cipic'` are
        supported. The former is default for legacy purpose. The latter is
        recommended for new experiments.
    interp : bool
        Parameter to determine whether to restrict use to known HRTF values or 
        to use an interpolated HRTF for angles not in the source; set to 
        False by default

    Returns
    -------
    data_hrtf : array
        A 2D array ``shape=(2, n_samples)`` containing the convolved data.

    Notes
    -----
    CIPIC data downloaded from:

        http://earlab.bu.edu/databases/collections/cipic/Default.aspx.

    Additional documentation:

        http://earlab.bu.edu/databases/collections/cipic/documentation/hrir_data_documentation.pdf  # noqa

    The data were modified to suit our experimental needs. Below is the
    licensing information for the CIPIC data:

    **Copyright**

    Copyright (c) 2001 The Regents of the University of California. All Rights
    Reserved.

    **Disclaimer**

    THE REGENTS OF THE UNIVERSITY OF CALIFORNIA MAKE NO REPRESENTATION OR
    WARRANTIES WITH RESPECT TO THE CONTENTS HEREOF AND SPECIFICALLY DISCLAIM
    ANY IMPLIED WARRANTIES OR MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR
    PURPOSE.

    Further, the Regents of the University of California reserve the right to
    revise this software and/or documentation and to make changes from time to
    time in the content hereof without obligation of the Regents of the
    University of California to notify any person of such revision or change.

    Use of Materials
    The Regents of the University of California hereby grant users permission
    to reproduce and/or use materials available therein for any purpose-
    educational, research or commercial. However, each reproduction of any part
    of the materials must include the copyright notice, if it is present.

    In addition, as a courtesy, if these materials are used in published
    research, this use should be acknowledged in the publication. If these
    materials are used in the development of commercial products, the Regents
    of the University of California request that written acknowledgment of such
    use be sent to:

    CIPIC- Center for Image Processing and Integrated Computing University of
    California 1 Shields Avenue Davis, CA 95616-8553
    """
    fs = float(fs)
    angle = float(angle)
    known_sources = ['barb', 'cipic']
    known_fs = [24414, 44100]  # must be sorted
    if source not in known_sources:
        raise ValueError('Source "{0}" unknown, must be one of {1}'
                         ''.format(source, known_sources))
    if not isinstance(interp, bool):
        raise ValueError('interp must be bool')
    data = np.array(data, np.float64)
    data = _fix_audio_dims(data, n_channels=1).ravel()

    # Find out which sampling rate to get--first that is >= fs
    # Use the last, highest one whether it is high enough or not
    ge = [int(np.round(fs)) <= k for k in known_fs[:-1]] + [True]
    brir_fs = known_fs[ge.index(True)]

    brir, brir_fs, leftward = _get_hrtf(angle, source, brir_fs, interp)
    order = [1, 0] if leftward else [0, 1]
    if not np.allclose(brir_fs, fs, rtol=0, atol=0.5):
        from mne.filter import resample
        brir = [resample(b, fs, brir_fs) for b in brir]
    out = np.array([np.convolve(data, brir[o]) for o in order])
    return out

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],
                              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 * len(bp_oa), window='boxcar'),
                           len(bp_oa), window='boxcar')
    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)

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