Matlab實現IIR數字濾波器設計

低通濾波器的設計：

1、基本指標

　　低通濾波器顧名思義就是低頻率成分通過，高頻率成分截止，那么在設計一個低通濾波器時首先要明白想要截止多大的頻率。比如想截掉4Hz以上的信號，理想狀態下就是將4Hz以上的信號成分全部截止，4Hz以下的信號全部保留。然而事實上幾乎不存在這樣的濾波器，通常情況下總是在通過頻率和截止頻率之間存在一個過渡帶。通過頻率這部分稱為通帶，允許通過的最大頻率為通帶截止頻率ω_p，截至頻率這部分稱為阻帶，阻帶最小截至頻率為ω_s，通帶和阻帶之間的部分為過渡帶，也即ω_p~ω_p。通帶之間的波動稱為通帶波動δ_p，阻帶之間的波動稱為阻帶波動δ_s如圖所示：

2、Matlab函數介紹

聲明：上面一節的基本指標均為模擬濾波器狀態下的指標，Matlab不僅提供了模擬濾波器設計函數，同樣提供了數字濾波器的設計函數，只是參數之間需要一些轉換。

數字濾波器設計中，Wp和Ws為歸一化角頻率，如果工程要求設計低通數字濾波器滿足通帶截止頻率為fp = 40，阻帶截止頻率fs = 50，則fp，fs和Wp，Ws之間的轉換關系為：

Wp = 2*pi*fp/Fs;
Ws = 2*pi*fs/Fs;

截止頻率的定義是當輸出幅值響應下降到輸入幅值的-3dB (  )，也就是0.707(也就是  )時對應的頻率，也稱3db截頻

巴特沃斯數字濾波器函數：
[n,Wc] = buttord(Wp/pi,Ws/pi,Ap,As);
[b,a] = butter(n,Wc,'low');%低通濾波器

Ap和As的計算用上圖中的公式計算獲得。

function FilteredSignal = filtered(signal,Fs,fp,fs,Ap,As,FilterType)
% fp=40; fs=50; Ap=1; As=20;
%=============================================
% Fs:信號的頻率
% fp：濾波器的通帶頻率
% fs:濾波器的截止頻率
% Ap：通帶最大衰減
% As：阻帶最小衰減
% FilterType:濾波器類型，'low',低通，'high',高通
% FilteredSignal：濾波后的信號
%==============================================
Wp = 2*pi*fp/Fs;
Ws = 2*pi*fs/Fs;
[n,Wc] = buttord(Wp/pi,Ws/pi,Ap,As);
[b,a] = butter(n,Wc,FilterType);
omega=[Wp Ws]; 
h = freqz(b,a,omega); %Compute Ap and As of AF
fprintf('Ap= %.4f\n',-20*log10(abs(h(1))));
fprintf('As= %.4f\n',-20*log10(abs(h(2))));
FilteredSignal = filter(b,a,signal);

　當 fp=40; fs=50; Ap=1; As=20;時，通過上面的代碼設計的濾波器Ap= 0.9998，As= 20.0000，滿足要求，設計的濾波器如下圖：

　

切比雪夫I型數字濾波器函數：

[N,wc] = cheb1ord(Wp/pi,Ws/pi,Ap,As);
[b,a] = cheby1(N,Ap,wc,'low');

%% 切比雪夫1型低通濾波器
fp=40; fs=50; Ap=1; As=20;
Wp = 2*pi*fp/Fs;
Ws = 2*pi*fs/Fs;
[N,wc] = cheb1ord(Wp/pi,Ws/pi,Ap,As);
[b,a] = cheby1(N,Ap,wc,'low');
omega = [Wp,Ws];
h = freqz(b,a,omega);
fprintf('N = %.2f\n',N);
fprintf('Ap=%.4f\n',-20*log10(abs(h(1))));
fprintf('As=%.4f\n',-20*log10(abs(h(2))));
figure(7)
freqz(b,a,512,200);title('CBI Lowpass Filter')

　

N = 5.00；Ap=1.0000<=1；As=20.0093>20滿足定義的Ap和As，通帶剛好，阻帶有余量。

%% 橢圓低通濾波器
[N,wc] = ellipord(Wp/pi,Ws/pi,Ap,As);
[b,a] = ellip(N,Ap,As,wc,'low');
omega = [Wp,Ws];
h = freqz(b,a,omega);
fprintf('N = %.2f\n',N);
fprintf('Ap=%.4f\n',-20*log10(abs(h(1))));
fprintf('As=%.4f\n',-20*log10(abs(h(2))));
figure
freqz(b,a,512,200);title('elip Lowpass Filter')