前辈信号处理实验报告

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=n; */int

2、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,gd=VGA,gm

3、=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);/*调用函数直接线卷*/plotxy2

4、(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”);dcomplex x32=1,0,2,0,

5、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-DieJieSheQuFa*/*重叠舍去*/cl

6、rscr();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(i=0;i8;i+)yi=xx

7、i+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 convolution in Die

8、JieXianJiaFan”);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;/* for(i=0;i8;i+)print

9、f(“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次幂换而言之，这两种方法每次做圆卷的两个序列是输入序列补零(或补前一段的值）和系统序列补零至上述算式所算出序列长后的序列。下面分别详细介绍这两种方法。重叠相加：这种方法是将输入序列所分的每段尾部补零至所要求的序列长。系统函数依旧是补零至要求的序列长。圆卷满足交换律故以系统函数作为乘数。当系统函数翻转，周期延拓后，由于其主值序列除n=0处的值都已被翻转至负半

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

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

13、后的结果：实验程序：#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;double pi,omga1,omga2,

14、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*/elseif(f3f2 |f4f

15、3 |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*(omga4-omga1)/(omg

16、a2*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;double pi;pi=4*atan(1.0

17、);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,*b,namda;pi=4*atan

18、(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=1+U*L; c10=U-L;c11

19、=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,int L, double b2,dou

20、ble 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+nae+nbd)/z+(mbb+2mac

21、+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*c10*c12 +c11*c11+

22、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(psumR*psumR+psumI*psum

23、I);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 */double Btw20lgHs(dou

24、ble 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=orderBtw(bandType,db1,d

25、b2,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

26、)*/*to do here,please use printf() display Digital Filters coeaffence */*/for(i=0;iL;i+)for(j=0;j2;j+)for(k=0;k5;k+) printf(%3.5f ,(float)Hijk); printf(n); printf(n); printf(n);getch();printf(123-go!n);for(j=0;j5;j+) Z0j=Z10j=H00j; Z1j=Z11j=H01j;for(j=1;jL;j+) for(n=0;n9;n+) Z0n=0.0; Z1n=0.0; for(n=

27、0;n9;n+)/*多项式相乘*/ for(l=0;l=n;l+) if(n-l)5&l5) Z0n+=Z10l*Hj0n-l; for(n=0;n9;n+) for(l=0;l=n;l+) if(n-l)5&l5) Z1n+=Z11l*Hj1n-l; for(n=0;n9;n+) Z10n=Z0n; Z11n=Z1n; for(i=0;i2;i+)for(j=0;j9;j+) printf(%3.5f ,(float)Zij); printf(nn); getch();int gd=VGA,gm=1;double C;initgraph(&gd,&gm,c:xxgc );window2(lo

28、wpass,-1.,5.,1000.,-60.,hz,db,BLUE,BLUE);plotxy2(GREEN,Btw20lgHz,&fs,H,&L);getch();clearviewport();for(i=0;i25;i+)/*IIR*/ for(j=0;j=0) yyi+=Fangi-j*Z0j; for(j=1;j=0) yyi-=yyi-j*Z1j; yyi=yyi/Z10;window2(lowpass,-50.,10.,50,-10.,hz,db,BLUE,BLUE);plotxy2(RED,show1,Fang,25);getch();clearviewport();windo

29、w2(lowpass,-50.,10.,50,-10.,hz,db,BLUE,BLUE);plotxy2(RED,show1,yy,25);getch();clearviewport();for(k=0;kL;k+)/*IIR*/ for(i=0;i25;i+) yyi=0.0; for(i=0;i25;i+) for(j=0;j=0) yyi+=Fangi-j*Hk0j; for(j=1;j=0) yyi-=yyi-j*Hk1j; yyi=yyi/Hk10; for(i=0;i25;i+) Fangi=yyi; window2(lowpass,-50.,10.,50,-10.,hz,db,B

30、LUE,BLUE);plotxy2(RED,show1,yy,25);getch();/*(3)*/*please modify here draw Analong Buttworth Filter */*C=? */*plotxy2(RED,Btw20lgHs,&C,&order);*/*/setcolor(LIGHTRED);if(bandType=LOWPASS) rectangle2(0.0,0.0,f1,-db1);rectangle2(0.0,0.0,f2,-db2);else if(bandType=HIGHPASS) rectangle2(0.0,.0,f2,-db1);rec

31、tangle2(0.0,0.0,f1,-db2);else if(bandType=BANDSTOP) rectangle2(f1,0.0,f4,-db1);rectangle2(f2,0.0,f3,-db2);else if(bandType=BANDPASS) rectangle2(f2,0.0,f3,-db1);rectangle2(f1,0.0,f4,-db2);getch();closegraph();main()IIR(1,3.0,18.0,2000.,200.,300.,400.,500.);/*IIR(2,3.0,18.0,2000.,200.,300.,400.,500.);

32、*/*IIR(3,3.0,18.0,2000.,200.,300.,400.,500.);*/*IIR(4,3.0,18.0,2000.,200.,300.,400.,500.);*/getch();/*(4)*/*to do here, Realize IIR filter */*/实验心得：本实验利用分别过三个系统函数和过一个总系统函数对比以互相验证。红色的第一段为多项式相乘代码，后两段为IIR计算，后者是分别过三个系统函数本实验首先需要编写多项式相乘的函数，原理是线卷，多项式相乘就是线卷。两个多项式各项系数就是要线卷的两个序列，按升幂排列。之后需要编写IIR滤波器的通用函数。步骤是，首先

33、需要获得IIR系统函数的分子分母阶次，而后为分子分母分别编写循环计算输入序列与分子各项乘积后求和；输出序列与分母各项相乘后求和，再求总和，由于例程所求出的系统函数没有归一化，故需要再将所求输出除以分母的第一项系数，以得出每项输出。实验二十一 FIR数字滤波器设计实验目的：1、掌握用傅里叶级数法设计FIR滤波器的方法 2、了解窗函数在信号处理中的作用实验内容：实验截图：低通幅度（方型窗）：低通相位（方型窗）：带通幅度（方型窗）：带通相位（方型窗）：高通幅度（方型窗）：高通相位（方型窗）：带阻幅度（方型窗）：带阻相位（方型窗）：带通幅度（汉明窗）：带通相位（汉明窗）：高通幅度（汉明窗）：高通相位（

34、汉明窗）：带阻幅度（汉明窗）：带阻相位（汉明窗）：方波时域：方波频域：低通频域：方波过低通(重叠相加法)：方波过低通(fft直接圆卷)：方波过低通后频域：实验程序：#include#define N1 50 /*N1(L1=2x)M1*/#define M1 20#define L1 32double 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;main()initgd(c:xxgc );/*FIR experiment

35、of lowpass*/int flag,N=21;double b21,fl=200,fs=1000;dsnFIR(b,21,LOWPASS,RECTANG,fs,fl);flag=0;window2(LowPassAmp,-fs/2,2,fs/2,-2, hz,|H|, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();flag=1;window2(LowPassPhase,-fs/2,5,fs/2,-5, hz,FI, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch(

36、);/*FIR experiment of bandpass*/int flag,N=21;double b21,fl=150,fh=350,fs=1000;dsnFIR(b,21,BANDPASS,RECTANG,fs,fl,fh);flag=0;window2(BandPassAmp,-fs/2,2,fs/2,-2, hz,|H|, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();flag=1;window2(BandPassPhase,-fs/2,5,fs/2,-5, hz,FI, RED,BLUE);plotxy2(GRE

37、EN,drawFIR,&fs,b,&N,&flag);getch();/*FIR experiment of highpass,To DO here!*/int flag,N=21;double b21,fl=100,fh=250,fs=1000;dsnFIR(b,21,HIGHPASS,RECTANG,fs,fh);flag=0;window2(BandPassAmp,-fs/2,2,fs/2,-2, hz,|H|, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();flag=1;window2(BandPassPhase,-fs

38、/2,5,fs/2,-5, hz,FI, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();/*FIR experiment of bandstop,To Do here!*/int flag,N=21;double b21,fl=150,fh=350,fs=1000;dsnFIR(b,21,BANDSTOP,RECTANG,fs,fl,fh);flag=0;window2(BandPassAmp,-fs/2,2,fs/2,-2, hz,|H|, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&f

39、lag);getch();flag=1;window2(BandPassPhase,-fs/2,5,fs/2,-5, hz,FI, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();int flag,N=21;double b21,fl=200,fs=1000;dsnFIR(b,21,LOWPASS,HAMMING,fs,fl);flag=0;window2(LowPassAmp,-fs/2,2,fs/2,-2, hz,|H|, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getc

40、h();flag=1;window2(LowPassPhase,-fs/2,5,fs/2,-5, hz,FI, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();/*FIR experiment of bandpass*/int flag,N=21;double b21,fl=150,fh=350,fs=1000;dsnFIR(b,21,BANDPASS,HAMMING,fs,fl,fh);flag=0;window2(BandPassAmp,-fs/2,2,fs/2,-2, hz,|H|, RED,BLUE);plotxy2(GR

41、EEN,drawFIR,&fs,b,&N,&flag);getch();flag=1;window2(BandPassPhase,-fs/2,5,fs/2,-5, hz,FI, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();/*FIR experiment of highpass,To DO here!*/int flag,N=21;double b21,fl=100,fh=250,fs=1000;dsnFIR(b,21,HIGHPASS,HAMMING,fs,fh);flag=0;window2(BandPassAmp,-fs

42、/2,2,fs/2,-2, hz,|H|, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();flag=1;window2(BandPassPhase,-fs/2,5,fs/2,-5, hz,FI, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();/*FIR experiment of bandstop,To Do here!*/int flag,N=21;double b21,fl=150,fh=350,fs=1000;dsnFIR(b,21,BANDSTOP,HAMMING,fs,fl,fh);flag=0;window2(BandPassAmp,-fs/2,2,fs/2,-2, hz,|H|, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();flag=1;window2(BandPassPhase,-fs/2,5,fs/2,-5, hz,FI, RED,BLUE);plotxy2(GREEN,drawFIR,&fs,b,&N,&flag);getch();/*Realize FIR as convloution way of experiment 19t