用matlab编写的四种fir滤波器

fs=1000; % 采样频率 wp = 10*2*pi/1000; ws = 12*2*pi/1000; %通带截止频率10HZ，阻带起始频率12HZ，转换公式为wp=f1*2*pi/fs,ws=f2*2*pi/fs(记住即可） deltaw= ws - wp; % 过渡带宽Δω的计算 N = ceil(6.6*pi/ deltaw) ; % 按海明窗计算所需的滤波器阶数N,该式子可相关数字信号处理的课本 wdham = (hamming(N)); % 海明窗计算 Wn=(wp+ws)/2; % 计算截止频率 b=ideal_lp(Wn,N); %滤波器阶数为N，所以其冲击响应系数也为N个，其冲击响应系数序号从0开始直到N-1， %共N个，故fir1函数的第一个参数为N-1，fir1的返回值b为滤波器冲击响应系数 [H,w]=freqz(b,1); %freqz相当于离散傅里叶变换，函数freqz的调用格式为 %[H,w]=freqz(B,A，N），其中B和A分别为离散系统的系统函数分子、分母多项式的系数向量 %返回量H则包含了离散系统频响N个频率等分点的值（其中N为正整数）， %w则包含了范围内N个频率等分点。调用默认的N时，其值是512，由于 %fir滤波器分母为1，故第二个参数值为1，第三参数为默认值 db=20*log10(abs(H)); % 画频响曲线 figure(1); plot(w*1000/(2*pi),db);title('幅度响应（单位： dB）');grid %w*1000/(2*pi)再次将横轴单位变为HZ，而非弧度，与第二行对应来看 axis([0 50 -100 10]); xlabel('频率（单位：Hz）'); ylabel('分贝') set(gca,'XTickMode','manual','XTick',[0,10,12,50]) set(gca,'YTickMode','manual','YTick',[-50,0]) T=1024;%采样点数 n=0:T-1; t=n/fs; y=20*sin(2*pi*5*t)+30*cos(2*pi*6*t)+250*sin(2*pi*10*t)+100*sin(2*pi*25*t)+150*sin(2*pi*30*t); figure(2); plot(y); x=fft(y,T); %对输入信号进行快速傅里叶变换，求其频谱图 mag=abs(x); %求取Fourier变换的振幅 F=n*fs/T; figure(3); plot(F,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; Y=fftfilt(b,y);%调用fir滤波器 figure(4); plot(Y); X=fft(Y,T); mag=abs(X); %求取Fourier变换的振幅 F=n*fs/T; figure(5); plot(F,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; %*************************************** wp1=20*2*pi/fs; ws1=22*2*pi/fs; wp2=210*2*pi/fs; ws2=220*2*pi/fs; tr_width=ws1-wp1; %过度带宽 M=ceil(6.6*pi/tr_width)+1; %滤波器阶数 wc1=(ws1+wp1)/2; wc2=(ws2+wp2)/2; %通带截止频率 hd=ideal_lp(wc2,M)-ideal_lp(wc1,M); %通带滤波器的设计可以看成是两个低通滤波器想减 w_ham=(hamming(M))'; h1=hd .* w_ham; [H1,w1]=freqz(h1,1); db=20*log10(abs(H1)); n=0:M-1; figure(6); subplot(2,2,1);stem(n,hd);title('理想脉冲响应'); axis([0 1700 -0.1 0.1]);xlabel('n');ylabel('hd(n)'); subplot(2,2,2);stem(n,w_ham);title('海明窗'); axis([0 1700 0 1]);xlabel('n');ylabel('w_ham(n)'); subplot(2,2,3);stem(n,h1);title('实际脉冲响应'); axis([0 1700 -0.1 0.1]);xlabel('n');ylabel('h(n)'); figure(7); plot(w*fs/(2*pi),db);title('幅度响应（单位： dB）');grid %w*fs/(2*pi)再次将横轴单位变为HZ，而非弧度，与第二行对应来看 axis([-10 600 -200 0]); xlabel('频率（单位：Hz）'); ylabel('分贝') set(gca,'XTickMode','manual','XTick',[0,3,500,600]) set(gca,'YTickMode','manual','YTick',[-10,0]) n=0:T-1; t=n/fs; figure(8); plot(y); x1=fft(y,T); %对输入信号进行快速傅里叶变换，求其频谱图 mag=abs(x1); %求取Fourier变换的振幅 F1=n*fs/T; figure(9); plot(F1,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; Y1=fftfilt(h1,y);%调用fir滤波器 figure(10); plot(Y1); X1=fft(Y1,T); mag=abs(X1); %求取Fourier变换的振幅 F1=n*fs/T; figure(11); plot(F1,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on;