fir.m

Matlab设计FIR,IIR

close all,clear all,clc;

N=(40-7.95)/2.286/(0.15*pi);
N=round(N)+1                          %N取奇数
omega=[0.15*pi,0.3*pi,0.5*pi,0.65*pi];
wc1=(omega(1)+omega(2))/2;            %理想带通的截止频率
wc2=(omega(3)+omega(4))/2;
a=(N-1)/2;                            %满足线性相位的条件

for n=0:100                           %此处我只画出理想带通滤波器101个点的单位样值
    if n==a                           %响应，当然还可以取其它值
        hd(n+1)=(wc2-wc1)/pi;
    else
    hd(n+1)=(sin((n-a)*wc2)-sin((n-a)*wc1))/pi/(n-a);  %由频率特性求单位样值相应的公式
    end
end
n=0:100;
stem(n,hd);

ns=40;
b=0.5842*(ns-21)^0.4+0.07886*(ns-21)
for m=0:N-1
    b1=b*sqrt(1-(1-2*m/(N-1))^2);
    wk(m+1)=besl(b1)/besl(b);       %调用besl函数，求得kaiser窗函数
end
m=0:N-1;
figure;
stem(m,wk);                         %画出窗函数

for n=0:N-1
    h(n+1)=hd(n+1)*wk(n+1);         %加窗后的带通滤波器的单位样值响应
end
figure;
stem(m,h);


k=2*pi/(pi/40)+1;
H=zeros(1,k);                       %对数字带通滤波器频率特性初始化
q=0;
for w=-pi:pi/40:pi
    q=q+1;
    for n2=0:N-1
        H(q)=H(q)+h(n2+1)*exp(-j*w*n2);%计算数字滤波器的幅频特性
    end
end
figure,hold on,grid on;
H=20*log10(abs(H));                    %转换为对数幅频特性
w=-pi:pi/40:pi;
w=w/pi;
plot(w,H);