发布时间 : 星期二 文章数字信号处理上机指导 - 1dl更新完毕开始阅读a620bd4b852458fb770b56fc
[Ba,Aa,Wn]=shiyan41(mp,ms,rp,rs);
[b,a]=lp2lp(Ba,Aa,Wn); % 将归一化频率的低通滤波器转换成截止频
% 率为Wn的数字低通滤波器
[bz,az]=bilinear(b,a,fs); [H,W]=freqz(bz,az); figure(2);
plot(W*fs/(2*pi),20*log10(abs(H))); hold on
plot(W0,-rr,'r') Xlabel('频率/Hz')
ylabel('幅频特性/dB') grid
function [bz,az]=shiyan43(fp,fst,rp,rs,fs) %数字高通滤波器的生成
W0=[0,fst,fp,fs/2]; rr=[rs,rs,rp,rp]; %设计指标
wp=fp*2*pi/fs; ws=fst*2*pi/fs; %模拟指标数字化 C=2*fs; %频率预畸 hp=C*tan(wp/2); hs=C*tan(ws/2);
mp=hp; ms=hp^2/hs; %模拟高通指标转化为模拟低通指标 [Ba,Aa,Wn]=shiyan41(mp,ms,rp,rs);
[b,a]=lp2hp(Ba,Aa,mp); %将归一化频率的低通滤波器转换成截止
%频率为Wn的高通滤波器
[bz,az]=bilinear(b,a,fs); [H,W]=freqz(bz,az); figure(3)
plot(W*fs/(2*pi),20*log10(abs(H))); hold on
plot(W0,-rr,'r') Xlabel('频率/Hz')
ylabel('幅频特性/dB') grid
实验五 FIR滤波器的实现与应用
一、实验目的
1．掌握用窗函数法设计FIR数字滤波器的方法。 2．应用FIR数字滤波器实现对信号的处理。
二、实验内容
设调制信号x(t)=3(1+0.5cos(2π×103t)) cos(2π×106t)，试截取它的载波信号，并绘出它的时域波形及频谱图。 三、实验原理
（一）用窗函数法设计数字滤波器
一个截止频率为ω0的理想数字低通滤波器，其表达式如下：
?1Hd(?)???0|?|??0?0?|?|??
故其冲激响应序列为
hd(n)?12?????H(?)ej?nd???0?sinc(?0?n)
- 32 -
这个滤波器是物理不可实现的，加窗函数ω(n)将其截短，保留冲激响应hd(n)的中心部分，就可获得一个线性相位的FIR滤波器h(n)。即
h(n)=hd(n)·ω(n)
最简单的窗函数是矩形窗，其ω(n)=RN(n)，问题是在通带和阻带的边缘存在波动（吉布斯效应），该现象不能通过增加滤波器长度来消除，只能采用非线性窗如汉宁窗、海明窗等。
利用窗函数法设计FIR滤波器的步骤如下： 1．确定所要设计的数字滤波器的性能指标，确定理想的频率响应函数Hd(e)并计算hd(n)； 2．由最小阻带衰减选择窗的形状，由过渡带的宽度确定N的大小（参见教材的窗函数基本参数表）；
3．求得所设计的滤波器的单位抽样响应h(n)；
jω
4．计算H(e)=DTFT[h(n)]，检验是否满足设计指标。
在Matlab环境下，函数fir1就是基于这种窗函数方法的。给定滤波器的阶数和期望的理想滤波器的描述，它返回理想滤波器的加窗逆傅里叶变换。因而用窗函数法设计FIR滤波器的步骤简化如下：
1）确定所要设计的数字滤波器的性能指标； 2）由最小阻带衰减选择窗的形状； 3）由过渡带的宽度确定N的大小； 4）调用函数fir1生成FIR滤波器。
（二）应用FIR数字滤波器实现对信号的处理
此处主要是应用FIR数字滤波器实现对调制信号的解调，即利用FIR数字滤波器选出基带信号。
四、实验步骤
1．分析实验内容，设计实验方案 例如：
2．编程实现每一个子功能并调试之
例如：FIR数字带通滤波器的设计
FIR数字带通滤波器设计指标：通带截止频率ωp1、ωp2，阻带截止频率ωs1、ωs2，阻带最小衰减rs。
首先根据rs的大小选择窗形，计算过渡带ω=min(ωp1－ωs1，ωs2－ωp2)，由此确定N的大小，程序如下：
function [windowxing,jieshu]=shiyan51(rs,wguo)
在此基础上确定FIR数字带通滤波器的系数b，编写程序如下： function b=shiyan52(wp,ws,rs)
3．选择合理的滤波器参数，编程完成信号的截取。 五、实验要求
1．针对实验内容，首先进行理论分析，给出预测结果；
2．熟悉实验教材§6-5的内容；
3．编写每一个子功能的程序，上机运行；
4．研究滤波器参数值、采样频率、记录长度的大小对实验结果的影响。 六、程序示例
function yn=shiyan59
k=10; %记录长度为信号周期的k倍
jω
- 33 -
[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 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 - 34 - 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) 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 - 35 -