Python scipy.signal.welch() Examples

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 

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 

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

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 

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 

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

def get_highest_periodicity(x): #highest peaks of fft
    from scipy.signal import welch
    length = len(x)
    w = welch(x, scaling='spectrum', nperseg=length // 2)
    peak_ind = get_peaks(w[1], 5)
    periods = [np.round(1 / w[0][peak]) for peak in peak_ind]
    #filter periods. Remove too small and too high
    max_p = length // 2
    min_p = 3  # Do not show short term AR(1) - AR(3)
    return tuple(period for period in periods if min_p < period < max_p) 

def estimate_sigma(fluor, range_ff=(0.25, 0.5), method='logmexp'):
    """
    Estimate noise power through the power spectral density over the range of
    large frequencies

    Parameters
    ----------
    fluor : nparray
        One dimensional array containing the fluorescence intensities with
        one entry per time-bin.
    range_ff : 2-tuple, optional, nonnegative, max value <= 0.5
        range of frequency (x Nyquist rate) over which the spectrum is
        averaged. Default: (0.25, 0.5)
    method : {'mean', 'median', 'logmexp'}, optional
        method of averaging: Mean, median, exponentiated mean of logvalues.
        Default: 'logmexp'

    Returns
    -------
    sigma       : noise standard deviation
    """

    if len(fluor) < 256:
        nperseg = len(fluor)
    else:
        nperseg = 256

    ff, Pxx = welch(fluor, nperseg=nperseg)
    ind1 = ff > range_ff[0]
    ind2 = ff < range_ff[1]
    ind = np.logical_and(ind1, ind2)
    Pxx_ind = Pxx[ind]
    sigma = {
        'mean': lambda Pxx_ind: np.sqrt(np.mean(Pxx_ind / 2)),
        'median': lambda Pxx_ind: np.sqrt(np.median(Pxx_ind / 2)),
        'logmexp': lambda Pxx_ind: np.sqrt(
            np.exp(np.mean(np.log(Pxx_ind / 2))))
    }[method](Pxx_ind)

    return sigma 

def calc_bands_power(x, dt, bands):
    from scipy.signal import welch
    f, psd = welch(x, fs=1. / dt)
    power = {band: np.mean(psd[np.where((f >= lf) & (f <= hf))]) for band, (lf, hf) in bands.items()}
    return power 

def getcardcoeffs(cardiacwaveform, slicesamplerate, minhr=40.0, maxhr=140.0, smoothlen=101, debug=False, display=False):
    if len(cardiacwaveform) > 1024:
        thex, they = welch(cardiacwaveform, slicesamplerate, nperseg=1024)
    else:
        thex, they = welch(cardiacwaveform, slicesamplerate)
    initpeakfreq = np.round(thex[np.argmax(they)] * 60.0, 2)
    if initpeakfreq > maxhr:
        initpeakfreq = maxhr
    if initpeakfreq < minhr:
        initpeakfreq = minhr
    if debug:
        print('initpeakfreq:', initpeakfreq, 'BPM')
    freqaxis, spectrum = tide_filt.spectrum(tide_filt.hamming(len(cardiacwaveform)) * cardiacwaveform,
                                            Fs=slicesamplerate,
                                            mode='complex')
    # remove any spikes at zero frequency
    minbin = int(minhr // (60.0 * (freqaxis[1] - freqaxis[0])))
    maxbin = int(maxhr // (60.0 * (freqaxis[1] - freqaxis[0])))
    spectrum[:minbin] = 0.0
    spectrum[maxbin:] = 0.0

    # find the max
    ampspec = savgolsmooth(np.abs(spectrum), smoothlen=smoothlen)
    if display:
        figure()
        plot(freqaxis, ampspec, 'r')
        show()
    peakfreq = freqaxis[np.argmax(ampspec)]
    if debug:
        print('cardiac fundamental frequency is', np.round(peakfreq * 60.0, 2), 'BPM')
    normfac = np.sqrt(2.0) * tide_math.rms(cardiacwaveform)
    if debug:
        print('normfac:', normfac)
    return peakfreq 

def noise_welch(y, noise_range, noise_method):
    ff, Pxx = welch(y)
    mask0, mask1 = ff > noise_range[0], ff < noise_range[1]
    mask = np.logical_and(mask0, mask1)
    Pxx_ind = Pxx[mask]
    sn = {
        'mean': lambda x: np.sqrt(np.mean(x / 2)),
        'median': lambda x: np.sqrt(np.median(x / 2)),
        'logmexp': lambda x: np.sqrt(np.exp(np.mean(np.log(x / 2))))
    }[noise_method](Pxx_ind)
    return sn 

def _welch(x, **kwargs):
    return welch(x, **kwargs)[1] 

def display(self, data):
        """Make graphicsitem for spectrum figure.

        Parameters
        ----------
        data : ndarray
            1D vector containing the data only

        This function can be called by self.display_window (which reads the
        data for the selected channel) or by the mouse-events functions in
        traces (which read chunks of data from the user-made selection).
        """
        value = self.config.value
        self.scene = QGraphicsScene(value['x_min'], value['y_min'],
                                    value['x_max'] - value['x_min'],
                                    value['y_max'] - value['y_min'])
        self.idx_fig.setScene(self.scene)

        self.add_grid()
        self.resizeEvent(None)

        s_freq = self.parent.traces.data.s_freq
        f, Pxx = welch(data, fs=s_freq,
                       nperseg=int(min((s_freq, len(data)))))  # force int

        freq_limit = (value['x_min'] <= f) & (f <= value['x_max'])

        if self.config.value['log']:
            Pxx_to_plot = log(Pxx[freq_limit])
        else:
            Pxx_to_plot = Pxx[freq_limit]

        self.scene.addPath(Path(f[freq_limit], Pxx_to_plot),
                           QPen(QColor(LINE_COLOR), LINE_WIDTH)) 

def psd(x, fs=2, nfft=512, noverlap=None, window='hanning', color=None, style='solid', thickness=1, marker=None, filled=False, size=6, title=None, xlabel='Frequency (Hz)', ylabel='Power spectral density (dB/Hz)', xlim=None, ylim=None, width=None, height=None, legend=None, hold=False, interactive=None):
    """Plot power spectral density of a given time series signal.

    :param x: time series signal
    :param fs: sampling rate
    :param nfft: segment size (see `scipy.signal.welch <https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html>`_)
    :param noverlap: overlap size (see `scipy.signal.welch`_)
    :param window: window to use (see `scipy.signal.welch`_)
    :param color: line color (see `Bokeh colors`_)
    :param style: line style ('solid', 'dashed', 'dotted', 'dotdash', 'dashdot')
    :param thickness: line width in pixels
    :param marker: point markers ('.', 'o', 's', '*', 'x', '+', 'd', '^')
    :param filled: filled markers or outlined ones
    :param size: marker size
    :param title: figure title
    :param xlabel: x-axis label
    :param ylabel: y-axis label
    :param xlim: x-axis limits (min, max)
    :param ylim: y-axis limits (min, max)
    :param width: figure width in pixels
    :param height: figure height in pixels
    :param legend: legend text
    :param interactive: enable interactive tools (pan, zoom, etc) for plot
    :param hold: if set to True, output is not plotted immediately, but combined with the next plot

    >>> import arlpy.plot
    >>> import numpy as np
    >>> arlpy.plot.psd(np.random.normal(size=(10000)), fs=10000)
    """
    f, Pxx = _sig.welch(x, fs=fs, nperseg=nfft, noverlap=noverlap, window=window)
    Pxx = 10*_np.log10(Pxx+_np.finfo(float).eps)
    if xlim is None:
        xlim = (0, fs/2)
    if ylim is None:
        ylim = (_np.max(Pxx)-50, _np.max(Pxx)+10)
    plot(f, Pxx, color=color, style=style, thickness=thickness, marker=marker, filled=filled, size=size, title=title, xlabel=xlabel, ylabel=ylabel, xlim=xlim, ylim=ylim, maxpts=len(f), width=width, height=height, hold=hold, legend=legend, interactive=interactive) 

def test_nfft_too_short(self):
        assert_raises(ValueError, welch, np.ones(12), nfft=3, nperseg=4) 

def test_nondefault_noverlap(self):
        x = np.zeros(64)
        x[::8] = 1
        f, p = welch(x, nperseg=16, noverlap=4)
        q = np.array([0, 1./12., 1./3., 1./5., 1./3., 1./5., 1./3., 1./5., 1./6.])
        assert_allclose(p, q, atol=1e-12) 

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) 

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

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

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

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

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) 

def test_window_long_or_nd(self):
        with warnings.catch_warnings():
            warnings.simplefilter('ignore', UserWarning)
            assert_raises(ValueError, welch, np.zeros(4), 1, np.array([1,1,1,1,1]))
            assert_raises(ValueError, welch, np.zeros(4), 1,
                          np.arange(6).reshape((2,3))) 

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) 

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) 

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) 

def test_window_external(self):
        x = np.zeros(16)
        x[0] = 1
        x[8] = 1
        f, p = welch(x, 10, 'hanning', 8)
        win = signal.get_window('hanning', 8)
        fe, pe = welch(x, 10, win, 8)
        assert_array_almost_equal_nulp(p, pe)
        assert_array_almost_equal_nulp(f, fe) 

def test_empty_input(self):
        f, p = welch([])
        assert_array_equal(f.shape, (0,))
        assert_array_equal(p.shape, (0,))
        for shape in [(0,), (3,0), (0,5,2)]:
            f, p = welch(np.empty(shape))
            assert_array_equal(f.shape, shape)
            assert_array_equal(p.shape, shape) 

def test_short_data(self):
        x = np.zeros(8)
        x[0] = 1
        with warnings.catch_warnings():
            warnings.simplefilter('ignore', UserWarning)
            f, p = welch(x)

        f1, p1 = welch(x, nperseg=8)
        assert_allclose(f, f1)
        assert_allclose(p, p1) 

def _get_freq_psd_from_nn_intervals(nn_intervals: List[float], method: str = WELCH_METHOD,
                                    sampling_frequency: int = 4,
                                    interpolation_method: str = "linear",
                                    vlf_band: namedtuple = VlfBand(0.003, 0.04),
                                    hf_band: namedtuple = HfBand(0.15, 0.40)) -> Tuple:
    """
    Returns the frequency and power of the signal.

    Parameters
    ---------
    nn_intervals : list
        list of Normal to Normal Interval
    method : str
        Method used to calculate the psd. Choice are Welch's FFT or Lomb method.
    sampling_frequency : int
        Frequency at which the signal is sampled. Common value range from 1 Hz to 10 Hz,
        by default set to 7 Hz. No need to specify if Lomb method is used.
    interpolation_method : str
        Kind of interpolation as a string, by default "linear". No need to specify if Lomb
        method is used.
    vlf_band : tuple
        Very low frequency bands for features extraction from power spectral density.
    hf_band : tuple
        High frequency bands for features extraction from power spectral density.

    Returns
    ---------
    freq : list
        Frequency of the corresponding psd points.
    psd : list
        Power Spectral Density of the signal.
    """

    timestamp_list = _create_timestamp_list(nn_intervals)

    if method == WELCH_METHOD:
        # ---------- Interpolation of signal ---------- #
        funct = interpolate.interp1d(x=timestamp_list, y=nn_intervals, kind=interpolation_method)

        timestamps_interpolation = _create_interpolated_timestamp_list(nn_intervals, sampling_frequency)
        nni_interpolation = funct(timestamps_interpolation)

        # ---------- Remove DC Component ---------- #
        nni_normalized = nni_interpolation - np.mean(nni_interpolation)

        #  --------- Compute Power Spectral Density  --------- #
        freq, psd = signal.welch(x=nni_normalized, fs=sampling_frequency, window='hann',
                                 nfft=4096)

    elif method == LOMB_METHOD:
        freq, psd = LombScargle(timestamp_list, nn_intervals,
                                normalization='psd').autopower(minimum_frequency=vlf_band[0],
                                                               maximum_frequency=hf_band[1])
    else:
        raise ValueError("Not a valid method. Choose between 'lomb' and 'welch'")

    return freq, psd