shiyan59.m

本人提供的是利用MATLAB语言设计FIR滤波器的源代码

function yn=shiyan59
k=10;              %记录长度为信号周期的k倍
[xn,freq,N,fs,Tt]=shiyan50(k);
T=1/fs;
fp=[9800,10200];
fst=[9200,10800];
rs=60;
wp=(2*pi/fs).*fp;
ws=(2*pi/fs).*fst;
b=shiyan52(wp,ws,rs);  %生成带通数字滤波器 
yn=filter(b,1,xn);        %对信号xn进行滤波
kk=5; %kk<k
ynn=yn(kk*N/k+1:N);   %去掉yn的前kk个周期
knn=abs(fft(ynn));       %计算滤波后信号ynn的频谱
M=length(ynn);
t=kk*Tt:T:(k*Tt-T);
m=[0:M-1];
f=m.*(fs/M);
figure(9)
subplot(2,1,1)
plot(t,ynn);
Xlabel('时间/S');
legend('输出信号的波形');
subplot(2,1,2)
stem(f(1:M/2),knn(1:M/2));
Xlabel('频率/Hz');
legend('输出信号的频谱');

function [xn,freq,N,fs,Tt]=shiyan50(k)  %生成信号序列
fs=40000;            %设定抽样频率fs/Hz
Tt=0.0001;           %信号周期Tt/S
T=1/fs;
T0=k*Tt;            %设定记录长度T0/S
t=0:T:(T0-T);
xn=2.*(1+0.5.*cos(2000*pi*t)).*cos(20000*pi*t); %信号序列xn
xk=fft(xn);           %计算信号序列xn的频谱xk
N=length(xn);         %计算信号序列的长度
i=1; 
m=[0:N-1];
f=m.*(fs/N);          %将数字频率转换成模拟线性频率/Hz
freq=0;              % freq记录信号xn的频率分量
for m=1:1:(N+1)/2  
   if (abs(xk(m))>0.0001)
       freq(i)=m-1;
       i=i+1;
    end
end
figure(1)
subplot(2,1,1)
plot(t,xn);Xlabel('时间');legend('信号波形');
subplot(2,1,2)
stem(f(1:N/2),xk(1:N/2),'r');
Xlabel('频率/Hz');legend('信号频谱');
grid

function [windowxing,jieshu]=shiyan51(rs,wguo)
%根据rs、过渡带wguo的大小选择窗形windowxing，确定阶数jieshu
if rs>0
N(1)=ceil(4*pi/wguo);as(1)=21;str(1).window=boxcar(N(1)+1);   N(2)=ceil(8*pi/wguo);as(2)=25;str(2).window=triang(N(2)+1);    N(3)=ceil(8*pi/wguo);as(3)=44;str(3).window=hanning(N(3)+1);   N(4)=ceil(8*pi/wguo);as(4)=53;str(4).window=hamming(N(4)+1);  N(5)=ceil(12*pi/wguo);as(5)=74;str(5).window=blackman(N(5)+1);    N(6)=ceil(10*pi/wguo);as(6)=80;str(6).window=kaiser(N(6)+1,7.865);
    i=0;
    while i<6
         i=i+1;
         if as(i)<rs
             N(i)=1000000;
         end
     end
    [Nmin,m]=min(N)
    windowxing=str(m).window;
    jieshu=Nmin;
else 
    error('rs > 0')
end

function b=shiyan52(wp,ws,rs)
%生成滤波器
    wguo=min(abs(wp-ws));
    wn=(wp+ws)./2;
    [windowxing,N]=shiyan51(rs,wguo);
    b=fir1(N,wn/pi,windowxing);
figure(2)
[H,W]=freqz(b,1);
plot(W/pi,20*log10(abs(H)));
grid