实验8IIR数字滤波器设计综合实验(预习报告)

实验8IIR数字滤波器设计综合实验1．实验程序及运行结果信号产生函数：mstg.mfunctionst=mstg%产生信号序列向量st，并显示st的时域波形和频谱%st=mstg返回三路调幅信号相加形成的混合信号，长度N=1600N=1600;%N为信号st的长度。Fs=10000;T=1/Fs;Tp=N*T;%采样频率Fs=10kHz，Tp为采样时间t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=Fs/10;%第1路调幅信号的载波频率fc1=1000Hz,fm1=fc1/10;%第1路调幅信号的调制信号频率fm1=100Hzfc2=Fs/20;%第2路调幅信号的载波频率fc2=500Hzfm2=fc2/10;%第2路调幅信号的调制信号频率fm2=50Hzfc3=Fs/40;%第3路调幅信号的载波频率fc3=250Hz,fm3=fc3/10;%第3路调幅信号的调制信号频率fm3=25Hzxt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t);%产生第1路调幅信号xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);%产生第2路调幅信号xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);%产生第3路调幅信号st=xt1+xt2+xt3;%三路调幅信号相加fxt=fft(st,N);%计算信号st的频谱%====以下为绘图部分，绘制st的时域波形和幅频特性曲线====================subplot(2,1,1)plot(t,st);grid;xlabel('t/s');ylabel('s(t)');axis([0,Tp/8,min(st),max(st)]);title('(a)s(t)的波形')subplot(2,1,2)stem(f,abs(fxt)/max(abs(fxt)),'.');grid;title('(b)s(t)的频谱')axis([0,Fs/5,0,1.2]);xlabel('f/Hz');ylabel('幅度')函数说明：（1）该函数生成的三路时域信号分别为1()cos(21000)cos(2100)yttt2()cos(2500)cos(250)yttt3()cos(2250)cos(225)yttt（2）合成信号为123()()()()stytytyt（3）采样频率为10000sFHz，采样点数为1600N点，采样时间pT0.16NTs。（4）抑制载波单频调幅信号的数学表示式为0001()cos(2)cos(2)cos(2())cos(2())2cccstftftfftfft源程序：shzxhchlshiyan4_1%shzxhchlshiyan4_1%IIR数字滤波器设计及软件实现clearall,closeall,clc,clf;Fs=10000;T=1/Fs;%采样频率%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号stst=mstg;%低通滤波器设计与实现===========================================fp=280;fs=450;wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;%椭圆低通滤波器技术指标（低通滤波器的通、阻带边界频）[N,wp]=ellipord(wp,ws,rp,rs);%调用ellipord计算椭圆低通滤波器阶数N和通带截止频率wp[B,A]=ellip(N,rp,rs,wp);%调用ellip计算椭圆低通滤波器系统函数系数向量B和Ay1t=filter(B,A,st);%求低通滤波器的零状态响应[H,w]=freqz(B,A);%求低通滤波器的频响特性%绘制椭圆低通滤波器衰减特性曲线figure(2);mag=abs(H);db=20*log10((mag+eps)/max(mag));%eps浮点数的相对误差subplot(2,1,1);plot(w/pi,db);%绘制损耗函数曲线axis([0,1,-100,5]);gridon;title('低通滤波损耗函数曲线');xlabel('\omega/\pi');ylabel('幅度(dB)');set(gca,'Xtick',[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]);set(gca,'Ytick',[-100,-80,-60,-40,-20,0]);%绘制225-275Hz范围内时域信号波形N=1600;t=0:T:(N-1)*T;subplot(2,1,2);plot(t,y1t);title('低通滤波分离225-275Hz范围内时域信号波形');xlabel('t/s');ylabel('y_1(t)');axis([0,0.16,-1.1,1.1]);%带通滤波器设计与实现===========================================fpl=440;fpu=560;fsl=275;fsu=900;wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;[N,wp]=ellipord(wp,ws,rp,rs);[B,A]=ellip(N,rp,rs,wp);y2t=filter(B,A,st);%求帯通滤波器的零状态响应[H,w]=freqz(B,A);%求帯通滤波器的频响特性%绘制椭圆低通滤波器衰减特性曲线figure(3);mag=abs(H);db=20*log10((mag+eps)/max(mag));subplot(2,1,1);plot(w/pi,db);axis([0,1,-100,5]);gridon;title('带通滤波损耗函数曲线');xlabel('\omega/\pi');ylabel('幅度(dB)');set(gca,'Xtick',[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]);set(gca,'Ytick',[-100,-80,-60,-40,-20,0]);%绘制450-550Hz范围内时域信号波形N=1600;t=0:T:(N-1)*T;subplot(2,1,2);plot(t,y2t);title('带通滤波分离450-550Hz范围内时域信号波形');xlabel('t/s');ylabel('y_2(t)');axis([0,0.16,-1.1,1.1]);%高通滤波器设计与实现================================================fp=890;fs=600;wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;[N,wp]=ellipord(wp,ws,rp,rs);[B,A]=ellip(N,rp,rs,wp,'high');y3t=filter(B,A,st);[H,w]=freqz(B,A);%绘制椭圆高通滤波器衰减特性曲线figure(4);mag=abs(H);db=20*log10((mag+eps)/max(mag));subplot(2,1,1);plot(w/pi,db);axis([0,1,-100,5]);gridon;title('高通滤波损耗函数曲线');xlabel('\omega/\pi');ylabel('幅度(dB)');set(gca,'Xtick',[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]);set(gca,'Ytick',[-100,-80,-60,-40,-20,0]);%绘制900Hz以上时域信号波形N=1600;t=0:T:(N-1)*T;subplot(2,1,2);plot(t,y3t);xlabel('t/s');ylabel('y_3(t)');axis([0,0.16,-1.1,1.1]);title('高通滤波分离900-1100Hz范围内时域信号波形');程序运行结果：00.0020.0040.0060.0080.010.0120.0140.0160.0180.02-10123t/ss(t)(a)s(t)的波形020040060080010001200140016001800200000.51(b)s(t)的频谱f/Hz幅度图4-1三路调幅信号st的时域波形和幅频特性曲线00.10.20.30.40.50.60.70.80.91-100-80-60-40-200低通滤波损耗函数曲线/幅度(dB)00.020.040.060.080.10.120.140.16-1-0.500.51低通滤波分离225-275Hz范围内时域信号波形t/sy1(t)图4-2椭圆低通滤波器的幅频特性曲线和输出时域波形00.10.20.30.40.50.60.70.80.91-100-80-60-40-200带通滤波损耗函数曲线/幅度(dB)00.020.040.060.080.10.120.140.16-1-0.500.51带通滤波分离450-550Hz范围内时域信号波形t/sy2(t)图4-3椭圆带通滤波器的幅频特性曲线和输出时域波形00.10.20.30.40.50.60.70.80.91-100-80-60-40-200高通滤波损耗函数曲线/幅度(dB)00.020.040.060.080.10.120.140.16-1-0.500.51高通滤波分离900-1100Hz范围内时域信号波形t/sy3(t)图4-4椭圆高通滤波器的幅频特性曲线和输出时域波形程序运行结果分析：由图4-2~图4-4可见：三个分离滤波器指标参数选取正确，算耗函数曲线达到所给指标。分离出的三路信号1()yn，2()yn和3()yn的波形是抑制载波的单频调幅波。2．思考题简答（1）请阅读信号产生函数mstg，确定三路调幅信号的载波频率和调制信号频率。答：第一路调幅信号的载波频率为11000Hzcf，调制信号频率为1100Hzmf；第二路调幅信号的载波频率为2500Hzcf，调制信号频率为250Hzmf；第三路调幅信号的载波频率为3250Hzcf，调制信号频率为325Hzmf。（2）信号产生函数mstg中采样点数N=1600，对st进行N点FFT可以得到6根理想谱线。如果取N=1800，可否得到6根理想谱线？为什么？N=2000呢？请改变函数mstg中采样点数N的值，观察频谱图验证您的判断是否正确。答：因为信号st是周期序列，谱分析时要求观察时间为整数倍周期。所以，本题的一般解答方法是，先确定信号st的周期，在判断所给采样点数N对应的观察时间pTNT是否为st的整数个周期。但信号产生函数mstg产生的信号st共有6个频率成分，求其周期比较麻烦，故采用下面的方法解答。分析发现，st的每个频率成分都是25Hz的整数倍。采样频率Fs=10kHz=25×400Hz，即在25Hz的正弦波的1个周期中采样400点。所以，当N为400的整数倍时一定为st的整数个周期。因此，采样点数N=800和N=2000时，对st进行N点FFT可以得到6根理想谱线。如果取N=1000，不是400的整数倍，不能得到6根理想谱线。（3）修改信号产生函数mstg，给每路调幅信号加入载波成分，产生调幅（AM）信号，重复本实验，观察AM信号与抑制载波调幅信号的时域波形及其频谱的差别。答：AM信号表示式dm0dm()A+Acos(2)cos(2)AAcstftft按上式修改程序，观察运行结果。