前辈信号处理实验报告.docx

1、前辈信号处理实验报告信息工程实验报告 08022128 陈希爽实验十九 线性卷积计算实验目的 ：1.掌握两短序列线性卷积的直接计算方法 2.掌握DFT计算线性卷积 3.掌握一长序列与一短序列做线性卷积的方法实验内容：实验截图：直接卷积：Fft卷积：Fft卷积：重叠保留：重叠相加：实验程序：#includevoid convol(double *x,int N,double *h,int M,double *y,int L)/*y0.yL-1,x0.xN-1,h0.hM-1;general:L=M+N-1;as:0=k=M+N-2,0=n-k=M-1,0=k=0,K=n-M+1; k=N-1,K

2、=n; */int n,k,kk;for(n=0;n=M)?n-M+1:0,kk=(Nn)?n:N-1,yn=0;k=kk;k+) yn+=xk*hn-k;double show(double t,dcomplex *y,int L)if(t=0) return sqrt(y(int)t.r*y(int)t.r+y(int)t.i*y(int)t.i);else return 0;double show1(double t,double y,int L)if(t=0&t=L) return y(int)t;else return 0;main()int N=5,M=3,L=10,i,a;/*n

3、,gd=VGA,gm=1,w,key,kv;*/double x8=1,2,3,4,5,h8=3,2,1,y20;dcomplex f8=1,0,2,0,3,0,4,0,5,0;dcomplex r8=3,0,2,0,1,0,0,0;initgd(“c:xxgc “);window2(“function convolution(any key to next)”,-1,100,40,0,“t”,”f”,RED,BLUE);/*step 1 */clrscr();printf(“(1) directly convolution. n”);convol(h,M,x,N,y,L); /*调用函数直接

4、线卷*/plotxy2(RED,show1,y,L); /*for(i=0;i1;L=L1;fft_ccc(f,L,1);fft_ccc(r,L,1);for(i=0;iL;i+)fi=cmul(fi,ri);fft_ccc(f,L,-1);plotxy2(RED,show,&f,8);/*for(i=0;iL;i+)printf(“f%d=%f+j%f,n”,I,fi.r,fi.i);*/getch();clearviewport();/*step 3 ZhieJieFFT*/ /*利用FFT圆卷*/printf(“(3) 15 data fft convolutionn”);dcomple

5、x x32=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0,9,0,10,0,11,0,12,0,13,0,14,0,15,0;dcomplex h32=3,0,2,0,1,0,;fft_ccc(x,32,1);fft_ccc(h,32,1);for(i=0;i32;i+)xi=cmul(xi,hi);fft_ccc(x,32,-1);plotxy2(RED,show,&x,32);/*for(i=0;i20;i+)printf(“x%d=%f+j%f,n”,I,xi.r,xi.i);*/getch();clearviewport();/*step 3-2-DieJieSheQ

6、uFa*/ /*重叠舍去*/clrscr();printf(“(4) 15 data fft convolution in DieJieSheQuFan”);dcomplex xx15+3+(8-(3+15)%5)=0,0,0,0,0,0,1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0,9,0,10,0,11,0,12,0,13,0,14,0,15,0;dcomplex h8=3,0,2,0,1,0,y8,yy32;int j,yyj=0;/*printf(“%d=”,18%5);getch();*/fft_ccc(h,8,1);for(j=0;j=23-8;j=j+5)for

7、(i=0;i8;i+)yi=xxi+j;fft_ccc(y,8,1);for(i=0;i8;i+)yi=cmul(yi,hi);fft_ccc(y,8,-1);for(i=3;i8;i+,yyj+) yyyyj=yi;/*notie why i=3?*/* for(i=3;i8;i+)printf(“y%d=%f+%fjn”,I,yi.r,yi.i);*/plotxy2(RED,show,&yy,18);getch();clearviewport();/*step 3-2-DieJieXiangJiaFa */ /*重叠相加*/clrscr();printf(“(4) 15 data fft

8、convolution in DieJieXianJiaFan”);dcomplex xx15+3+(8-(3+15)%5)=1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0,9,0,10,0,11,0,12,0,13,0,14,0,15,0;dcomplex h8=3,0,2,0,1,0,y8,yy32,y18;int j,yyj=0;fft_ccc(h,8,1);/*please you ansewer why below 23-8*/for(j=0;j=23-8;j=j+5) for(i=0;i=5)yi.r=0.0;yi.i=0.0; else yi=xxi+j;/* f

9、or(i=0;i8;i+)printf(“y#%d=%f+%fj,n”,I,yi.r,yi.i); */fft_ccc(y,8,1);for(i=0;i=5) for(i=5;i8;i+) yi-5=cadd(y1i,yi-5); for(i=0;i5;i+,yyj+)yyyyj.i=yi.i;yyyyj.r=yi.r;for(i=0;i8;i+)y1i.i=yi.i;y1i.r=yi.r;/* for(i=0;i8;i+)printf(“y$%d=%f+%fj,n”,I,yi.r,yi.i);*/plotxy2(RED,show,&yy,18);getch(); closegraph();实

10、验心得：基本原理：这两种方法的原理都是当做N点圆卷时，做圆卷的两序列未补零时其两序列长之和小于等于N-1时，两序列的圆卷结果等于线卷结果原理说明：无论是重叠相加还是重叠保留，都是将一次圆卷分为几次完成，而每次圆卷的序列长= 输入序列分段后每段的序列长+系统序列长1+补零至2的n次幂换而言之，这两种方法每次做圆卷的两个序列是输入序列补零(或补前一段的值）和系统序列补零至上述算式所算出序列长后的序列。下面分别详细介绍这两种方法。重叠相加：这种方法是将输入序列所分的每段尾部补零至所要求的序列长。系统函数依旧是补零至要求的序列长。圆卷满足交换律故以系统函数作为乘数。当系统函数翻转，周期延拓后，由于其主

11、值序列除n=0处的值都已被翻转至负半轴，因此不难想象其在原主值序列的补零处将出现因周期延拓后而产生的除n=0外的系统主值序列。这里其与重叠舍去法是相同的。但由于输入序列尾部补零，故其原补零处出现的非零值没有对圆卷造成干扰，相反由于其在最初几个值做圆卷时，没有算上前一段尾部的值(这部分都被补零了)故需要把前一段的尾部(相应的这部分没有算上后端的头部)与当前段相加，再将各段相接即可。重叠舍去：这种方法将前面一段的尾部补至本段的头部以保证每段输入序列的序列长为要求的序列长。系统函数依旧是补零至要求的序列长。圆卷满足交换律故以系统函数作为乘数。当系统函数翻转，周期延拓后，由于其主值序列除n=0处的值都

12、已被翻转至负半轴，因此不难想象其在原主值序列的补零处将出现因周期延拓后而产生的除n=0外的系统主值序列。也就是说，在原补零处的非零值移出相乘范围之前，其值都将是混叠的值(也就是此时圆卷不等于线卷值)。那么如前面所写，其非零值的数量共为系统函数序列长减一，故其每段圆卷只需将这些值舍去，再将各段拼接即可。实验二十IIR数字滤波器设计实验目的：掌握用双线性变换法设计IIR数字滤波器的方法掌握由低通到高通，带通，带阻的频率变换方法，掌握巴特沃茨滤波器设计的原理实验介绍如实验指导书实验内容：实验截图：多项式相乘结果：123-go之前的数据为3级串联的系统分别的系数，之后为3级系统合并的结果方波时域：方波

13、过合并系统的结果：方波分别过3个系统后的结果：实验程序：#include/*IIR Digital Filter Desgin Process*/double show1(double t,double y,int L) if(t=0&t=L) return y(int)t; else return 0;int orderBtw(int,double,double,double,double,double,.);int orderBtw(bandType,db1,db2,fs,f1,f2,f3,f4)int bandType;double db1,db2,fs,f1,f2,f3,f4;doub

14、le pi,omga1,omga2,omga3,omga4,p1,p2,p;int N;if(bandType4 |db10 |db2db1 |f10 |f2f1 |fsf2) printf(input argument error in orderBtw()n); getch(); exit(1);pi=M_PI;/*4*atan(1.0);*/*for low_pass and high_pass*/omga1=tan(pi*f1/fs);omga2=tan(pi*f2/fs);if(bandType3)p=omga2/omga1;/*for band_pass and band_stop

15、*/elseif(f3f2 |f4f3 |fsf4) printf(Input argument error in orderBtw()n); getch();exit(1);omga3=tan(pi*f3/fs);omga4=tan(pi*f4/fs); if(bandType=3) p1=(omga1*omga1-omga3*omga2)/(omga1*(omga3-omga2); p2=(omga4*omga4-omga3*omga2)/(omga4*(omga3-omga2); p1=-p1; p=(p1p2)?p1:p2; else if(bandType=4) p1=-omga2*

16、(omga4-omga1)/(omga2*omga2-omga4*omga1); p2=-omga3*(omga4-omga1)/(omga3*omga3-omga4*omga1); p2=-p2; p=(p1p2)?p1:p2; p1=pow(10,db1/10)-1;p2=pow(10,db2/10)-1;p=0.5*log10(p2/p1)/log10(p);N=ceil(p);return N;int getBtw(double b2,int order)/*notice:please declare double dN,N=filter_order */int k,M,L;doubl

17、e pi;pi=4*atan(1.0);M=order/2;L=M;if(order!=M*2)bM0=0;bM1=1;L+;for(k=0;kM;k+) bk0=1;bk1=(-2)*cos(2*k+order+1)*pi/(2*order);return L;double getc23(double 23, int,int,double,double,double,double,.);double getc23(c,bandType,N,db1,fs,f1,f2,f3,f4)int bandType,N;double c23,db1,fs,f1,f2,f3,f4;double fc,pi,

18、*b,namda;pi=4*atan(1.0);if(bandType=LOWPASS) c00=1;c01=-c00;c02=0.0; c10=tan(pi*f1/fs);c11=c10;c12=0.0; fc=f1;else if(bandType=HIGHPASS) c00=tan(pi*f2/fs);c01=c00;c02=0.0; c10=1;c11=-1;c12=0.0; fc=f2;else if(bandType=BANDPASS) double U,L; U=tan(pi*f3/fs);L=tan(pi*f2/fs); c00=1+U*L;c01=2*(U*L-1);c02=

19、1+U*L; c10=U-L;c11=0;c12=L-U; fc=f2;else if(bandType=BANDSTOP) double U,L; U=tan(pi*f4/fs);L=tan(pi*f1/fs); c10=1+U*L;c11=2*(U*L-1);c12=1+U*L; c00=U-L;c01=0;c02=L-U; fc=f1;namda=pow(10,db1/10)-1;namda=pow(1/namda,0.5/N);c10=namda*c10;c11=namda*c11;c12=namda*c12;return fc;void BtwAFtoDF(double H25,in

20、t L, double b2,double c23 )/*jiang s=(a+b/z+c/zz)/(d+e/z+f/zz) dairuh(s)=1/(m*ss+n*s+1) de xinshi zhongdedao:h(1/z)=dd+2de/z+(ee+2df)/zz+2ef/zzz+ff/zzzz/maa+2ab/z+(bb+2ac)/zz+2bc/zzz+cc/zzzz+dd+2de/z+(ee+2df)/zz+2ef/zzz+ff/zzzz+nad+(ae+bd)/z+(cd+af+be)/zz+(bf+ce)/zzz+cf/zzzz=/maa+dd+nad+(2mab+2de+na

21、e+nbd)/z+(mbb+2mac+2df+ee+n(cd+af+be)/zz+(2mbc+2ef+nbf+nce)/zzz+(mcc+ff+ncf)/zzzz*/int i;for(i=0;iL;i+)Hi00=c10*c10;Hi01=2*c10*c11;Hi02=c11*c11+2*c10*c12;Hi03=2*c11*c12;Hi04=c12*c12;Hi10=bi0*c00*c00+c10*c10+bi1*c00*c10;Hi11=2*bi0*c00*c01+2*c10*c11+bi1*(c00*c11+c01*c10);Hi12=bi0*(c01*c01+2*c00*c02)+2

22、*c10*c12 +c11*c11+bi1*(c02*c10+c00*c12+c01*c11);Hi13=2*bi0*c01*c02+2*c11*c12 +bi1*(c01*c12+c02*c11);Hi14=bi0*c02*c02+c12*c12+bi1*c02*c12;double pzValue(double *a,int N,double x)int i;double psumI,psumR;psumR=a0;psumI=0.0;for(i=1;iN;i+)psumR=psumR+ai*cos(x*i); psumI=psumI-ai*sin(x*i);return(sqrt(psum

23、R*psumR+psumI*psumI);double Btw20lgHz(double f,double *fs,double H25,int *L)double wT,sum;int i;wT=8*f*atan(1.0)/fs0;for(sum=0,i=0;iL0;i+) sum=sum+20*log10(pzValue(&(Hi00),5,wT)/pzValue(&(Hi10),5,wT);return sum;/*(1)*/*To do here,according to book p74(20-1), draw Analog Buttworth Filters figure */do

24、uble Btw20lgHs(double f,double *C,int *N)/*program your function here*/*/void IIR(int,double,double,double,double,double,.);void IIR(bandType,db1,db2,fs,f1,f2,f3,f4)int bandType;double db1,db2,fs,f1,f2,f3,f4;int order,L,i,NK,j,l,n,k;double fc,b102,c23,H1025,d3;double Z29,Z129,Fang25,yy25;order=order

25、Btw(bandType,db1,db2,fs,f1,f2,f3,f4);for(i=0;i9;i+) Z0i=0.0; Z1i=0.0; Z10i=0.0; Z11i=0.0;for(i=0;i5;i+) Fangi=5.0; Fangi+5=-5.0; Fangi+10=5.0; Fangi+15=0.0; Fangi+20=0.0; yyi=0.0; yyi+5=0.0; yyi+10=0.0; yyi+15=0.0; yyi+20=0.0;L=getBtw(b,order);/*fc=*/getc23(c,bandType,order,db1,fs,f1,f2,f3,f4);BtwAFtoDF(H,L,b,c);/*(2)*/*to do here,please use printf() display Digital Filters coeaffence */*/for(i=0;iL;i+)for(j=0;j2;j