低通濾波器實現過程解析
首先來看一個低通濾波的例子：給定一組波形如下：
如果想要獲取低通0-0.5Hz的波形：需要呼叫低通濾波器可得到：
很簡單的一個例子：程式見：點選開啟連結
低通濾波器實現過程：
主函式：
|
y=lowp(x,f1,f3,rp,rs,Fs) |
X:輸入訊號；Fs: 取樣頻率 f1:通帶截止頻率，rp:通帶衰減 f3:阻帶截止頻率，rs:阻帶衰減 |
Step1：模擬頻率/數字頻率轉換
|
wp=2*pi*f1/Fs; ws=2*pi*f3/Fs; |
Step2：求Chebyshev Type I filter order
|
[n,wn]=cheb1ord(wp/pi,ws/pi,rp,rs); N:濾波器階數 Wn：濾波器截止頻率 |
wp/pi:歸一化通帶截止頻率 ws/pi：歸一化阻帶截止頻率 rp：通帶衰減 rs:阻帶衰減 |
Step3：求濾波器分子、分母系數
|
[bz1,az1]=cheby1(n,rp,wp/pi); bz1：分子係數矩陣； az1：分母系數矩陣 |
N：階數 Rp：通帶衰減 wp/pi：歸一化通帶衰減頻率 |
Step3.1：get analog, pre-warped frequencies
|
if ~analog, fs = 2; u = 2*fs*tan(pi*Wn/fs); else u = Wn; end |
Step3.2:convert to low-pass prototype estimate
|
if btype == 1 % lowpass Wn = u; elseif btype == 2 % bandpass Bw = u(2) - u(1); Wn = sqrt(u(1)*u(2)); % center frequency elseif btype == 3 % highpass Wn = u; elseif btype == 4 % bandstop Bw = u(2) - u(1); Wn = sqrt(u(1)*u(2)); % center frequency end |
Step3.3:Get N-th order Chebyshev type-I lowpass analog prototype
|
[z,p,k] = cheb1ap(n, r); Z:濾波器零點 P:濾波器極點 K:係數 |
N:濾波器階數 r:通帶衰減 |
包含：asinh，exp，flipud，complex，prod |
|
|
[a,b,c,d] = zp2ss(z,p,k); 狀態方程 x = Ax + Bu y = Cx + Du A,b,c,d:狀態方程的係數 |
z,p,k：同上 |
包含：parse_input， zp2tf，tf2ss，isfinite cplxpair，poly. 待處理：zp2tf，tf2ss 異常處理：Try...catch |
零極點轉化為狀態方程 |
Step3.4:Transform to lowpass, bandpass, highpass, or bandstop of desired Wn
|
[a,b,c,d] = lp2lp(a,b,c,d,Wn); |
a,b,c,d：同上 Wn:截止頻率 |
Step3.5 ：Use Bilinear transformation to find discrete equivalent:
|
[a,b,c,d] = bilinear(a,b,c,d,fs); |
雙線性變換 |
|
den = poly(a); |
分子 |
|
num = cheb1num(btype,n,Wn,Bw,analog,den,r); |
分母 |
Step4：差分方程求濾波後的序列
|
y=filter(bz1,az1,x); |
bz1：分子係數矩陣； az1：分母系數矩陣 X：輸入序列 |
-------------------------------------------------------------------------------------
matlab實現是相當簡單的，問題難得是如果要用C語言實現肯定是要費一番功夫的。
under-edit