Blog from: https://www.cnblogs.com/LXP-Never/p/10886622.html
Technical specifications of filter
- w p w_p wp: passband cut-off frequency
- w s w_s ws: stopband cut-off frequency
- δ p δ_p δ p: passband fluctuation
- δ s δ_s δ s: stopband fluctuation
Relationship among filter coefficient, impulse response and frequency response: filter transformation, time-domain convolution is equal to frequency-domain product, and filtering operation is represented as convolution of input signal and filter impulse response in time-domain; In the frequency domain, the filter operation is expressed as the product of the Fourier transform of the input signal and the Fourier transform of the impulse response.
Low pass filter: only signals below a certain frequency are allowed to pass through the filter without attenuation, and the frequency at its boundary is called cut-off frequency.
picture
High pass filter: it is specified that the signal above the set critical frequency (fc) can pass normally, while the signal below the set critical frequency (fc) is blocked and attenuated.
Bandpass filter: a filter that allows a specific frequency signal to pass through, blocking and attenuating the signals at the upper and lower frequencies of the frequency band. The band-pass filter can be regarded as the result of the synergistic action of high pass and low-pass filters. The cut-off frequencies of high pass and low-pass filters can be used as the lower limit frequency and upper limit frequency of the passband (points L and U in Fig. 3). The main parameters are the center frequency and bandwidth. The difference between the upper limit frequency and the lower limit frequency is the bandwidth of the filter.
Band stop filter: a filter that can pass through most frequency components but attenuate some frequency components to a very low level, which is opposite to the concept of band-pass filter. notch filter is a special band stop filter, and its stop band range is very small.
Notch filter: it can quickly attenuate the input signal at a certain frequency point to achieve the filtering effect of hindering the passage of this frequency signal. It is mainly used to filter the signal of unnecessary frequency on the circuit.
For FIR filter, the filter coefficient is the impulse response. Therefore, for FIR filter, the FFT transformation of the coefficient is the frequency response curve of the filter
Butterworth filter
Frequency domain characteristics of butterworth low pass filter
N
N
N: Order of filter
ω
c
ω_c
ω c: 3dB cutoff
Python implementation
scipy.signal.butter(N, Wn, btype='low', analog=False, output='ba')
input
- N N N: Order of filter
-
W
n
W_n
Wn: cut off frequency. For Butterworth filter, this is the gain falling to the passband
1
2
\frac{1}{\sqrt{2}}
2
1 - point (i.e. "- 3db point").
- For digital filters: W n W_n Wn , and f s f_s The units of fs are the same H z H_z Hz (fs is the sampling rate), in MATLAB W n W_n Wn , is the normalized cut-off frequency (0,1), W n W_n Wn = cut-off frequency / signal frequency, (signal frequency = half of the sampling rate, Nyquist sampling theorem). Normalized frequency = 1, cut-off frequency w c w_c wc represents sampling frequency and 0.5 represents Nyquist frequency
- For analog filters: W n W_n Wn is the angular frequency (radians/second)
- btype: filter type {'lowpass',' highpass', 'bandpass',' bandstop '}
- Analog: if True, the analog filter is returned; otherwise, the digital filter is returned
output
- b. A: filter coefficient, a is denominator and B is numerator.
scipy.signal.freqs(b, a, worN=200, plot=None)
Calculate the frequency response H(w) of the analog filter.
parameter
- b. A: numerator and denominator of linear filter,
- Word: optional. If it is None, 200 frequencies around the interesting part of the response curve are calculated. If it is an integer, so many frequencies are calculated.
return
- w: Calculate the angular frequency of h
- h: Frequency response (amplitude at this frequency)
scipy.signal.freqz(b, a=1, worN=512, whole=False, plot=None)
Calculate the frequency response of the digital filter.
parameter
- b. A: numerator and denominator of linear filter
- Word: if it is None (default), 512 frequencies with equal intervals around the unit circle are calculated. If it is an integer, so many frequencies are calculated. If it's array_like, the response (in radians / sample) of a given frequency is calculated.
return
- w: Calculate the frequency of h in the same unit as fs. By default, W is normalized to [0,1) * π [0,1) * π range (radians/sample). For example, for a system with a sampling frequency of 1000Hz, 300Hz corresponds to 300 / 500 = 0.6. To convert the normalized frequency to radians on the unit circle, multiply the normalized value by π. 0.5 in the figure represents 0.5 π (rad) , that is, w = 0.5 π. From w = 2 * pi * f / fs, f = w * fs / 2pi, that is, f = 0.5 * fs / 2. Because fs = 122.88MHz, the cut-off frequency f = 30.72MHz.
- h: Frequency response, expressed in complex number
scipy.signal.lfilter(b, a, x, axis=-1, zi=None)
Use IIR or FIR filters to filter the data along one dimension. Use digital filters to filter the data sequence x.
input
- b. A: numerator and denominator, i.e. filter coefficient
- x: Input data
Return: output of digital filter
from scipy.signal import butter, lfilter from scipy import signal import numpy as np import matplotlib.pyplot as plt b, a = signal.butter(4, 100, 'low', analog=True) # Design an N-order digital or analog Butterworth filter and return the filter coefficients w, h = signal.freqs(b, a) # The frequency response of the filter is calculated according to the coefficients. w is the angular frequency and h is the frequency response plt.semilogx(w, 20 * np.log10(abs(h))) plt.title('Butterworth filter frequency response') plt.xlabel('Frequency [radians / second]') plt.ylabel('Amplitude [dB]') plt.margins(0, 0.1) plt.grid(which='both', axis='both') plt.axvline(100, color='green') # cutoff frequency plt.show()
Extract narrowband speech signal
For the speech with sampling rate of 16000Hz and Nyquist frequency of 8000Hz, the speech with frequency higher than 4000Hz is filtered through Butterworth low-pass filter to extract low-frequency speech. When the sampling rate is the same, the frequency of the filtered signal is only half of the original.
import librosa import numpy as np from scipy.signal import butter, lfilter, freqz import matplotlib.pyplot as plt def butter_lowpass(cutoff, fs, order=5): # Cutoff: cutoff frequency # fs sampling rate nyq = 0.5 * fs # signal frequency normal_cutoff = cutoff / nyq # Normal cut-off frequency = cut-off frequency / signal frequency b, a = butter(order, normal_cutoff, btype='lowpass', analog=False) return b, a def butter_lowpass_filter(data, cutoff, fs, order=5): b, a = butter_lowpass(cutoff, fs, order=order) y = lfilter(b, a, data) return y # Filter requirements. order = 10 fs = 16000 # Sampling rate, Hz cutoff = 4000 # The desired cut-off frequency of the filter, Hz # Get the filter coefficients so that we can check its frequency response. b, a = butter_lowpass(cutoff, fs, order) # Plot frequency response w, h = freqz(b, a) plt.subplot(3, 1, 1) plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b') plt.axvline(cutoff, color='k') plt.xlim(0, 0.5*fs) plt.title("Lowpass Filter Frequency Response") plt.xlabel('Frequency [Hz]') data, wav_fs = librosa.load("./48k/p225_001.wav", sr=16000, mono=True) # 48000--->16000 y = butter_lowpass_filter(data, cutoff, fs, order) plt.subplot(3, 1, 2) plt.specgram(data, Fs=16000, scale_by_freq=True, sides='default') plt.xlabel('Time [sec]') plt.subplot(3, 1, 3) plt.specgram(y, Fs=16000, scale_by_freq=True, sides='default') plt.xlabel('Time [sec]') plt.show()
Chebyshev I shape filter
Frequency domain characteristics of CB I low pass filter
N
N
N: Order of filter
ε ε ε: Passband ripple
ω
c
ω_c
ω c: passband cutoff
Features: the passband is equal fluctuation, and the stopband is monotonous
scipy.signal.cheby1(N, rp, Wn, btype='low', analog=False, output='ba')
Chebyshev type I digital and analog filters, design N-order digital or analog Chebyshev type I filters and return filter coefficients
Parameters:
- N: Order of filter
- rp: the maximum ripple allowed in the passband is lower than the unit gain, in decibels, positive
- Wn: for digital filters: wn should be normalized to (0,1), WN = cut-off frequency / signal frequency, (signal frequency = half of the sampling rate, Nyquist sampling theorem) for analog filters: wn is angular frequency, radian / sample, rad/s
- btype: filter type {'lowpass',' highpass', 'bandpass',' bandstop '}
- Analog: if True, the analog filter is returned; otherwise, the digital filter is returned.
- output: the default is "ba", which outputs numerator and denominator
return
b. A: filter coefficient, a is denominator and B is numerator.
import numpy as np from scipy import signal import matplotlib.pyplot as plt b, a = signal.cheby1(4, 5, 100, 'low', analog=True) w, h = signal.freqs(b, a) plt.semilogx(w, 20 * np.log10(abs(h))) plt.title('Chebyshev Type I frequency response (rp=5)') plt.xlabel('Frequency [radians / second]') plt.ylabel('Amplitude [dB]') plt.margins(0, 0.1) plt.grid(which='both', axis='both') plt.axvline(100, color='green') # cutoff frequency plt.axhline(-5, color='green') # rp plt.show()
Chebyshev II shape filter
Frequency domain characteristics of CB II low pass filter
N
N
N: Order of filter
ε
ε
ε: Stopband ripple
ω
c
ω_c
ω c: stopband cutoff
Features: the passband is monotonic, and the stopband is equifluctuating
scipy.signal.cheby2(N, rs, Wn, btype='low', analog=False, output='ba')
Chebyshev II Digital and analog filters, design N-order digital or analog Chebyshev II filters and return filter coefficients.
Parameters:
- N: Order of filter
- rs: minimum attenuation required for stopband, in decibels, positive
- Wn: for digital filters: wn should be normalized to (0,1), WN = cut-off frequency / signal frequency, (signal frequency = half of sampling rate, Nyquist sampling theorem) for analog filters: wn is angular frequency, radian / sample, rad/s
- btype: filter type {'lowpass',' highpass', 'bandpass',' bandstop '}
- Analog: if True, the analog filter is returned; otherwise, the digital filter is returned.
- output: the default is "ba", which outputs numerator and denominator
return
- b. A: filter coefficient, a is denominator and B is numerator.
rom scipy import signal import numpy as np import matplotlib.pyplot as plt b, a = signal.cheby2(4, 40, 100, 'low', analog=True) w, h = signal.freqs(b, a) plt.semilogx(w, 20 * np.log10(abs(h))) plt.title('Chebyshev Type II frequency response (rs=40)') plt.xlabel('Frequency [radians / second]') plt.ylabel('Amplitude [dB]') plt.margins(0, 0.1) plt.grid(which='both', axis='both') plt.axvline(100, color='green') # cutoff frequency plt.axhline(-40, color='green') # rs plt.show()
Elliptic low pass filter
Frequency domain characteristics of elliptic analog low-pass filter
Features: both passband and stopband fluctuate equally
scipy.signal.ellip(N, rp, rs, Wn, btype='low', analog=False, output='ba')
Elliptic digital and analog filters, design N-order digital or analog elliptic filters and return filter coefficients.
Parameters:
- N: Order of filter
- rp: the maximum ripple allowed in the passband is lower than the unit gain, in decibels, positive
- rs: minimum attenuation required for stopband, in decibels, positive
- Wn: for digital filters: wn should be normalized to (0,1), WN = cut-off frequency / signal frequency, (signal frequency = half of sampling rate, Nyquist sampling theorem) for analog filters: wn is angular frequency, radian / sample, rad/s
- btype: filter type {'lowpass',' highpass', 'bandpass',' bandstop '}
- Analog: if True, the analog filter is returned; otherwise, the digital filter is returned.
- output: the default is "ba", which outputs numerator and denominator
return
- b. A: filter coefficient, a is denominator and B is numerator.
import numpy as np from scipy import signal import matplotlib.pyplot as plt b, a = signal.ellip(4, 5, 40, 100, 'low', analog=True) w, h = signal.freqs(b, a) plt.semilogx(w, 20 * np.log10(abs(h))) plt.title('Elliptic filter frequency response (rp=5, rs=40)') plt.xlabel('Frequency [radians / second]') plt.ylabel('Amplitude [dB]') plt.margins(0, 0.1) plt.grid(which='both', axis='both') plt.axvline(100, color='green') # cutoff frequency plt.axhline(-40, color='green') # rs plt.axhline(-5, color='green') # rp plt.show()
Down sampling method
Interpolation method for down sampling
In Volodymyr Kuleshov's paper, the anti aliasing filter is used to down sample the speech signal, and then the down sampled signal is up sampled to the same length through cubic spline interpolation
from scipy.signal import decimate import librosa import numpy as np import matplotlib.pyplot as plt from scipy import interpolate def upsample(x_lr, r): """ Up sampling, remove the noise of voice waveform every other step r Three points, and then use the cubic spline interpolation method to make up the removed points. If you have a chance, you can draw a picture :param x_lr: Audio data :param r: Number before spline interpolation :return: Audio signal after spline interpolation """ x_lr = x_lr.flatten() # Put x_lr arrays are folded into one-dimensional arrays x_hr_len = len(x_lr) * r i_lr = np.arange(x_hr_len, step=r) i_hr = np.arange(x_hr_len) f = interpolate.splrep(i_lr, x_lr) # Spline interpolation coefficient x_sp = interpolate.splev(i_hr, f) # Given the nodes and coefficients represented by the spline, return the spline value at the node return x_sp yt, wav_fs = librosa.load("./48k/p225_001.wav", sr=16000, mono=True) x_lr = decimate(yt, 2) # After the anti aliasing filter is applied, the signal is down sampled to obtain low-resolution audio, and the down sampling factor is scale=2 print(len(yt)) print(len(x_lr)) plt.subplot(2, 1, 1) plt.specgram(yt, Fs=16000, scale_by_freq=True, sides='default') x_lr = upsample(x_lr, 2) # Up sampling plt.subplot(2, 1, 2) plt.specgram(x_lr, Fs=16000, scale_by_freq=True, sides='default') plt.show()
Resampling (signal.resample) -- down sampling
Resampling method is used to down sample speech
scipy.signal.resample(x, num, t=None, axis=0, window=None)
Resample x to num samples along a given axis using the Fourier method. Because the Fourier method is used, it is assumed that the signal is periodic.
Parameters:
- x: Array to resample
- num: number of samples of the resampled signal
return: - resample_x: Resampling returned array
import librosa from scipy import signal import numpy as np import matplotlib.pyplot as plt y, wav_fs = librosa.load("./48k/p225_001.wav", sr=16000, mono=True) f = signal.resample(y, len(y)//2) f = signal.resample(f, len(y)) plt.subplot(2,1,1) plt.specgram(y, Fs=16000, scale_by_freq=True, sides='default') plt.subplot(2,1,2) plt.specgram(f, Fs=16000, scale_by_freq=True, sides='default') plt.show()
librosa.core.resample resampling (down sampling)
Teacher Ling Zhenhua's down sampling method is the same as that above
import librosa from scipy import signal import numpy as np import matplotlib.pyplot as plt y, wav_fs = librosa.load("./48k/p225_001.wav", sr=16000, mono=True) audio8k = librosa.core.resample(y, wav_fs, wav_fs/2) # Lower sampling rate 16000 -- > 8000 audio8k = librosa.core.resample(audio8k, wav_fs/2, wav_fs) # The upper sampling rate is 8000 -- > 16000, and the high-frequency part is not recovered plt.subplot(2,1,1) plt.specgram(y, Fs=16000, scale_by_freq=True, sides='default') plt.subplot(2,1,2) plt.specgram(audio8k, Fs=16000, scale_by_freq=True, sides='default') plt.show()
Sampling under librosa.load
Use librosa.load to down sample, and then up sample without restoring the frequency.
import librosa import matplotlib.pyplot as plt y_16k, fs_16k = librosa.load("./48k/p225_001.wav", sr=16000, mono=True) y_8k, fs_8k = librosa.load("./48k/p225_001.wav", sr=8000, mono=True) librosa.output.write_wav('./8k_sample.wav', y_8k, sr=8000) # Write down the samples y_8k, fs_8k = librosa.load("./8k_sample.wav", sr=16000, mono=True) # You can't make up for what you lose plt.subplot(2, 1, 1) plt.specgram(y_16k, Fs=16000, scale_by_freq=True, sides='default') plt.xlabel('16k') plt.subplot(2, 1, 2) plt.specgram(y_8k, Fs=16000, scale_by_freq=True, sides='default') plt.xlabel('8k') plt.show()
reference
Professor Chen Houjin, Beijing Jiaotong University (digital signal processing)
Signal processing (scipy.signal)
Meaning of Normalized Frequency (* π rad/sample)
Differences between high pass filter, low pass filter, band-pass filter and notch filter