Python scipy.signal.spectrogram() Examples

def __init__(self, method='HouseDetector', duration=None):

        self.method = method

        if method == 'HouseDetector':
            self.freq_band1 = (5, None)
            self.freq_band2 = (0.2, None)
            self.spectrogram = {'dur': 1,
                                'overlap': 0.5,
                                'detrend': 'linear'}
            self.det_thresh = 1.2
            self.det_thresh_end = 1.1
            self.min_interval = 10
            self.duration = (3, 30)

        else:
            raise ValueError('Unknown method')
            
        if duration is None:
            self.duration = (self.min_dur, self.max_dur)       
        else:
            self.duration = duration 

def createSpec(data):
    fs=256
    lowcut=117
    highcut=123

    y=butter_bandstop_filter(data, lowcut, highcut, fs, order=6)
    lowcut=57
    highcut=63
    y=butter_bandstop_filter(y, lowcut, highcut, fs, order=6)
    
    cutoff=1
    y=butter_highpass_filter(y, cutoff, fs, order=6)
    
    Pxx=signal.spectrogram(y, nfft=256, fs=256, return_onesided=True, noverlap=128)[2]    
    Pxx = np.delete(Pxx, np.s_[117:123+1], axis=0)
    Pxx = np.delete(Pxx, np.s_[57:63+1], axis=0)
    Pxx = np.delete(Pxx, 0, axis=0)
    
    result=(10*np.log10(np.transpose(Pxx))-(10*np.log10(np.transpose(Pxx))).min())/(10*np.log10(np.transpose(Pxx))).ptp()
    return result

# Creazione spettrogramma e visualizzazione con la libreria matplotlib 

def show_plot(
    key, segment_times, sample_freqs, spec, duration, wav_data, vad_feat
):
    """This function plots the vad against the signal and the spectrogram.

    Args:
        segment_times: the time intervals acting as the x axis
        sample_freqs: the frequency bins acting as the y axis
        spec: the spectrogram
        duration: duration of the wave file
        wav_data: the wave data
        vad_feat: VAD features
    """

    import matplotlib.pyplot as plt
    import matplotlib.mlab as mlb

    plt.subplot(3, 1, 1)
    plt.pcolormesh(segment_times, sample_freqs, 10 * np.log10(spec), cmap="jet")
    plt.ylabel("Frequency [Hz]")
    plt.xlabel("Time [sec]")

    plt.subplot(3, 1, 2)
    axes = plt.gca()
    axes.set_xlim([0, duration])
    tmp_axis = np.linspace(0, duration, wav_data.shape[0])
    plt.plot(tmp_axis, wav_data / np.abs(np.max(wav_data)))
    plt.xlabel("Time [sec]")

    plt.subplot(3, 1, 3)
    axes = plt.gca()
    axes.set_xlim([0, duration])
    tmp_axis = np.linspace(0, duration, vad_feat.shape[0])
    plt.plot(tmp_axis, vad_feat)
    plt.xlabel("Time [sec]")

    plt.savefig("plots/" + key, bbox_inches="tight") 

def specgram(x, fs=2, nfft=None, noverlap=None, colormap='Plasma256', clim=None, clabel='dB', title=None, xlabel='Time (s)', ylabel='Frequency (Hz)', xlim=None, ylim=None, width=None, height=None, hold=False, interactive=None):
    """Plot spectrogram of a given time series signal.

    :param x: time series signal
    :param fs: sampling rate
    :param nfft: FFT size (see `scipy.signal.spectrogram <https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.spectrogram.html>`_)
    :param noverlap: overlap size (see `scipy.signal.spectrogram`_)
    :param colormap: named color palette or Bokeh ColorMapper (see `Bokeh palettes`_)
    :param clim: color axis limits (min, max), or dynamic range with respect to maximum
    :param clabel: color axis label
    :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 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.specgram(np.random.normal(size=(10000)), fs=10000, clim=30)
    """
    f, t, Sxx = _sig.spectrogram(x, fs=fs, nperseg=nfft, noverlap=noverlap)
    Sxx = 10*_np.log10(Sxx+_np.finfo(float).eps)
    if isinstance(clim, float) or isinstance(clim, int):
        clim = (_np.max(Sxx)-clim, _np.max(Sxx))
    image(Sxx, x=(t[0], t[-1]), y=(f[0], f[-1]), title=title, colormap=colormap, clim=clim, clabel=clabel, xlabel=xlabel, ylabel=ylabel, xlim=xlim, ylim=ylim, width=width, height=height, hold=hold, interactive=interactive) 

def computePitch(cPitchTrackName, afAudioData, f_s, afWindow=None, iBlockLength=4096, iHopLength=2048):
    
    #mypackage = __import__(".Pitch" + cPitchTrackName, package="pyACA")
    hPitchFunc = getattr(pyACA, "Pitch" + cPitchTrackName)

    # pre-processing
    afAudioData = ToolPreprocAudio(afAudioData, iBlockLength)

    if isSpectral(cPitchTrackName):
        # compute window function for FFT
        if afWindow is None:
            afWindow = ToolComputeHann(iBlockLength)

        assert(afWindow.shape[0] == iBlockLength), "parameter error: invalid window dimension"

        # in the real world, we would do this block by block...
        [f_k, t, X] = spectrogram(afAudioData,
                                  f_s,
                                  afWindow,
                                  iBlockLength,
                                  iBlockLength - iHopLength,
                                  iBlockLength,
                                  False,
                                  True,
                                  'spectrum')

        # we just want the magnitude spectrum...
        X = np.sqrt(X / 2)

        # compute instantaneous pitch chroma
        f = hPitchFunc(X, f_s)

    if isTemporal(cPitchTrackName):
        [f, t] = hPitchFunc(afAudioData, iBlockLength, iHopLength, f_s)

    return (f, t)


#######################################################
# helper functions 

def computeFeature(cFeatureName, afAudioData, f_s, afWindow=None, iBlockLength=4096, iHopLength=2048):
 
    #mypackage = __import__(".Feature" + cFeatureName, package="pyACA")
    hFeatureFunc = getattr(pyACA, "Feature" + cFeatureName)

    # pre-processing
    afAudioData = ToolPreprocAudio(afAudioData, iBlockLength)

    if isSpectral(cFeatureName):
        # compute window function for FFT
        if afWindow is None:
            afWindow = ToolComputeHann(iBlockLength)

        assert(afWindow.shape[0] == iBlockLength), "parameter error: invalid window dimension"

        # in the real world, we would do this block by block...
        [f, t, X] = spectrogram(afAudioData,
                                f_s,
                                afWindow,
                                iBlockLength,
                                iBlockLength - iHopLength,
                                iBlockLength,
                                False,
                                True,
                                'spectrum')

        # we just want the magnitude spectrum...
        X = np.sqrt(X / 2)

        # compute instantaneous feature
        v = hFeatureFunc(X, f_s)

    if isTemporal(cFeatureName):
        [v, t] = hFeatureFunc(afAudioData, iBlockLength, iHopLength, f_s)
        # [v, t] = hFeatureFunc(afAudioData, iBlockLength, iHopLength, f_s, np.array([2, 3]))

    return (v, t)


#######################################################
# helper functions 

def __init__(self, **kwargs):
        self.breaks = kwargs['breaks']
        self.spectrogram = kwargs.get('spectrogram', False)
        self.output_dir = kwargs.get('output_dir', None)

        super().__init__(**kwargs) 

def power_spectrum(signal: np.ndarray,
                   fs: int,
                   window_width: int,
                   window_overlap: int) -> (np.ndarray, np.ndarray, np.ndarray):
    """
    Computes the power spectrum of the specified signal.
    
    A periodic Hann window with the specified width and overlap is used.
    
    Parameters
    ----------
    signal: numpy.ndarray
        The input signal
    fs: int
        Sampling frequency of the input signal
    window_width: int
        Width of the Hann windows in samples
    window_overlap: int
        Overlap between Hann windows in samples

    Returns
    -------
    f: numpy.ndarray
        Array of frequency values for the first axis of the returned spectrogram
    t: numpy.ndarray
        Array of time values for the second axis of the returned spectrogram
    sxx: numpy.ndarray
        Power spectrogram of the input signal with axes [frequency, time]
    """
    f, t, sxx = spectrogram(x=signal,
                            fs=fs,
                            window=hann(window_width, sym=False),
                            noverlap=window_overlap,
                            mode="magnitude")

    return f, t, (1.0 / window_width) * (sxx ** 2) 

def mel_spectrum(power_spectrum: np.ndarray,
                 mel_fbank: np.ndarray = None,
                 fs: int = None,
                 window_width: int = None,
                 n_filt: int = 40) -> np.ndarray:
    """
    Computes a Mel spectrogram from the specified power spectrogram.
    
    Optionally, precomputed Mel filter banks can be passed to this function, in which case the n_filt, fs, and 
    window_width parameters are ignored. If precomputed Mel filter banks are used, the caller has to ensure that they
    have correct shape.
    
    Parameters
    ----------
    power_spectrum: numpy.ndarray
        The power spectrogram from which a Mel spectrogram should be computed
    mel_fbank: numpy.ndarray, optional
        Precomputed Mel filter banks
    fs: int
        Sampling frequency of the signal from which the power spectrogram was computed. Ignored if precomputed Mel
        filter banks are used.
    window_width: int
        Window width in samples that was used to comput the power spectrogram. Ignored if precomputed Mel filter banks 
        are used.
    n_filt: int
        Number of Mel filter banks to use. Ignored if precomputed Mel filter banks are used.

    Returns
    -------
    numpy.ndarray
        Mel spectrogram computed from the specified power spectrogram
    """
    if mel_fbank is None:
        _, mel_fbank = mel_filter_bank(fs, window_width, n_filt)

    filter_banks = np.dot(mel_fbank, power_spectrum)

    return filter_banks 

def power_to_db(spectrum: np.ndarray,
                clip_below: float = None,
                clip_above: float = None) -> np.ndarray:
    """
    Convert a spectrogram to the Decibel scale.
    
    Optionally, frequencies with amplitudes below or above a certain threshold can be clipped.
    
    Parameters
    ----------
    spectrum: numpy.ndarray
        The spectrogram to convert
    clip_below: float, optional
        Clip frequencies below the specified amplitude in dB
    clip_above: float, optional
        Clip frequencies above the specified amplitude in dB

    Returns
    -------
    numpy.ndarray
        The spectrogram on the Decibel scale
    """
    # there might be zeros, fix them to the lowest non-zero power in the spectrogram
    epsilon = np.min(spectrum[np.where(spectrum > 0)])

    sxx = np.where(spectrum > epsilon, spectrum, epsilon)
    sxx = 10 * np.log10(sxx / np.max(sxx))

    if clip_below is not None:
        sxx = np.maximum(sxx, clip_below)

    if clip_above is not None:
        sxx = np.minimum(sxx, clip_above)

    return sxx 

def compute_scv(sig, fs, window='hann', nperseg=None, noverlap=0, outlier_pct=None):
    """Compute the spectral coefficient of variation (SCV) at each frequency.

    Parameters
    -----------
    sig : 1d array
        Time series of measurement values.
    fs : float
        Sampling rate, in Hz.
    window : str or tuple or array_like, optional, default: 'hann'
        Desired window to use. See scipy.signal.get_window for a list of available windows.
        If array_like, the array will be used as the window and its length must be nperseg.
    nperseg : int, optional
        Length of each segment, in number of samples.
        If None, and window is str or tuple, is set to 1 second of data.
        If None, and window is array_like, is set to the length of the window.
    noverlap : int, optional, default: 0
        Number of points to overlap between segments.
    outlier_pct : float, optional
        Percentage of the windows with the lowest and highest total log power to discard.
        Must be between 0 and 100.

    Returns
    -------
    freqs : 1d array
        Frequencies at which the measure was calculated.
    scv : 1d array
        Spectral coefficient of variation.

    Notes
    -----
    White noise should have a SCV of 1 at all frequencies.

    Examples
    --------
    Compute the spectral coefficient of variation of a simulated time series:

    >>> from neurodsp.sim import sim_combined
    >>> sig = sim_combined(n_seconds=10, fs=500,
    ...                    components={'sim_powerlaw': {}, 'sim_oscillation' : {'freq': 10}})
    >>> freqs, scv = compute_scv(sig, fs=500)
    """

    # Compute spectrogram of data
    nperseg, noverlap = check_spg_settings(fs, window, nperseg, noverlap)
    freqs, _, spg = spectrogram(sig, fs, window, nperseg, noverlap)

    if outlier_pct is not None:
        spg = discard_outliers(spg, outlier_pct)

    scv = np.std(spg, axis=-1) / np.mean(spg, axis=-1)

    return freqs, scv 

def computeKey(afAudioData, f_s, afWindow=None, iBlockLength=4096, iHopLength=2048):

    # compute window function for FFT
    if afWindow is None:
        afWindow = ToolComputeHann(iBlockLength)

    assert(afWindow.shape[0] == iBlockLength), "parameter error: invalid window dimension"

    # key names
    cKeyNames = np.array(['C Maj', 'C# Maj', 'D Maj', 'D# Maj', 'E Maj', 'F Maj', 'F# Maj', 'G Maj', 'G# Maj', 'A Maj', 'A# Maj', 'B Maj',
                         'c min', 'c# min', 'd min', 'd# min', 'e min', 'f min', 'f# min', 'g min', 'g# min', 'a min', 'a# min', 'b min'])

    # template pitch chroma (Krumhansl major/minor), normalized to a sum of 1
    t_pc = np.array([[6.35, 2.23, 3.48, 2.33, 4.38, 4.09, 2.52, 5.19, 2.39, 3.66, 2.29, 2.88],
                    [6.33, 2.68, 3.52, 5.38, 2.60, 3.53, 2.54, 4.75, 3.98, 2.69, 3.34, 3.17]])
    t_pc = t_pc / t_pc.sum(axis=1, keepdims=True)

    # pre-processing
    afAudioData = ToolPreprocAudio(afAudioData, iBlockLength)

    # in the real world, we would do this block by block...
    [f, t, X] = spectrogram(afAudioData,
                            f_s,
                            afWindow,
                            iBlockLength,
                            iBlockLength - iHopLength,
                            iBlockLength,
                            False,
                            True,
                            'spectrum')

    #  scale the same as for matlab
    X = np.sqrt(X / 2)

    # compute instantaneous pitch chroma
    v_pc = FeatureSpectralPitchChroma(X, f_s)

    # average pitch chroma
    v_pc = v_pc.mean(axis=1)
    # compute manhattan distances for modes (major and minor)
    d = np.zeros(t_pc.shape)
    v_pc = np.concatenate((v_pc, v_pc), axis=0).reshape(2, 12)
    for i in range(0, 12):
        d[:, i] = np.sum(np.abs(v_pc - np.roll(t_pc, i, axis=1)), axis=1)

    # get unwrapped key index
    iKeyIdx = d.argmin()

    cKey = cKeyNames[iKeyIdx]

    return (cKey) 

def computeNoveltyFunction(cNoveltyName, afAudioData, f_s, afWindow=None, iBlockLength=4096, iHopLength=512):

    # compute window function for FFT
    if afWindow is None:
        afWindow = ToolComputeHann(iBlockLength)

    assert(afWindow.shape[0] == iBlockLength), "parameter error: invalid window dimension"

    #mypackage = __import__(".Novelty" + cNoveltyName, package="pyACA")
    hNoveltyFunc = getattr(pyACA, "Novelty" + cNoveltyName)

    # initialization
    fLengthLpInS = 0.3
    iLengthLp = np.max([2, math.ceil(fLengthLpInS * f_s / iHopLength)])

    # pre-processing
    afAudioData = ToolPreprocAudio(afAudioData, iBlockLength)

    # in the real world, we would do this block by block...
    [f, t, X] = spectrogram(afAudioData,
                            f_s,
                            afWindow,
                            iBlockLength,
                            iBlockLength - iHopLength,
                            iBlockLength,
                            False,
                            True,
                            'spectrum')

    #  scale the same as for matlab
    X = np.sqrt(X / 2)

    # novelty function
    d = hNoveltyFunc(X, f_s)

    # smooth novelty function
    b = np.ones(10) / 10
    d = filtfilt(b, 1, d)
    d[d < 0] = 0

    # compute threshold
    b = np.ones(iLengthLp) / iLengthLp
    G_T = .5 * np.mean(d[np.arange(1, d.shape[0])]) + filtfilt(b, 1, d)

    # find local maxima above the threshold
    iPeaks = find_peaks(d - G_T, height=0)

    return (d, t, iPeaks[0]) 

def aggregate(self, ldframe, **kwargs):
        """Run a collection of annotators over the input material.

        If output_dir is not none, produces PNG files of the spectrograms for
        each group in the desired output location. If spectrogram is set to
        True, will return the numeric spectrograms. Otherwise returns an
        empty output.
        """
        _check_data_exists(ldframe, ["meta", "audio", "audiometa"])

        if self.output_dir is not None:
            _check_out_dir(self.output_dir)
            use("template")

        dta = ldframe['audio']['data'].values
        rate = ldframe['audiometa']['rate'].values[0]

        saved_times = []
        saved_specs = []

        for stime, audio, i in _audio_chunks(
            self.breaks, dta, rate, ldframe['meta']['fps']
        ):

            frequencies, times, spec = spectrogram(audio, fs=rate)

            if self.output_dir is not None:
                opath = join(self.output_dir, "frame-{0:06d}.png".format(i))

                pcolormesh(times + int(stime), frequencies, 10 * log10(spec))
                xlabel("Time (seconds)")
                ylabel("Frequency")
                savefig(opath)
                close()

            if self.spectrogram:
                saved_times.extend(times + stime)
                saved_specs.extend([transpose(spec)])
                print(spec.shape)

        if self.spectrogram:
            return {
                'times': saved_times,
                'spectrogram': vstack(saved_specs)
            }

        return None 

def mel_filter_bank(fs: int,
                    window_width: int,
                    n_filt: int = 40) -> (np.ndarray, np.ndarray):
    """
    Computes Mel filter banks for the specified parameters.
    
    A power spectrogram can be converted to a Mel spectrogram by multiplying it with the filter bank. This method exists
    so that the computation of Mel filter banks does not have to be repeated for each computation of a Mel spectrogram.
    
    The coefficients of Mel filter banks depend on the sampling frequency and the window width that were used to 
    generate power spectrograms.
    
    Parameters
    ----------
    fs: int
        The sampling frequency of signals
    window_width: int
        The window width in samples used to generate spectrograms
    n_filt: int
        The number of filters to compute

    Returns
    -------
    f: numpy.ndarray
        Array of Hertz frequency values for the filter banks
    filters: numpy.ndarray
        Array of Mel filter bank coefficients. The first axis corresponds to different filters, and the second axis
        corresponds to the original frequency bands
    """
    n_fft = window_width
    low_freq_mel = 0
    high_freq_mel = (2595 * np.log10(1 + (fs / 2) / 700))  # Convert Hz to Mel
    mel_points = np.linspace(low_freq_mel, high_freq_mel, n_filt + 2)  # Equally spaced in Mel scale
    hz_points = (700 * (10 ** (mel_points / 2595) - 1))  # Convert Mel to Hz
    bin = np.floor((n_fft + 1) * hz_points / fs)

    fbank = np.zeros((n_filt, int(np.floor(n_fft / 2 + 1))))

    for m in range(1, n_filt + 1):
        f_m_minus = int(bin[m - 1])  # left
        f_m = int(bin[m])  # center
        f_m_plus = int(bin[m + 1])  # right

        for k in range(f_m_minus, f_m):
            fbank[m - 1, k] = (k - bin[m - 1]) / (bin[m] - bin[m - 1])
        for k in range(f_m, f_m_plus):
            fbank[m - 1, k] = (bin[m + 1] - k) / (bin[m + 1] - bin[m])

    return hz_points[1:n_filt + 1], fbank