模擬和數字低通濾波器的MATLAB實現
轉自：http://blog.sina.com.cn/s/blog_79ecf6980100vcrf.html
低通濾波器引數：Fs=8000，fp=2500，fs=3500，Rp=1dB，As=30dB，其他濾波器可以通過與低通之間的對映關係實現。
%%模擬濾波器
%巴特沃斯——濾波器設計
wp=2*pi*2500;ws=2*pi*3500;Rp=1;As=30;
[N,wc]=buttord(wp,ws,Rp,As,'s')%計算率波器的階數和3dB截止頻率
[B,A]=butter(N,wc,'s');%計算濾波器系統函式分子分母多項式
fk=0:800/512:8000;wk=2*pi*fk;
Hk=freqs(B,A,wk);
figure
plot(fk/1000,20*log10(abs(Hk)));
grid on,xlabel('頻率（kHz）'),ylabel('幅度（dB）')
title('巴特沃斯模擬濾波器')
axis([0,4,-35,5])
%%
%切比雪夫I——濾波器設計
wp=2*pi*2500;ws=2*pi*3500;Rp=1;As=30;
[N1,wp1]=cheb1ord(wp,ws,Rp,As,'s')%計算切比雪夫濾波器的階數和通帶邊界頻率
[B1,A1]=cheby1(N1,Rp,wp1,'s'
fk=0:800/512:8000;wk=2*pi*fk;
Hk=freqs(B1,A1,wk);figure,
plot(fk/1000,20*log10(abs(Hk)));
grid on,xlabel('頻率（kHz）'),ylabel('幅度（dB）')
title('切比雪夫I模擬濾波器')
axis([0,4,-35,5])
%%
%切比雪夫II——濾波器設計
wp=2*pi*2500;ws=2*pi*3500;Rp=1;As=30;
[N2,wso]=cheb2ord(wp,ws,Rp,As,
[B2,A2]=cheby2(N1,Rp,wso,'s');%計算濾波器系統函式分子分母多項式
fk=0:800/512:8000;wk=2*pi*fk;
Hk=freqs(B1,A1,wk);figure,
plot(fk/1000,20*log10(abs(Hk)));
grid on,xlabel('頻率（kHz）'),ylabel('幅度（dB）')
title('切比雪夫II模擬濾波器')
axis([0,4,-35,5])
%%
%橢圓——濾波器設計
wp=2*pi*2500;ws=2*pi*3500;Rp=1;As=30;
[N,wpo]=ellipord(wp,ws,Rp,As,'s')%計算濾波器的階數和通帶邊界頻率
[B,A]=ellip(N,Rp,As,wpo,'s');%計算濾波器系統函式分子分母多項式
fk=0:800/512:8000;wk=2*pi*fk;
Hk=freqs(B1,A1,wk);figure,
plot(fk/1000,20*log10(abs(Hk)));
grid on,xlabel('頻率（kHz）'),ylabel('幅度（dB）')
axis([0,4,-35,5]),title('橢圓模擬濾波器')
%%
%數字濾波器
%脈衝響應法濾波器設計
fp=2500;fs=3500;Fs=8000;
wp=2*fp/Fs;ws=2*fs/Fs;%求歸一化數字通帶截止頻率,求歸一化數字阻帶起始頻率
deltaw=ws-wp;%求過渡頻寬
N0=ceil(6.6/deltaw);%求視窗長度
N=N0+mod(N0+1,2); %確保視窗長度 N為奇數
n=N-1;%求出濾波器的階數 n
wn=(ws+wp)/2; %求濾波器的截止頻率
b=fir1(n,wn)%利用 fir1 函式求出濾波器的係數
[Hk,w] = freqz(b,1); % 計算頻率響應
mag = abs(Hk); % 求幅頻特性
db = 20*log10(mag/max(mag)); % 化為分貝值
dw =pi/512; %關於pi歸一化
Rp = -(min(db(1:wp*pi/dw+1))) % 檢驗通帶波動
As = -(max(db(ws*pi/dw+1:512))) % 檢驗最小阻帶衰減
figure,plot(0:pi/511:pi,db),grid on
axis([0,4.0,-80,5]),title('數字濾波器——脈衝響應法')
%%
%fir1窗函式法
fp=2500;fs=3500;Fs=8000;rs=30;
wp=2*fp*pi/Fs;ws=2*fs*pi/Fs;%求歸一化數字通帶截止頻率,求歸一化數字阻帶起始頻率
Bt=ws-wp;%求過渡頻寬
alpha=0.5842*(rs-21)^0.4+0.07886*(rs-21);%計算kaiser窗的控制引數
M=ceil((rs-8)/2.285/Bt);%求出濾波器的階數
wc=(ws+wp)/2/pi; %求濾波器的截止頻率並關於pi歸一化
hk=fir1(M,wc,kaiser(M+1,alpha))%利用 fir1 函式求出濾波器的係數
[Hk,w] = freqz(hk,1); % 計算頻率響應
mag = abs(Hk); % 求幅頻特性
db = 20*log10(mag/max(mag)); % 化為分貝值
db1=db';
figure,plot(0:pi/511:pi,db1),grid on
axis([0,4.0,-80,5]),title('數字濾波器——fir1窗函式法')
%%
%頻率取樣法
fp=2500;fs=3500;Fs=8000;rs=30;
wp=2*fp*pi/Fs;ws=2*fs*pi/Fs;%求歸一化數字通帶截止頻率,求歸一化數字阻帶起始頻率
Bt=ws-wp;%求過渡頻寬
m=1;alpha=0.5842*(rs-21)^0.4+0.07886*(rs-21);%計算kaiser窗的控制引數
N=ceil(m+1)*2*pi/Bt;%求出濾波器的階數
N=N+mod(N+1,2);
Np=fix(wp/(2*pi/N));
Ns=N-2*Np-1;
Hk=[ones(1,Np+1),zeros(1,Ns),ones(1,Np)];
wc=(ws+wp)/2/pi; %求濾波器的截止頻率並關於pi歸一化
hk=fir1(M,wc,kaiser(M+1,alpha))%利用 fir1 函式求出濾波器的係數
[Hk,w] = freqz(hk,1); % 計算頻率響應
mag = abs(Hk); % 求幅頻特性
db = 20*log10(mag/max(mag)); % 化為分貝值
db1=db';
figure,plot(0:pi/511:pi,db1),grid on
axis([0,4.0,-80,5]),title('數字濾波器——頻率取樣法')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%利用等波紋最佳逼近法設計FIR數字濾波器
Fs=8000;f=[2500,3500];m=[1,0];
rp=1;rs=30;
delta1=(10^(rp/20)-1)/(10^(rp/20)+1);delta2=10^(-rs/20);
rip=[delta1,delta2];
[M,fo,mo,w]=remezord(f,m,rip,Fs);%邊界頻率為模擬頻率時必須加入取樣頻率
M=M+1;%估算的M直達不到要求，家1後滿足要求
hn=remez(M,fo,mo,w);
[Hk,w] = freqz(hn,1); % 計算頻率響應
mag = abs(Hk); % 求幅頻特性
db = 20*log10(mag/max(mag)); % 化為分貝值
db1=db';
figure,plot(0:pi/511:pi,db1),grid on
axis([0,4.0,-80,5]),title('數字濾波器——等波紋最佳逼近法')