Source code for maad.sound.filter

#!/usr/bin/env python
""" 
Collection of functions to filter audio signal in time and frequency
domains. 
"""
#
# Authors:  Juan Sebastian ULLOA <lisofomia@gmail.com>
#           Sylvain HAUPERT <sylvain.haupert@mnhn.fr>
#
# License: New BSD License

# =============================================================================
# Load the modules
# =============================================================================
# Import external modules
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import (
    sosfiltfilt, convolve, iirfilter, get_window, 
    kaiserord, firwin, lfilter
    )
from skimage import filters

# Import internal modules 
from maad.util import plot2d

#%%
# =============================================================================
# public functions
# =============================================================================
[docs] def select_bandwidth (x, fs, fcut, forder, fname ='butter', ftype='bandpass', rp=None, rs=None): """ Use a lowpass, highpass, bandpass or bandstop filter to process a 1d signal with an iir filter. Parameters ---------- x : array_like 1d vector of scalars to be filtered fs : scalar sampling frequency fcut : array_like A scalar or length-2 sequence giving the critical frequencies. forder : int The order of the filter. fname : {`butter`, 'cheby1`, `cheby2`, `ellip`, `bessel`}, optional, default is 'butter' The type of IIR filter to design: - Butterworth : `butter` - Chebyshev I : 'cheby1` - Chebyshev II : `cheby2` - Cauer/elliptic: `ellip` - Bessel/Thomson: `bessel` ftype : {`bandpass`, `lowpass`, `highpass`, `bandstop`}, optional, default is `bandpass` The type of filter. rp : float, optional For Chebyshev and elliptic filters, provides the maximum ripple in the passband. (dB) rs : float, optional For Chebyshev and elliptic filters, provides the minimum attenuation in the stop band. (dB) Returns ------- y : array_like The filtered output with the same shape and phase as x See Also -------- fir_filter1d Lowpass, highpass, bandpass or bandstop a 1d signal with an Fir filter Examples -------- >>> import maad Load and display the spectrogram of a sound waveform >>> w, fs = maad.sound.load('../data/cold_forest_daylight.wav') >>> Sxx_power,tn,fn,_ = maad.sound.spectrogram(w,fs) >>> Sxx_dB = maad.spl.power2dBSPL(Sxx_power, gain=42) # convert into dB SPL >>> fig_kwargs = {'vmax': Sxx_dB.max(),\ 'vmin':0,\ 'extent':(tn[0], tn[-1], fn[0], fn[-1]),\ 'figsize':(4,13),\ 'title':'Power spectrogram density (PSD)',\ 'xlabel':'Time [sec]',\ 'ylabel':'Frequency [Hz]',\ } >>> ax, fig = maad.util.plot2d(Sxx_dB, **fig_kwargs) >>> print('The energy in the spectrogram is {} dB SPL'.format(maad.util.add_dB(maad.util.add_dB(Sxx_dB)))) The energy in the spectrogram is 90.68996875800728 dB SPL Filter the waveform : keep the bandwidth between 6-10kHz >>> w_filtered = maad.sound.select_bandwidth(w,fs,fcut=(6000,10000), forder=5, fname ='butter', ftype='bandpass') >>> Sxx_power_filtered,tn,fn,_ = maad.sound.spectrogram(w_filtered,fs) >>> Sxx_dB_filtered = maad.spl.power2dBSPL(Sxx_power_filtered, gain=42) # convert into dB SPL >>> ax, fig = maad.util.plot2d(Sxx_dB_filtered, **fig_kwargs) >>> print('The energy in the spectrogram is {} dB SPL'.format(maad.util.add_dB(maad.util.add_dB(Sxx_dB_filtered)))) The energy in the spectrogram is 72.75045485363799 dB SPL """ sos = iirfilter(N=forder, Wn=np.asarray(fcut)/(fs/2), btype=ftype,ftype=fname, rp=rp, rs=rs, output='sos') # use sosfiltfilt insteasd of sosfilt to keep the phase of y matches x y = sosfiltfilt(sos, x) return y
#%%
[docs] def fir_filter(x, kernel, axis=0): """ Filter a signal using a 1d finite impulse response filter. This function uses a digital filter based on convolution of 1d kernel over a vector or along an axis of a matrix. Parameters ---------- x : array_like 1d vector or 2d matrix of scalars to be filtered kernel : array_like or tuple Pass directly the kernel (1d vector of scalars) Or pass the arguments in a tuple to create a kernel. Arguments are: - window : string, float, or tuple. The type of window to create. boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann, - (kaiser, beta), - (gaussian, standard deviation), - (general_gaussian, power, width), - (slepian, width), - (dpss, normalized half-bandwidth), - (chebwin, attenuation), - (exponential, decay scale), - (tukey, taper fraction) - N : length of the kernel Examples: - kernel = ('boxcar', 9) - kernel = (('gaussian', 0.5), 5) - kernel = [1 3 5 7 5 3 1] axis : int Determine along which axis is performed the filtering in case of 2d matrix axis = 0 : vertical axis = 1 : horizontal Returns ------- y : array_like The filtered output with the same shape and phase as x See Also -------- select_bandwidth Lowpass, highpass, bandpass or bandstop a 1d signal with an iir filter Examples -------- >>> import maad Load and display the spectrogram of a sound waveform >>> w, fs = maad.sound.load('../data/cold_forest_daylight.wav') >>> Sxx_power, tn, fn, ext = maad.sound.spectrogram(w,fs) >>> Lxx = maad.spl.power2dBSPL(Sxx_power, gain=42) # convert into dB SPL >>> fig_kwargs = {'vmax': Lxx.max(),\ 'vmin':0,\ 'extent': ext,\ 'figsize': (4,13),\ 'title': 'Power spectrogram density (PSD)',\ 'xlabel': 'Time [sec]',\ 'ylabel': 'Frequency [Hz]',\ } >>> ax, fig = maad.util.plot2d(Lxx,**fig_kwargs) >>> print('The energy in the spectrogram is {} dB SPL'.format(maad.util.add_dB(maad.util.add_dB(Lxx)))) The energy in the spectrogram is 90.68996875800728 dB SPL Smooth the waveform (lowpass) >>> w_filtered = maad.sound.fir_filter(w, kernel=(('gaussian', 2), 5)) >>> Sxx_power_filtered,tn,fn,_ = maad.sound.spectrogram(w_filtered,fs) >>> Lxx_filtered = maad.spl.power2dBSPL(Sxx_power_filtered, gain=42) # convert into dB SPL >>> ax, fig = maad.util.plot2d(Lxx_filtered,**fig_kwargs) >>> print('The energy in the spectrogram is {} dB SPL'.format(maad.util.add_dB(maad.util.add_dB(Lxx_filtered)))) The energy in the spectrogram is 89.51776241848223 dB SPL Smooth the spectrogram, frequency by frequency (blurr) >>> Lxx_blurr = maad.sound.fir_filter(Lxx, kernel=(('gaussian', 1), 5), axis=1) >>> ax, fig = maad.util.plot2d(Lxx_blurr,**fig_kwargs) >>> print('The energy in the spectrogram is {} dB SPL'.format(maad.util.add_dB(maad.util.add_dB(Lxx_blurr)))) The energy in the spectrogram is 88.83611122906812 dB SPL """ if isinstance(kernel,tuple) : if len(kernel) ==1 : win = get_window(kernel[0]) if len(kernel) ==2 : win = get_window(kernel[0], kernel[1]) elif isinstance(kernel,list) or isinstance(kernel,np.ndarray): win = kernel if x.ndim == 1: # Trick to avoid strange values at the beginning and end of the convolve # signal. Add mirror values x = np.insert(x, 0, np.flipud(x[0:len(win)//2]), axis=0) x = np.insert(x, len(x), np.flipud(x[-(len(win)//2):]), axis=0) # convolve and normalized the result by the sum of the kernel y = convolve(x, win, mode='same') / sum(win) y = y[len(win)//2:-(len(win)//2)] elif x.ndim == 2: if axis ==1 : x = x.transpose() # Trick to avoid strange values at the beginning and end of the convolve # signal. Add mirror values x = np.insert(x, 0, np.flipud(x[0:len(win)//2]), axis=0) x = np.insert(x, x.shape[0], np.flipud(x[-(len(win)//2):]), axis=0) y = np.zeros(x.shape) for i in np.arange (x.shape[1]): y[:,i] = convolve(x[:,i], win, mode='same') / sum(win) y = y[len(win)//2:-(len(win)//2),:] if axis ==1 :y = y.transpose() return y
#%%
[docs] def sinc(s, cutoff, fs, atten=80, transition_bw=0.05, bandpass=True): """ Filter 1D signal with a Kaiser-windowed filter. Parameters ---------- s : ndarray input 1D signal cutoff : ndarray upper and lower frequencies (min_f, max_f) in Hertz atten : float attenuation in dB transition_bw : float transition bandwidth in percent default 5% of total band bandpass : bool bandpass (True) or bandreject (False) filter, default is bandpass Returns ------- s_filt : array signal filtered Examples -------- >>> import maad >>> s, fs = maad.sound.load('../data/spinetail.wav') >>> import numpy as np >>> tn = np.arange(0,len(s))/fs >>> s_filt_7_12kHz = maad.sound.sinc(s, cutoff=[7000,12000], fs=fs, atten=80, transition_bw=0.8) >>> s_filt_4_8kHz = maad.sound.sinc(s, cutoff=[4500,8000], fs=fs, atten=80, transition_bw=0.8) >>> import matplotlib.pyplot as plt >>> fig, (ax0, ax1,ax2) = plt.subplots(3,1, sharex=True, squeeze=True) >>> ax0, fig = maad.util.plot1d(tn,s,ax=ax0, figtitle='original') >>> ax1, fig = maad.util.plot1d(tn,s_filt_7_12kHz,ax=ax1, figtitle='Kaiser-windowed filter 7-12kHz') >>> a2 , fig = maad.util.plot1d(tn,s_filt_4_8kHz,ax=ax2, figtitle='Kaiser-windowed filter 4.5-8kHz') >>> fig.tight_layout() """ width = (cutoff[1] - cutoff[0]) * transition_bw numtaps, beta = kaiserord(atten, width/(0.5*fs)) np.ceil(numtaps-1) // 2 * 2 + 1 # round to nearest odd to have Type I filter taps = firwin( numtaps, cutoff, window=('kaiser', beta), scale=False, nyq=0.5*fs, pass_zero=not(bandpass) ) s_filt = lfilter(taps, 1, s) return s_filt
#%%
[docs] def smooth (Sxx, std: float=1, verbose=False, display = False, savefig=None, **kwargs): """ Smooth a spectrogram with a gaussian filter. Parameters ---------- Sxx : 2d ndarray of scalars Spectrogram (or image) std : scalar, optional, default is 1 Standard deviation of the gaussian kernel used to smooth the spectrogram. The larger is the number, the smoother will be the image and the longer it takes. Standard values should fall between 0.5 and 3. verbose : boolean, optional, default is False Print messages display : boolean, optional, default is False Display the signal if True savefig : string, optional, default is None Root filename (with full path) is required to save the figures. Postfix is added to the root filename. kwargs, optional. This parameter is used by plt.plot and savefig functions - savefilename : str, optional, default :'_spectro_after_noise_subtraction.png' Postfix of the figure filename - extent : list of scalars [left, right, bottom, top], optional, default: None The location, in data-coordinates, of the lower-left and upper-right corners. If `None`, the image is positioned such that the pixel centers fall on zero-based (row, column) indices. - title : string, optional, default : 'Spectrogram' title of the figure - xlabel : string, optional, default : 'Time [s]' label of the horizontal axis - ylabel : string, optional, default : 'Amplitude [AU]' label of the vertical axis - cmap : string or Colormap object, optional, default is 'gray' See https://matplotlib.org/examples/color/colormaps_reference.html in order to get all the existing colormaps examples: 'hsv', 'hot', 'bone', 'tab20c', 'jet', 'seismic', 'viridis'... - vmin, vmax : scalar, optional, default: None `vmin` and `vmax` are used in conjunction with norm to normalize luminance data. Note if you pass a `norm` instance, your settings for `vmin` and `vmax` will be ignored. - extent : scalars (left, right, bottom, top), optional, default: None The location, in data-coordinates, of the lower-left and upper-right corners. If `None`, the image is positioned such that the pixel centers fall on zero-based (row, column) indices. - dpi : integer, optional, default is 96 Dot per inch. For printed version, choose high dpi (i.e. dpi=300) => slow For screen version, choose low dpi (i.e. dpi=96) => fast - format : string, optional, default is 'png' Format to save the figure ... and more, see matplotlib Returns ------- im_out: smothed or blurred image Examples -------- >>> import maad Load audio recording and convert it into spectrogram >>> s, fs = maad.sound.load('../data/cold_forest_daylight.wav') >>> Sxx,tn,fn,ext = maad.sound.spectrogram (s, fs, tcrop=(5,10), fcrop=(0,10000)) Convert linear spectrogram into dB >>> Sxx_dB = maad.util.power2dB(Sxx) +96 Smooth the spectrogram >>> Sxx_dB_std05 = maad.sound.smooth(Sxx_dB, std=0.5) >>> Sxx_dB_std10 = maad.sound.smooth(Sxx_dB, std=1) >>> Sxx_dB_std15 = maad.sound.smooth(Sxx_dB, std=1.5) Plot spectrograms >>> import matplotlib.pyplot as plt >>> fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1) >>> ax1, fig = maad.util.plot2d(Sxx_dB, ax=ax1, extent=ext, title='original', vmin=10, vmax=70) >>> ax2, fig = maad.util.plot2d(Sxx_dB_std05, ax=ax2, extent=ext, title='smooth (std=0.5)', vmin=10, vmax=70) >>> ax3, fig = maad.util.plot2d(Sxx_dB_std10, ax=ax3, extent=ext, title='smooth (std=1)', vmin=10, vmax=70) >>> ax4, fig = maad.util.plot2d(Sxx_dB_std15, ax=ax4, extent=ext, title='smooth (std=1.5)', vmin=10, vmax=70) >>> fig.set_size_inches(7,9) >>> fig.tight_layout() """ if verbose: print(72 * '_') print('Smooth the image with a gaussian filter (std = %.1f)' %std) # use scikit image (faster than scipy) Sxx_out = filters.gaussian(Sxx,std) # Display if display : ylabel =kwargs.pop('ylabel','Frequency [Hz]') xlabel =kwargs.pop('xlabel','Time [sec]') cmap =kwargs.pop('cmap','gray') vmin=kwargs.pop('vmin',np.percentile(Sxx_out,0.1)) vmax=kwargs.pop('vmax',np.percentile(Sxx_out,99.9)) extent=kwargs.pop('extent',None) if extent is None : ylabel = "pseudofrequency [points]" xlabel = "pseudotime [points]" fig, (ax1, ax2) = plt.subplots(2, 1) plot2d (Sxx, ax=ax1, extent=extent, title=('Orignal Spectrogram'), ylabel = ylabel, xlabel = xlabel, vmin=vmin, vmax=vmax, cmap=cmap, **kwargs) plot2d (Sxx_out, ax=ax2, extent=extent, title='Blurred Spectrogram (std='+str(std)+')', ylabel = ylabel, xlabel = xlabel, vmin=vmin, vmax=vmax, cmap=cmap, **kwargs) # SAVE FIGURE if savefig is not None : dpi =kwargs.pop('dpi',96) format=kwargs.pop('format','png') filename=kwargs.pop('filename','_spectro_blurr') filename = savefig+filename+'.'+format print('\n''save figure : %s' %filename) fig.savefig(filename, bbox_inches='tight', dpi=dpi, format=format, **kwargs) return Sxx_out
if __name__ == "__main__": import doctest doctest.testmod()