fir.m

在Matlab环境下实现的FIR滤波器

%由模拟滤波器指标求理想低通数字滤波器的参数
clear
wp=0.2*pi;              %wp为通带截止频率
ws=0.4*pi;              %ws为阻带下限频率
A=-40;                  %A为最大边瓣峰值

%首先得角频率
%归一化过渡带
b=ws-wp;
k=1;

%根据最小阻带衰减来确定窗函数类型
if  A>=-13 
    k=1;disp('矩形窗');
elseif  A>=-27
    k=2;disp('三角窗');
elseif   A>=-32
    k=3;disp('汉宁窗');
elseif   A>=-43
    k=4;disp('汉明窗');
else
    k=5;disp('布莱克曼窗');
end

% 求滤波器阶数
switch k
    case 1,             %矩形窗
        N=ceil(1.8*pi/b);
        for n=1:1:N+1
            w(n)=1;
        end
    case 2,             %三角窗
        N=ceil(6.1*pi/b);
        for n=1:1:N+1
            if(n<(N+1)/2)
                w(n)=2*(n-1)/N;
            else
                w(n)=2*(N-n+1)/N;
            end
        end
    case 3,             %汉宁窗
        N=ceil(6.2*pi/b);
        for n=1:1:N+1
            w(n)=0.5-0.5*cos(2*pi*(n-1)/N);
        end            
    case 4,             %汉明窗
        N=ceil(6.6*pi/b);
        for n=1:1:N+1
            w(n)=0.54-0.46*cos(2*pi*(n-1)/N);
        end
    case 5,             %布莱克曼窗
        N=ceil(11*pi/b);
        for n=1:1:N+1
            w(n)=0.42-0.5*cos(2*pi*(n-1)/N)+0.08*cos(4*pi*(n-1)/N);
        end
end

%求截至频率            
wc=(ws+wp)/2;
n=0:N;
a=(N-1)/2;
na=n-a+eps*((n-a)==0);
%计算理想低通单位脉冲响应
hdn=sin(wc*na)/pi./na;
%加窗截取
y=hdn.*w;

%绘制单位冲激响应
figure(1);
stem(n,y,'.');
xlabel('n');
ylabel('h(n)');
grid;

%窗函数的频谱及相位
M=40;
h=freqz(w,1,M);
figure(2);
subplot(2,1,1);
n=1:N+1;
stem(n,w);grid on;
subplot(2,1,2);
f=(1:1:M)/M
plot(f,abs(h/max(abs(h))),'-');grid on;