Python signal.welch方法代码示例

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def rmsmap(fbin):
    """
    Computes RMS map in time domain and spectra for each channel of Neuropixel probe

    :param fbin: binary file in spike glx format (will look for attached metatdata)
    :type fbin: str or pathlib.Path
    :return: a dictionary with amplitudes in channeltime space, channelfrequency space, time
     and frequency scales
    """
    if not isinstance(fbin, spikeglx.Reader):
        sglx = spikeglx.Reader(fbin)
    rms_win_length_samples = 2 ** np.ceil(np.log2(sglx.fs * RMS_WIN_LENGTH_SECS))
    # the window generator will generates window indices
    wingen = dsp.WindowGenerator(ns=sglx.ns, nswin=rms_win_length_samples, overlap=0)
    # pre-allocate output dictionary of numpy arrays
    win = {'TRMS': np.zeros((wingen.nwin, sglx.nc)),
           'nsamples': np.zeros((wingen.nwin,)),
           'fscale': dsp.fscale(WELCH_WIN_LENGTH_SAMPLES, 1 / sglx.fs, one_sided=True),
           'tscale': wingen.tscale(fs=sglx.fs)}
    win['spectral_density'] = np.zeros((len(win['fscale']), sglx.nc))
    # loop through the whole session
    for first, last in wingen.firstlast:
        D = sglx.read_samples(first_sample=first, last_sample=last)[0].transpose()
        # remove low frequency noise below 1 Hz
        D = dsp.hp(D, 1 / sglx.fs, [0, 1])
        iw = wingen.iw
        win['TRMS'][iw, :] = dsp.rms(D)
        win['nsamples'][iw] = D.shape[1]
        # the last window may be smaller than what is needed for welch
        if last - first < WELCH_WIN_LENGTH_SAMPLES:
            continue
        # compute a smoothed spectrum using welch method
        _, w = signal.welch(D, fs=sglx.fs, window='hanning', nperseg=WELCH_WIN_LENGTH_SAMPLES,
                            detrend='constant', return_onesided=True, scaling='density', axis=-1)
        win['spectral_density'] += w.T
        # print at least every 20 windows
        if (iw % min(20, max(int(np.floor(wingen.nwin / 75)), 1))) == 0:
            print_progress(iw, wingen.nwin)
    return win 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def check_signal_power_signal1_below_signal2(signals1, signals2, freq_range, fs):
    '''
    Check that spectrum power of signal1 is below the one of signal 2 in the range freq_range
    '''
    f1, pow1 = ss.welch(signals1, fs, nfft=1024)
    f2, pow2 = ss.welch(signals2, fs, nfft=1024)

    below = True

    for (p1, p2) in zip(pow1, pow2):

        r1_idxs = np.where((f1 > freq_range[0]) & (f1 <= freq_range[1]))
        r2_idxs = np.where((f2 > freq_range[0]) & (f2 <= freq_range[1]))

        sump1 = np.sum(p1[r1_idxs])
        sump2 = np.sum(p2[r2_idxs])

        if sump1 >= sump2:
            below = False
            break

    return below 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def welchogram(fs, wav, nswin=NS_WIN, overlap=OVERLAP, nperseg=NS_WELCH):
    """
    Computes a spectrogram on a very large audio file.

    :param fs: sampling frequency (Hz)
    :param wav: wav signal (vector or memmap)
    :param nswin: n samples of the sliding window
    :param overlap: n samples of the overlap between windows
    :param nperseg: n samples for the computation of the spectrogram
    :return: tscale, fscale, downsampled_spectrogram
    """
    ns = wav.shape[0]
    window_generator = dsp.WindowGenerator(ns=ns, nswin=nswin, overlap=overlap)
    nwin = window_generator.nwin
    fscale = dsp.fscale(nperseg, 1 / fs, one_sided=True)
    W = np.zeros((nwin, len(fscale)))
    tscale = window_generator.tscale(fs=fs)
    detect = []
    for first, last in window_generator.firstlast:
        # load the current window into memory
        w = np.float64(wav[first:last]) * _get_conversion_factor()
        # detection of ready tones
        a = [d + first for d in _detect_ready_tone(w, fs)]
        if len(a):
            detect += a
        # the last window may not allow a pwelch
        if (last - first) < nperseg:
            continue
        # compute PSD estimate for the current window
        iw = window_generator.iw
        _, W[iw, :] = signal.welch(w, fs=fs, window='hanning', nperseg=nperseg, axis=-1,
                                   detrend='constant', return_onesided=True, scaling='density')
        if (iw % 50) == 0:
            window_generator.print_progress()
    window_generator.print_progress()
    # the onset detection may have duplicates with sliding window, average them and remove
    detect = np.sort(np.array(detect)) / fs
    ind = np.where(np.diff(detect) < 0.1)[0]
    detect[ind] = (detect[ind] + detect[ind + 1]) / 2
    detect = np.delete(detect, ind + 1)
    return tscale, fscale, W, detect 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_psd():
    n_channels, n_times = data.shape
    _data = data[None, ...]
    # Only test output shape when `method='welch'` or `method='multitaper'`
    # since it is actually just a wrapper for MNE functions:
    psd_welch, _ = power_spectrum(sfreq, _data, psd_method='welch')
    psd_multitaper, _ = power_spectrum(sfreq, _data, psd_method='multitaper')
    psd_fft, freqs_fft = power_spectrum(sfreq, _data, psd_method='fft')
    assert_equal(psd_welch.shape, (1, n_channels, n_times // 2 + 1))
    assert_equal(psd_multitaper.shape, (1, n_channels, n_times // 2 + 1))
    assert_equal(psd_fft.shape, (1, n_channels, n_times // 2 + 1))

    # Compare result obtained with `method='fft'` to the Scipy's result
    # (implementation of Welch's method with rectangular window):
    expected_freqs, expected_psd = signal.welch(data, sfreq,
                                                window=signal.get_window(
                                                    'boxcar', data.shape[-1]),
                                                return_onesided=True,
                                                scaling='density')
    assert_almost_equal(expected_freqs, freqs_fft)
    assert_almost_equal(expected_psd, psd_fft[0, ...]) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def process(self):
        self.update_counter += 1
        if self.update_counter == self.update_rate:
            self.update_counter = 0
        else:
            return


        if self.welch_button.checkState():
            f, C = welch(self.input.buffer, fs=250, nperseg=self.window_size, scaling='spectrum')
        else:
            C = np.fft.rfft(self.input.buffer[-self.window_size:] * self.window)
            C = abs(C)
            #C = C*C
        C = C[self.lo_index: self.hi_index]

        if self.logarithm:
            C = np.log(C)

        # Roll down one and replace leading edge with new data
        self.waterfallImgArray = np.roll(self.waterfallImgArray, -1, axis=0)
        self.waterfallImgArray[-1] = C 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_real_onesided_even(self):
        x = np.zeros(16)
        x[0] = 1
        x[8] = 1
        f, p = welch(x, nperseg=8)
        assert_allclose(f, np.linspace(0, 0.5, 5))
        assert_allclose(p, np.array([0.08333333, 0.15277778, 0.22222222,
            0.22222222, 0.11111111])) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_real_onesided_odd(self):
        x = np.zeros(16)
        x[0] = 1
        x[8] = 1
        f, p = welch(x, nperseg=9)
        assert_allclose(f, np.arange(5.0)/9.0)
        assert_allclose(p, np.array([0.15958226, 0.24193954, 0.24145223,
            0.24100919, 0.12188675])) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_real_twosided(self):
        x = np.zeros(16)
        x[0] = 1
        x[8] = 1
        f, p = welch(x, nperseg=8, return_onesided=False)
        assert_allclose(f, fftpack.fftfreq(8, 1.0))
        assert_allclose(p, np.array([0.08333333, 0.07638889, 0.11111111,
            0.11111111, 0.11111111, 0.11111111, 0.11111111, 0.07638889])) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_real_spectrum(self):
        x = np.zeros(16)
        x[0] = 1
        x[8] = 1
        f, p = welch(x, nperseg=8, scaling='spectrum')
        assert_allclose(f, np.linspace(0, 0.5, 5))
        assert_allclose(p, np.array([0.015625, 0.028645833333333332,
            0.041666666666666664, 0.041666666666666664, 0.020833333333333332])) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_complex(self):
        x = np.zeros(16, np.complex128)
        x[0] = 1.0 + 2.0j
        x[8] = 1.0 + 2.0j
        f, p = welch(x, nperseg=8)
        assert_allclose(f, fftpack.fftfreq(8, 1.0))
        assert_allclose(p, np.array([0.41666667, 0.38194444, 0.55555556,
            0.55555556, 0.55555556, 0.55555556, 0.55555556, 0.38194444])) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_detrend_linear(self):
        x = np.arange(10, dtype=np.float64)+0.04
        f, p = welch(x, nperseg=10, detrend='linear')
        assert_allclose(p, np.zeros_like(p), atol=1e-15) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_detrend_external(self):
        x = np.arange(10, dtype=np.float64)+0.04
        f, p = welch(x, nperseg=10,
                detrend=lambda seg: signal.detrend(seg, type='l'))
        assert_allclose(p, np.zeros_like(p), atol=1e-15) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_detrend_external_nd_m1(self):
        x = np.arange(40, dtype=np.float64)+0.04
        x = x.reshape((2,2,10))
        f, p = welch(x, nperseg=10,
                detrend=lambda seg: signal.detrend(seg, type='l'))
        assert_allclose(p, np.zeros_like(p), atol=1e-15) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_detrend_external_nd_0(self):
        x = np.arange(20, dtype=np.float64)+0.04
        x = x.reshape((2,1,10))
        x = np.rollaxis(x, 2, 0)
        f, p = welch(x, nperseg=10, axis=0,
                detrend=lambda seg: signal.detrend(seg, axis=0, type='l'))
        assert_allclose(p, np.zeros_like(p), atol=1e-15) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import welch [as 别名]
def test_nd_axis_m1(self):
        x = np.arange(20, dtype=np.float64)+0.04
        x = x.reshape((2,1,10))
        f, p = welch(x, nperseg=10)
        assert_array_equal(p.shape, (2, 1, 6))
        assert_allclose(p[0,0,:], p[1,0,:], atol=1e-13, rtol=1e-13)
        f0, p0 = welch(x[0,0,:], nperseg=10)
        assert_allclose(p0[np.newaxis,:], p[1,:], atol=1e-13, rtol=1e-13)