Speech signal preprocessing -- digital filter

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").
1. 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
2. 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()