北京工业大学信号处理工程应用训练

北京工业大学专业：通信工程学生姓名：刘莹莹指导教师：席大林完成时间：2020年6月7日通信系统工程应用训练报告目录训练十一 DFT性质研究1训练十二 DFT及抽样定理研究13训练十三 数字滤波器制作20训练十四 IIR数字滤波器设计与实现25训练十五 线性卷积计算46训练十六 FIR数字滤波器设计与实现55 训练十一 DFT性质研究验证dft函数正确性设置原始输入信号为x8=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0，将输入信号x8进行DFT正变换，dft(X,x,8,1)，输出保存在X8，如下：可以看到，输入信号x(n)已经变换到频域X(k)，且仍为8位。再对X8进行DFT反变换，dft(x,X,8,-1)，重新得到x8，观察得到的输出与原始输入数据是否相同。结果如下：可以看到，输出的x8取值仍为x8=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0，证明经过DFT正反变换后，信号能够恢复原始信号。根据帕塞瓦尔定理，应有时域、频域总能量相等：。经过计算，时域、频域能量和分别为，证明时域、频域能量和相同，符合帕塞瓦尔定理。 综上，证明DFT变换正确。A、补0效应研究原数组：x8=1,0,2,0,3,0,4,0,5,0,6,0,7,8,0示例程序中补0后数组为：x216=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0补0方式我使用的补0方式为：for(i=8;i13;i+)x2i=COMPLEX(0,0);补0后数组为：x213=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0,0,0,0,0,0,0,0,0,0,0结果分析与图在时域中，信号长度增加，由于所增加的项均为零，波形仍与未补0时相同未补零时的信号时域图补5个零后的信号时域图 补8个零后的信号时域图经过DFT变换后，X(k)长度也会随着x(n)长度的增加而增加，且增加的值非零未在末端补零时，信号频谱图在末端补5个零时，信号频谱图 在末端补8个零时，信号频谱图可以看到，经过补0，经过DFT变换的频谱与未补零时形状基本相同，只是在长度上进行扩展，且补零数量越多，扩展越长。可以理解为经过补0效应，增加了频域采样频率，但是由于信号未增加新的信息，因此不能提高物理分辨率。在能量上，补5/8个零时，信号能量时域、频域能量和如下：时域能量和、频域能量和始终相等，符合帕塞瓦尔定理，且能量与未插值时的相同。B、插值效应研究原数组：x8=1,0,2,0,3,0,4,0,5,0,6,0,7,8,0示例程序中插值后数组为：x316=1,0,8,0,2,0,7,0,3,0,6,0,4,0,5,0,5,0,4,0,6,0,3,0,7,0,2,0,8,0,1,0插值方式我使用的插值方式为：for(i=0;i16;i=i+2)x3i=COMPLEX(1+i/2,0);x3i+1=COMPLEX(i*0.5+2.5,0);插值后数组为：x16=1,0,3,0,2,0,4,0,3,0,5,0,4,0,6,0,5,0,7,0,6,0,8,0,7,0,9,0,8,0,10,0结果分析与图(1)在示例程序中，在x8=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0中反向插入原序列，使原序列变为x316=1,0,8,0,2,0,7,0,3,0,6,0,4,0,5,0,5,0,4,0,6,0,3,0,7,0,2,0,8,0,1,0，再进行DFT变换，观察频谱，对比时域、频域能量和。反向插值后，时域、频域图可以看到，反向插值后，信号频谱有了很大的直流分量，且近乎左右对称。从三维频谱图上可以看出，高频、低频部分实际上是共轭反对称：反向插值后，三维频域图。符合帕塞瓦尔定理，且能量是未插值时的2倍。(2)在x8=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0中插入序列3,0,4,0,5,0,6,0,7,0,8,0,9,0,10,0，使原序列变为x316= 1,0,3,0,2,0,4,0,3,0,5,0,4,0,6,0,5,0,7,0,6,0,8,0,7,0,9,0,8,0,10,0，再进行DFT变换，观察频谱，对比时域、频域能量和。插值后，时域、频域图可以看到，插值后，信号频谱有了很大的直流分量，且近乎左右对称。，符合帕塞瓦尔定理。(3)在x8=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0中正向插入原序列，使原序列分别变为x216= 1,0,1,0,2,0,2,0,3,0,3,0,4,0,4,0,5,0,5,0,6,0,6,0,7,0,7,0,8,0,8,0，再进行DFT变换，观察频谱，对比时域、频域能量和。正向插值后，时域、频域图可以看到，正向插值后，信号频谱有了很大的直流分量，且近乎左右对称。符合帕塞瓦尔定理，且能量是未插值时的2倍。C、插0效应研究原数组：x8=1,0,2,0,3,0,4,0,5,0,6,0,7,8,0示例程序中插0后数组为：x416=1,0,0,0,2,0,0,0,3,0,0,0,4,0,0,0,5,0,0,0,6,0,0,0,7,0,0,0,8,0,0,0插0方式我使用的插0方式为：for(i=0;i16;i=i+3)x4i=COMPLEX(1+i/2,0);x4i+1=COMPLEX(2+i/2,0);x4i+2=COMPLEX(0,0);插0后数组为：x412=1,0,2,0,0,0,3,0,4,0,0,0,5,0,6,0,0,0,7,0,8,0,0,0结果分析与图(1) 在示例程序中，在x8=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0中，每隔一点，插入1个0值，使原序列分别变为x116= 1,0,0,0,2,0,0,0,3,0,0,0,4,0,0,0,5,0,0,0,6,0,0,0,7,0,0,0,8,0,0,0，再进行DFT变换，观察频谱，对比时域、频域能量和。插0前，时域、频域图 插0后，时域、频域图可以看到，插0后的频谱是对原始信号频谱的周期延拓。画出三维图像，可以更直观地看出周期延拓关系： 未插入零/插入一个零后的三维频谱图通过对插零后图像进行DFT运算，可以证明插零后的DFT是原信号DFT的周期延拓。符合帕塞瓦尔定理，且能量与未插值时的相同。（2） 在x8=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0中，每隔两点，插入1个0值，使原序列变为x416= 1,0,2,0,0,0,3,0,4,0,0,0,5,0,6,0,0,0,7,0,8,0,0,0，再进行DFT变换，观察频谱，对比时域、频域能量和。插0后，时域、频域图符合帕塞瓦尔定理源程序：/ 11yy.cpp : Defines the entry point for the console application./#include stdafx.h#includeD:xhclgcyyx_math.cpp#includeD:xhclgcyyx_graph.cppvoid plotgri2(COLORREF gridcolor,COLORREF linecolor,COMPLEX p,int N)int i;HPEN pen1=CreatePen(PS_SOLID,1,gridcolor),oldpen=(HPEN)SelectObject(win3.hdc,pen1);HPEN pen2=CreatePen(PS_SOLID,1,linecolor);for(i=0;iN;i+)line2(i,0,i,abs(pi);SelectObject(win3.hdc,pen2);moveto2(0,p0.r);for(i=0;iN;i+)lineto2(i,abs(pi);SelectObject(win2.hdc,oldpen);DeleteObject(pen1);DeleteObject(pen2);void plotgri3(COLORREF gridcolor,COLORREF linecolor,COMPLEX p,int N)int i;HPEN pen1=CreatePen(PS_SOLID,1,gridcolor),oldpen=(HPEN)SelectObject(win3.hdc,pen1);HPEN pen2=CreatePen(PS_SOLID,1,linecolor);for(i=0;iN;i+)line3(i,0,0,i,pi.r,pi.i);SelectObject(win3.hdc,pen2);moveto3(0,p0.r,p0.i);for(i=0;iN;i+)lineto3(i,pi.r,pi.i);SelectObject(win2.hdc,oldpen);DeleteObject(pen1);DeleteObject(pen2);void main()int i;double sumT,sumF;COMPLEX x8,/1,0,2,0,3,0,4,0,5,0,6,0,7,8,0X8,x213,/=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0,0,0,0,0,0,0,0,0,0,0X216,x316,/=1,0,1,0,2,0,2,0,3,0,3,0,4,0,4,0,5,0,5,0,6,0,6,0,7,0,7,0,8,0,8,0X316,x412,/=1,0,2,0,0,0,3,0,4,0,0,0,5,0,6,0,0,0,7,0,8,0,0,0X416;/给待变换的复数数组赋值：for(i=0;i8;i+)xi=COMPLEX(i+1,0);Xi=COMPLEX(0,0);for(i=0;i8;i+)x2i=COMPLEX(i+1,0);for(i=8;i13;i+)x2i=COMPLEX(0,0);for(i=0;i16;i=i+2)x3i=COMPLEX(1+i/2,0);x3i+1=COMPLEX(1+i/2,0);for(i=0;i16;i=i+3)x4i=COMPLEX(1+i/2,0);x4i+1=COMPLEX(2+i/2,0);x4i+2=COMPLEX(0,0);/第1步：验证dft函数正确性dft(X,x,8,1);for(i=0;i8;i+)printf(X%d=%f+%fn,i,Xi.r,Xi.i);getch();dft(x,X,8,-1);for(i=0;i8;i+)printf(x%d=%f+%fn,i,xi.r,xi.i);getch();for(sumT=0,sumF=0,i=0;i8;i+)sumT=sumT+xi.r*xi.r;sumF=sumF+Xi.r*Xi.r+Xi.i*Xi.i;printf(时域能量和=%f,频域能量和=%fn,sumT,sumF/8.0);window2(函数图形显示,-20,40,20,-20,t,f(t);xy2(BLUE);plotgri2(BLUE,RED,X,8);getch();frame2(win2.xstr,win2.ystr);xy2(BLUE);plotgri2(BLUE,RED,x,8);getch();window3(周期信号频谱图,-1,-12,-12,20,12,12,N,r,i);xyz3(BLUE);plotgri3(BLUE,RED,X,8);getch();frame3();xyz3(BLUE);plotgri3(BLUE,RED,x,8);getch();/第2步：补0效应dft(X2,x2,13,1);for(i=0;i13;i+)printf(X2%d=%f+%fn,i,X2i.r,X2i.i);getch();dft(x2,X2,13,-1);for(i=0;i13;i+)printf(x2%d=%f+%fn,i,x2i.r,x2i.i);getch();for(sumT=0,sumF=0,i=0;i13;i+)sumT=sumT+x2i.r*x2i.r;sumF=sumF+X2i.r*X2i.r+X2i.i*X2i.i;printf(时域能量和=%f,频域能量和=%fn,sumT,sumF/16.0);window2(函数图形显示,-20,40,20,-20,t,f(t);xy2(BLUE);plotgri2(BLUE,RED,X2,13);getch();frame2(win2.xstr,win2.ystr);xy2(BLUE);plotgri2(BLUE,RED,x2,13);getch();window3(周期信号频谱图,-1,-12,-12,20,12,12,N,r,i);xyz3(BLUE);plotgri3(BLUE,RED,X,13);getch();frame3();xyz3(BLUE);plotgri3(BLUE,RED,x,13);getch();/第3步：插值效应dft(X3,x3,16,1);for(i=0;i16;i+)printf(X3%d=%f+%fn,i,X3i.r,X3i.i);getch();dft(x3,X3,16,-1);for(i=0;i16;i+)printf(x3%d=%f+%fn,i,x3i.r,x3i.i);getch();for(sumT=0,sumF=0,i=0;i16;i+)sumT=sumT+x3i.r*x3i.r;sumF=sumF+X3i.r*X3i.r+X3i.i*X3i.i;printf(时域能量和=%f,频域能量和=%fn,sumT,sumF/16.0);window2(函数图形显示,-20,20,20,-20,t,f(t);xy2(BLUE);plotgri2(BLUE,RED,X3,16);getch();frame2(win2.xstr,win2.ystr);xy2(BLUE);plotgri2(BLUE,RED,x3,16);getch();window3(周期信号频谱图,-1,-12,-12,20,12,12,N,r,i);xyz3(BLUE);plotgri3(BLUE,RED,X3,16);getch();frame3();xyz3(BLUE);plotgri3(BLUE,RED,x3,16);getch();/第4步：插0效应dft(X4,x4,12,1);for(i=0;i12;i+)printf(X4%d=%f+%fn,i,X4i.r,X4i.i);getch();dft(x4,X4,12,-1);for(i=0;i12;i+)printf(x4%d=%f+%fn,i,x4i.r,x4i.i);getch();for(sumT=0,sumF=0,i=0;i12;i+)sumT=sumT+x4i.r*x4i.r;sumF=sumF+X4i.r*X4i.r+X4i.i*X4i.i;printf(时域总和=%f,频域总和=%fn,sumT,sumF/12.0);window2(函数图形显示,-20,40,20,-20,t,f(t);xy2(BLUE);plotgri2(BLUE,RED,X4,12);getch();frame2(win2.xstr,win2.ystr);xy2(BLUE);plotgri2(BLUE,RED,x4,12);getch();window3(周期信号频谱图,-1,-12,-12,20,12,12,N,r,i);xyz3(BLUE);plotgri3(BLUE,RED,X4,12);getch();frame3();xyz3(BLUE);plotgri3(BLUE,RED,x4,12);getch(); 训练十二 DFT及抽样定理研究D、给定单频信号抽样1. 对给定信号x(t)=sin(2fct),fc=50,N=264进行抽样，抽样频率分布为100Hz，110Hz，200Hz，230Hz，250Hz。单频信号抽样示意图关于能量泄露的问题： 1.计算时域能量、频域能量。观察是否满足Et=Ef=2|X50|2/N，若抽样频率正确，即无泄漏。 fs=110Hz fs=200Hz Fs=230Hz fs=250Hz 抽样定理指出，若信号x(t)的最高频率为fc,当抽样频率fs2fc时，可由抽样信号x(nTs)完全恢复原信号x(t)。本题中，fs100Hz，符合采样定理条件。但在fs=100Hz时，由于采样点均是原信号x(t)取值为0的点，导致采样信号值始终为0，所以时域能量、频域能量始终为0。其他采样频率时，时域、频域能量始终相等，约为132，符合帕塞瓦尔定理。实验结果：抽样频率为110Hz，200Hz时，没有发生能量泄漏，抽样频率为230Hz，250Hz时，发生能量泄漏。 计算是否发生能量泄漏的依据是，计算频率为原信号频率fc=50Hz处的频谱能量，E=2|X50|2/N。而经过长度为N的DFT变换，对应f=50Hz的k点的计算公式为k=N*fc/fs。当抽样频率为110Hz和200Hz时，k分别等于120和66，因此只需计算第120点和第66点的能量，就可以计算X50处的能量，不会发生泄漏。但当抽样频率为230Hz，250Hz时，k不能取到整数值，只能取到近似值，能量会泄漏到周围频率上，会产生能量泄漏，这种泄露是由于非整周期采样引起的。 2. 对x(t)=sin(2fct),fc=50进行fs=250hz的抽样，然后分别进行250点的DFT和补零至256点的FFT，并计算时域能量、频域能量。观察是否满足Et=Ef=2|X50|2/N，若抽样频率正确，则无泄漏。 250点DFT 256点FFT 由实验结果可知，经过DFT、FFT前后均符合帕塞瓦尔定理，即时域、频域能量一致，均为原始信号能量125。但在用DFT计算X50处的能量时，计算结果与原始信号能量相同，即未发生能量泄漏；在用FFT计算X50处的能量时，计算结果低于原始信号能量，即发生了能量泄漏。 因为取f=50Hz所对应的k值时，k=N*fc/fs。对250点DFT，k=250*50/250=50，计算k=50这点的能量即可；对256点FFT，k=256*50/250，不能取到整数，因此能量会泄漏到周围的频率区间，即发生能量泄漏。 250点DFT 256点FFT 由图可知，DFT在f=50Hz(k=50)的能量即是频谱上的全部能量。而FFT由于补零之后，k不能取到整数，能量一定会泄漏到周围的频率区间，f=50Hz的能量并不是频谱上的全部能量，而是有其他频率分量，所以只计算f=50Hz处的能量时，会低于频谱上的全部能量。结论：所以对于DFT，只要满足k=N*fc/fs可以取到整数，即可实现无能量泄漏。对于FFT，因为需要补零，所以能量一定会泄漏。E、画出单频信号抽样后的频谱图 F、标出频谱图上的单频信号位置 Fs=110Hz未发生能量泄漏 Fs=200Hz未发生能量泄漏 F=50Hz对应k=120 F=50Hz对应k=66 fs=230Hz发生能量泄漏 fs=250Hz发生能量泄漏 F=50Hz 无法对应整数k F=50Hz 无法对应整数k 对X(k)信号进行反变换，再次得到x(n)，与原始信号相对比，发现无论这五种抽样频率下，均能恢复原始信号。说明采样频率fs大于信号最高频率fc时，能不失真地恢复原信号，符合帕塞瓦尔定理。源程序：/ 12yy.cpp : Defines the entry point for the console application./#include stdafx.h#includeD:xhclgcyyx_math.cpp#includeD:xhclgcyyx_graph.cppCOMPLEX x264,Y264,X50,X264,x1264;/1,0,2,0,3,0,4,0,5,0,6,0,7,8,0,X做完DFT之后的,Y是取模之后的DFT,x1为反DFTvoid plotgri2(COLORREF gridcolor,COLORREF linecolor,COMPLEX p,int N)int i;HPEN pen1=CreatePen(PS_SOLID,1,gridcolor),oldpen=(HPEN)SelectObject(win3.hdc,pen1);HPEN pen2=CreatePen(PS_SOLID,1,linecolor);/for(i=0;iN;i+)line2(i,0,i,pi.r);SelectObject(win3.hdc,pen2);for(i=0;iN;i+)moveto2(i,0);lineto2(i,pi.r);SelectObject(win2.hdc,oldpen);DeleteObject(pen1);DeleteObject(pen2);void chouyang(int N,double fs,int fc) /对正弦信号进行抽样int i;double t=0,dt;dt=1.0/fs;for(i=0;iN;i+) t=i*dt;xi=COMPLEX(sin(2*M_PI*fc*t),0);void main(void)int N=264;/数据长度int i,fc=50,fs=110,down,up;/fs可变fs=100hz;fs=110hz;fs=200hz;fs=230hz;fs=250hzdouble sumT,sumF,dt,aa,bb,bb1,EX50;for(i=0;iN;i+)Xi=COMPLEX(0,0);chouyang(N,fs,fc);/给正弦波抽样dft(X,x,264,1);/输出X，输入x，264个数，1：正变换window2(函数图形显示,-10,60,10,-10,i,Xi);/抽样出来的xn显示window2(函数图形显示,-1,2,264,-2,i,xi);xy2(GREEN);/画xy轴。plotgri2(BLUE,RED,x,264);/宏函数工具方式,注意x是形式上的宏参数，它指明后面函数的自变量。getch();frame2();/对做完DFT的Xn取模之后显示for(i=0;i264;i+)Yi=abs(Xi);window2(函数图形显示,-1,100,264,-1,i,Xi);xy2(GREEN);/画xy轴。plotgri2(BLUE,BLUE,Y,264);/宏函数工具方式,注意x是形式上的宏参数，它指明后面函数的自变量。getch();frame2();for(sumT=0,sumF=0,i=0;i264;i+)sumT=sumT+xi.r*xi.r;sumF=sumF+Xi.r*Xi.r+Xi.i*Xi.i;printf(时域能量和=%f,频域能量和=%fn,sumT,sumF/264.0);/求是否有频谱泄露aa=50*N/fs; /aa:50HZ处对应的点；down=floor(aa);/向下取整X50=Ydown;EX50=pow(abs(X50),2)*2/N;printf(X50运算%fn,EX50);getchar();/求反DFT看是否可以恢复出原函数：dft(x1,X,264,-1);window2(函数图形显示,-1,2,264,-2,i,x1i);xy2(GREEN);/画xy轴。plotgri2(BLUE,RED,x1,264);/宏函数工具方式,注意x是形式上的宏参数，它指明后面函数的自变量。getch();frame2();/验证DFT及其频谱泄露情况N=250;fs=250;for(i=0;iN;i+)Xi=COMPLEX(0,0);chouyang(N,fs,fc);/对信号进行抽样，抽样频率为250hz，抽样点数为250点window2(函数图形显示,-1,2,250,-2,i,xi);xy2(RED);/画xy轴。plotgri2(BLUE,BLUE,x,250);/宏函数工具方式,注意x是形式上的宏参数，它指明后面函数的自变量。getch();frame2();dft(X,x,250,1);/做250点DFTfor(i=0;i250;i+)Yi=abs(Xi);window2(函数图形显示,-1,100,250,-1,i,Xi);xy2(RED);/画xy轴。plotgri2(BLUE,BLUE,Y,250);/宏函数工具方式,注意x是形式上的宏参数，它指明后面函数的自变量。getch();frame2();for(sumT=0,sumF=0,i=0;i250;i+)sumT=sumT+xi.r*xi.r;sumF=sumF+Xi.r*Xi.r+Xi.i*Xi.i;printf(时域能量和=%f,频域能量和=%fn,sumT,sumF/250.0);/看是否有泄漏aa=50*N/fs; /aa:50HZ处对应的点；down=floor(aa);/向下取整X50=Ydown;EX50=pow(abs(X50),2)*2/N;printf(X50=%2.fn,X50);printf(X50运算%fn,EX50);getchar();/验证反DFT可不可以恢复原信号dft(x,X,250,-1);window2(函数图形显示,-1,1,250,-1,i,Xi);xy2(RED);/画xy轴。plotgri2(BLUE,BLUE,x,250);/宏函数工具方式,注意x是形式上的宏参数，它指明后面函数的自变量。getch();frame2();/验证FFTN=256;fs=250;chouyang(N,fs,fc);for(i=250;iN;i+)xi=COMPLEX(0,0);for(sumT=0,sumF=0,i=0;i256;i+)sumT=sumT+xi.r*xi.r;/求时域能量window2(函数图形显示,-1,1,256,-1,i,xi);xy2(RED);/画xy轴。plotgri2(BLUE,BLUE,x,256);/宏函数工具方式,注意x是形式上的宏参数，它指明后面函数的自变量。getch();frame2();fft(x,256,1);for(i=0;i256;i+)Yi=abs(xi);window2(函数图形显示,-1,100,256,-100,i,Xi);xy2(RED);/画xy轴。plotgri2(BLUE,BLUE,Y,256);/宏函数工具方式,注意x是形式上的宏参数，它指明后面函数的自变量。getch();frame2();for(sumF=0,i=0;i256;i+)sumF=sumF+xi.r*xi.r+xi.i*xi.i;printf(时域能量和=%f,频域能量和=%fn,sumT,sumF/256.0);aa=50*N/fs; /aa:50HZ处对应的点；down=floor(aa);/向下取整X50=Ydown;EX50=pow(abs(X50),2)*2/N;printf(X50=%2.fn,X50);printf(X50运算%fn,EX50);getchar(); 训练十三 数字滤波器制作1设一阶数字滤波器的系统函数为： 系统输入为以下三种数字序列： （1）矩形序列： （2）三角序列： （3）正弦序列： 2设二阶数字滤波器的系统函数为： 其输入为单位阶跃信号 3、编写通用数字滤波器函数工具。通用数字滤波器函数工具程序：void IIRT(COMPLEX input,COMPLEX output,double a,double b,int N,int Ne)/N为数据点数，Ne为0状态最大负输入，a为Y输出的系数，b为x输入的系数；int i,n;for(n=0;nN;n+)for(i=1;iNe;i+)if(ni) outputn+=ai*0;/零状态输入所有负值都为0else outputn=outputn+ai*output(n-i);for(i=0;iNe;i+) if(ni)outputn+=bi*0;elseoutputn=outputn+bi*input(n-i);源程序：/ 13yy.cpp : Defines the entry point for the console application./#include stdafx.h#includeD:xhclgcyyx_math.cpp#includeD:xhclgcyyx_graph.cppdouble Y40;void plotgri2(COLORREF gridcolor,COLORREF linecolor,double p,int N)int i;HPEN pen1=CreatePen(PS_SOLID,1,gridcolor),oldpen=(HPEN)SelectObject(win3.hdc,pen1);HPEN pen2=CreatePen(PS_SOLID,1,linecolor);for(i=0;iN;i+)line2(i,0,i,pi);SelectObject(win2.hdc,pen2);moveto2(0,p0);SelectObject(win2.hdc,oldpen);DeleteObject(pen1);for(i=0;iN;i+)lineto2(i,pi);DeleteObject(pen2);void f1(double a0,double b0,double b1,double X,int N)/通用一阶数字滤波器函数工具int n;for(n=0;nN;n+)Yn=0; for(n=0;nN;n+)Yn=b0*Xn+b1*Xn-1+a0*Yn-1;void f2(double a0,double a1,double b0,double b1,double b2,double X,int N)/通用二阶数字滤波器函数工具int n;for(n=0;nN;n+)Yn=0; for(n=0;nN;n+)Yn=b0*Xn+b1*Xn-1+b2*Xn-2+a0*Yn-1+a1*Yn-2;int main(int argc, char* argv)int i,n;double X140,X240,X340,X440;/ONE/X1-1=0;for(i=0;i12;i+)X1i=1;for(i=12;i38;i+)X1i=0;window2(X1,-10,-1,30,5,x,y);xy2(BLUE);plotgri2(GREEN,RED,X1,30);getch();f1(0.4,1,0.2,X1,40);plotgri2(GREEN,BLUE,Y,30);getch();/TWO/X2-1=0;for(i=0;i6;i+)X2i=(double)i/5;for(i=6;i11;i+)X2i=(10-(double)i)/5;for(i=11;i38;i+)X2i=0;window2(X2,-10,-1,20,5,x,y);xy2(BLUE);plotgri2(GREEN,RED,X2,30);getch();f1(0.4,1,0.2,X2,40);plotgri2(GREEN,BLUE,Y,30);getch();/THREE/X3-1=0;for(i=0;i18;i+)X3i=sin(double)i*6.28/18);for(i=18;i38;i+)X3i=0;window2(X3,-1,-5,30,5,x,y);xy2(BLUE);plotgri2(GREEN,RED,X3,30);getch();f1(0.4,1,0.2,X3,40);plotgri2(GREEN,BLUE,Y,30);getch();/FOUR/X40=1;for(i=-19;i0;i+)X4i=0;for(i=1;i40;i+)X4i=0;window2(X4,-1,-2,20,5,x,y);xy2(BLUE);plotgri2(GREEN,RED,X4,30);getch();f2(1.5,-0.5,0,1,0,X4,40);plotgri2(GREEN,BLUE,Y,30);getch();return 0;训练十四 IIR数字滤波器设计与实现根据实验指导书要求：(1)令:f1=300hz,f2=400hz,fs=1000hz,a1=3dB,a2=35dB。 设计LP DF(2)令:f1=300hz,f2=400hz,fs=1000hz,a1=3dB,a2=35dB。 设计HP DF(3)令:f1=200hz,f2=300hz,f3=400hz,f4=500hz,fs=2000hz,a1=3dB, a2=40dB。设计BF DF(4) 参数同3，设计BS DF程序：void main()IIR(1,3.0,18.0,2000.,200.,300.,400.,500.);IIR(2,3.0,18.0,2000.,200.,300.,400.,500.);IIR(3,3.0,18.0,2000.,200.,300.,400.,500.);IIR(4,3.0,18.0,2000.,200.,300.,400.,500.);IIRLv(1,3.0,18.0,2000.,200.,300.,400.,500.);/对方波进行低通滤波（自定义函数）;IIRLv2(2,3.0,18.0,2000.,300.,400.,500.,600.);/对多次正弦谐波进行高通滤波（自定义函数）;G、自定指标，并按指标设计，按格式打印出系统H数组转移函数H数组，级连数L（1） 低通 f1=300hz,f2=400hz,fs=1000hz,a1=3dB,a2=35dB。（2）高通 f1=300hz,f2=400hz,fs=1000hz,a1=3dB,a2=35dB。 （3） 带通 f1=200hz,f2=300hz,f3=400hz,f4=500hz,fs=2000hz,a1=3dB,a2=40dB。（4） 带阻 f1=200hz,f2=300hz,f3=400hz,f4=500hz,fs=2000hz,a1=3dB,a2=40dB。H、对H数组分母常数系数归一void JILian_(void)/实现L级级联int i,j,n; window2(函数图形显示,-1,2,200,-2,i,xi);xy2(RED);/画xy轴。plotgri2(