濾波器的技術指標
$w_p$：通帶截止頻率
$w_s$：阻帶截止頻率
$\delta_p$：通帶波動
$\delta_s$：阻帶波動
衰減單位是db
濾波器系數、脈沖響應、頻率響應的關系：濾波器變換，時域卷積等於頻域乘積，濾波操作在時域表現為輸入信號與濾波器脈沖響應的卷積；從頻域上看濾波器操作表現為，輸入信號的傅立葉變換和脈沖響應的傅立葉變換做乘積。
低通濾波器：只允許某一頻率以下的信號無衰減地通過濾波器，其分界處的頻率稱為截止頻率。<p id="lowpass" style="#lowpass:hover{}">圖片</p>
高通濾波器：規定高於設定臨界值頻率（fc）的信號能正常通過，而低於設定臨界值頻率(fc)的信號則被阻隔和衰減。
帶通濾波器：允許特定頻率信號通過的濾波器，阻隔和衰減該頻帶上下頻率的信號。帶通濾波器可以看做是由高通和低通濾波器協同作用的結果，高通和低通濾波器的截止頻率可作為通帶的下限頻率和上限頻率（圖3中的L和U點）。其主要參數是中心頻率和帶寬，上限頻率和下限頻率之差就是濾波器的帶寬。
帶阻濾波器：能通過大多數頻率分量，但將某些范圍的頻率分量衰減到極低水平的濾波器，與帶通濾波器的概念相對。其中陷波濾波器（notch filter）是一種特殊的帶阻濾波器，它的阻帶范圍極小。
陷波濾波器：可以在某一個頻率點迅速衰減輸入信號，以達到阻礙此頻率信號通過的濾波效果，主要用在電路上濾除不需要的頻率的信號。
對於FIR濾波器，濾波器系數即為脈沖響應，因此，對於FIR濾波器，系數的FFT變換即為濾波器的頻率響應曲線
巴特沃斯濾波器
butterworth低通濾波器的頻域特性
$$|H(jw)|^2=\frac{1}{1+(\frac{\omega }{\omega _c})^{2N}}$$
N：濾波器的階數
$\omega _c$：3dB截頻
典型BW低通濾波器的幅度響應
特點：在通帶的頻率響應曲線最平滑
Python實現
scipy.signal.butter(N, Wn, btype='low', analog=False, output='ba')
輸入
- N：濾波器的階數
- Wn：截止頻率，對於巴特沃斯濾波器，這是增益下降到通帶$\frac{1}{\sqrt 2}$的點（即“-3db點”）。
- 對於數字濾波器：Wn與fs的單位相同Hz(fs是采樣率)，在MATLAB中Wn是歸一化截止頻率(0,1)，Wn=截止頻率/信號頻率，（信號頻率=采樣率的一半，奈奎斯特采樣定理）。歸一化頻率=1，截止頻率$w_c$代表采樣頻率，0.5代表奈奎斯特頻率
- 對於模擬濾波器：Wn是角頻率（弧度/秒，radians/second）
- btype：濾波器的類型{'lowpass'，'highpass'，'bandpass'，'bandstop'}
- analog：如果為True，則返回模擬濾波器，否則返回數字濾波器。
輸出
- b,a：濾波器系數， a為分母，b為分子。
scipy.signal.freqs(b, a, worN=200, plot=None)
計算模擬濾波器的頻率響應H(w)。
$$H(w)=\frac{b[0]*(jw)**M+b[1]*(jw)**(M-1)+...+b[M]}{ a[0]*(jw)**N + a[1]*(jw)**(N-1) + ... + a[N]}$$
參數
- b,a：線性濾波器的分子和分母，
- worN：可選，如果為None，則計算響應曲線的有趣部分周圍的200個頻率。如果是一個整數，則計算那么多頻率。
返回
- w：計算h的角頻率
- h：頻率響應(在該頻率值得振幅)
scipy.signal.freqz(b, a=1, worN=512, whole=False, plot=None)
計算數字濾波器的頻率響應。
$$H(e^{jw})=\frac{B(e^{jw})}{A(e^{jw})}=\frac{b(1)+b(2)e^{-jw}+...+b(n+1)e^{-jwn}}{a(1)+a(2)e^{-jw}+...+a(m+1)e^{-jwm}}$$
參數
- b,a：線性濾波器的分子和分母
- worN：如果為None（默認值），則計算在單位圓周圍等間隔的512個頻率。如果是一個整數，則計算那么多頻率。如果是array_like，則計算給定頻率的響應（以弧度/樣本為單位）。
返回
- w：計算h的頻率，單位與fs相同。默認情況下，w被歸一化為$[0，1)*\pi$范圍(radians/sample)，例如對於一個采樣頻率為1000Hz的系統，300Hz則對應300/500=0.6。若要將歸一化頻率轉換為單位圓上的弧度，則將歸一化值乘以π即可。圖中的0.5表示的是0.5π(rad)，即為w = 0.5π， 由 w = 2 * pi * f / fs 得到 f = w * fs / 2pi，即 f = 0.5 * fs / 2 ,因為 fs = 122.88MHz，所以截止頻率 f = 30.72MHz。
- h：頻率響應，以復數表示
scipy.signal.lfilter(b, a, x, axis=-1, zi=None)
使用IIR或FIR濾波器沿一維過濾數據。使用數字濾波器過濾數據序列x。
輸入
- b,a：分子和分母，即濾波器系數
- x：輸入數據
返回：數字濾波器的輸出
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) # 設計N階數字或模擬Butterworth濾波器並返回濾波器系數 w, h = signal.freqs(b, a) # 根據系數計算濾波器的頻率響應，w是角頻率，h是頻率響應 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()
提取窄帶語音信號
對采樣率為16000Hz，奈奎斯特頻率為8000Hz的語音，通過巴特沃斯低通濾波器，濾除高於4000Hz頻率的語音，提取低頻語音。過濾出的信號，在采樣率相同的情況下，頻率只有原來的一半。
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：截止頻率 # fs 采樣率 nyq = 0.5 * fs # 信號頻率 normal_cutoff = cutoff / nyq # 正常截止頻率=截止頻率/信號頻率 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 # 采樣率, Hz cutoff = 4000 # 濾波器的期望截止頻率，Hz # 得到濾波器系數，這樣我們就可以檢查它的頻率響應。 b, a = butter_lowpass(cutoff, fs, order) # 繪制頻率響應 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()
切比雪夫I形狀濾波器
CB I型低通濾波器的頻域特性
$$|H(jw)|^2=\frac{1}{1+\varepsilon ^2C_N^2(\frac{w}{w_c})}$$
N：濾波器的階數
$\varepsilon$：通帶波紋
$\omega _c$：通帶截頻
圖 CB I型低通濾波器的幅度響應
特點：通帶是等波動的，阻帶是單調的
scipy.signal.cheby1(N, rp, Wn, btype='low', analog=False, output='ba')
Chebyshev I型數字和模擬濾波器，設計N階數字或模擬Chebyshev I型濾波器並返回濾波器系數。
參數：
- N：濾波器的階數
- rp：通帶中允許的最大紋波低於單位增益，以分貝為單位，正數
- Wn：對於數字濾波器：Wn應該歸一化為(0,1)，Wn=截止頻率/信號頻率，（信號頻率=采樣率的一半，奈奎斯特采樣定理）對於模擬濾波器：Wn是角頻率，弧度/樣本，rad/s
- btype：濾波器的類型{'lowpass'，'highpass'，'bandpass'，'bandstop'}
- analog：如果為True，則返回模擬濾波器，否則返回數字濾波器。
- output：默認“ba”,輸出分子和分母
返回
- b,a：濾波器系數， a為分母，b為分子。
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()
切比雪夫II形狀濾波器
CB II型低通濾波器的頻域特性
$$|H(jw)|^2=1-\frac{1}{1+\varepsilon ^2C_N^2(\frac{w}{w_c})}$$
N：濾波器的階數
$\varepsilon$：阻帶波紋
$\omega _c$：阻帶截頻
圖 CB II型低通濾波器的幅度響應
特點：通帶是單調的，阻帶是等波動的
scipy.signal.cheby2(N, rs, Wn, btype='low', analog=False, output='ba')
Chebyshev II型數字和模擬濾波器，設計N階數字或模擬Chebyshev II型濾波器並返回濾波器系數。
參數：
- N：濾波器的階數
- rs：阻帶所需最小衰減，以分貝為單位，正數
- Wn：對於數字濾波器：Wn應該歸一化為(0,1)，Wn=截止頻率/信號頻率，（信號頻率=采樣率的一半，奈奎斯特采樣定理）對於模擬濾波器：Wn是角頻率，弧度/樣本，rad/s
- btype：濾波器的類型{'lowpass'，'highpass'，'bandpass'，'bandstop'}
- analog：如果為True，則返回模擬濾波器，否則返回數字濾波器。
- output：默認“ba”,輸出分子和分母
返回
- b,a：濾波器系數， a為分母，b為分子。
from 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()
橢圓低通濾波器
橢圓模擬低通濾波器的頻域特性
圖 橢圓低通濾波器的幅度相應
特點：通帶和阻帶都等波動
scipy.signal.ellip(N, rp, rs, Wn, btype='low', analog=False, output='ba')
橢圓數字和模擬濾波器，設計N階數字或模擬橢圓濾波器並返回濾波器系數。
參數：
- N：濾波器的階數
- rp：通帶中允許的最大紋波低於單位增益，以分貝為單位，正數
- rs：阻帶所需最小衰減，以分貝為單位，正數
- Wn：對於數字濾波器：Wn應該歸一化為(0,1)，Wn=截止頻率/信號頻率，（信號頻率=采樣率的一半，奈奎斯特采樣定理）對於模擬濾波器：Wn是角頻率，弧度/樣本，rad/s
- btype：濾波器的類型{'lowpass'，'highpass'，'bandpass'，'bandstop'}
- analog：如果為True，則返回模擬濾波器，否則返回數字濾波器。
- output：默認“ba”,輸出分子和分母
返回
- b,a：濾波器系數， a為分母，b為分子。
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()
下采樣方法
插值方法進行下采樣
Volodymyr Kuleshov的論文中使用抗混疊濾波器對語音信號進行下采樣，再通過三次樣條插值把下采樣信號上采樣到相同的長度。
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): """ 上采樣，每隔一步去掉語音波形的r個點，然后用三次樣條插值的方法把去掉的點補回來，有機會可以畫圖看看 :param x_lr: 音頻數據 :param r: 樣條插值前個數 :return: 樣條插值后的音頻信號 """ x_lr = x_lr.flatten() # 把x_lr數組折疊成一維的數組 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) # 樣條曲線插值系數 x_sp = interpolate.splev(i_hr, f) # 給定樣條表示的節點和系數，返回在節點處的樣條值 return x_sp yt, wav_fs = librosa.load("./48k/p225_001.wav", sr=16000, mono=True) x_lr = decimate(yt, 2) # 應用抗混疊濾波器后對信號進行下采樣，獲得低分辨率音頻，下采樣因子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) # 上采樣 plt.subplot(2, 1, 2) plt.specgram(x_lr, Fs=16000, scale_by_freq=True, sides='default') plt.show()
重采樣(signal.resample)——下采樣
利用重新采樣的方法對語音進行下采樣
scipy.signal.resample(x, num, t=None, axis=0, window=None)
沿給定軸使用傅立葉方法重新采樣x到num個樣本。因為使用傅立葉方法，所以假設信號是周期性的。
參數：
- x：要重采樣的數組
- num：重采樣信號的樣本數
返回：
- resample_x：重新采樣返回的數組
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重采樣(下采樣)
凌振華老師的下采樣方法和上面的一樣
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) # 下采樣率 16000-->8000 audio8k = librosa.core.resample(audio8k, wav_fs/2, wav_fs) # 上采樣率 8000-->16000，並不恢復高頻部分 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()
librosa.load下采樣
用librosa.load想下采樣，再不恢復頻率的情況下上采樣。
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) # 把下采樣的寫好 y_8k, fs_8k = librosa.load("./8k_sample.wav", sr=16000, mono=True) # 失去的就補不回來了 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()
