前辈信号处理实验报告.docx
《前辈信号处理实验报告.docx》由会员分享,可在线阅读,更多相关《前辈信号处理实验报告.docx(43页珍藏版)》请在冰豆网上搜索。
前辈信号处理实验报告
信息工程实验报告
08022128陈希爽
实验十九线性卷积计算
实验目的:
1.掌握两短序列线性卷积的直接计算方法
2.掌握DFT计算线性卷积
3.掌握一长序列与一短序列做线性卷积的方法
实验内容:
实验截图:
直接卷积:
Fft卷积:
Fft卷积:
重叠保留:
重叠相加:
实验程序:
#include
voidconvol(double*x,intN,double*h,intM,double*y,intL)
/*y[0]..y[L-1],x[0]..x[N-1],h[0]..h[M-1];general:
L=M+N-1;
as:
0<=k<=M+N-2,0<=n-k<=M-1,0<=k<=N-1
so:
K>=0,K>=n-M+1;
k<=N-1,K<=n;*/
{
intn,k,kk;
for(n=0;nfor(k=(n>=M)?
n-M+1:
0,kk=(N>n)?
n:
N-1,y[n]=0;k<=kk;k++)
y[n]+=x[k]*h[n-k];
}
doubleshow(doublet,dcomplex*y,intL)
{
if(t=0)returnsqrt(y[(int)t].r*y[(int)t].r+y[(int)t].i*y[(int)t].i);
elsereturn0;
}
doubleshow1(doublet,doubley[],intL)
{
if(t>=0&&t<=L)returny[(int)t];
elsereturn0;
}
main()
{
intN=5,M=3,L=10,i,a;/*n,gd=VGA,gm=1,w,key,kv;*/
doublex[8]={1,2,3,4,5},h[8]={3,2,1},y[20];
dcomplexf[8]={1,0,2,0,3,0,4,0,5,0};
dcomplexr[8]={3,0,2,0,1,0,0,0};
initgd(“c:
\xxgc\“);
window2(“functionconvolution(anykeytonext)”,
-1,100,40,0,
“t”,”f”,RED,BLUE);
/***step1***/
clrscr();
printf(“
(1)directlyconvolution...\n”);
convol(h,M,x,N,y,L);/*调用函数直接线卷*/
plotxy2(RED,show1,y,L);
/*for(i=0;igetch();
clearviewport();
/***step2***//*利用FFT圆卷*/
printf(“
(2)fftconvolution\n”);
for(i=M+N-1,L=1;i!
=0;){i=i>>1;L=L<<1;}
fft_ccc(f,L,1);
fft_ccc(r,L,1);
for(i=0;ifft_ccc(f,L,-1);
plotxy2(RED,show,&f,8);
/*for(i=0;igetch();
clearviewport();
/***step3ZhieJieFFT***//*利用FFT圆卷*/
printf(“(3)15datafftconvolution\n”);
{dcomplexx[32]={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};
dcomplexh[32]={3,0,2,0,1,0,};
fft_ccc(x,32,1);
fft_ccc(h,32,1);
for(i=0;i<32;i++)x[i]=cmul(x[i],h[i]);
fft_ccc(x,32,-1);
plotxy2(RED,show,&x,32);
/*for(i=0;i<20;i++)printf(“x[%d]=%f+j%f,\n”,I,x[i].r,x[i].i);*/
getch();
}
clearviewport();
/***step3-2-DieJieSheQuFa***//*重叠舍去*/
clrscr();
printf(“(4)15datafftconvolutioninDieJieSheQuFa\n”);
{dcomplexxx[15+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};
dcomplexh[8]={3,0,2,0,1,0},y[8],yy[32];
intj,yyj=0;
/*printf(“%d===”,18%5);getch();*/
fft_ccc(h,8,1);
for(j=0;j<=23-8;j=j+5){
for(i=0;i<8;i++)y[i]=xx[i+j];
fft_ccc(y,8,1);
for(i=0;i<8;i++)y[i]=cmul(y[i],h[i]);
fft_ccc(y,8,-1);
for(i=3;i<8;i++,yyj++)
{
yy[yyj]=y[i];
}
/***notiewhyi=3?
?
?
****/
/*for(i=3;i<8;i++)printf(“y[%d]=%f+%fj\n”,I,y[i].r,y[i].i);*/
}
plotxy2(RED,show,&yy,18);
getch();
}
clearviewport();
/***step3-2-DieJieXiangJiaFa***//*重叠相加*/
clrscr();
printf(“(4)15datafftconvolutioninDieJieXianJiaFa\n”);
{
dcomplexxx[15+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};
dcomplexh[8]={3,0,2,0,1,0},y[8],yy[32],y1[8];
intj,yyj=0;
fft_ccc(h,8,1);
/***pleaseyouansewerwhybelow23-8****/
for(j=0;j<=23-8;j=j+5){
for(i=0;i<8;i++){
if(i>=5){y[i].r=0.0;y[i].i=0.0;}
elsey[i]=xx[i+j];}
/*for(i=0;i<8;i++)printf(“y#[%d]=%f+%fj,\n”,I,y[i].r,y[i].i);*/
fft_ccc(y,8,1);
for(i=0;i<8;i++)y[i]=cmul(y[i],h[i]);
fft_ccc(y,8,-1);
if(j>=5)
{
for(i=5;i<8;i++)
{
y[i-5]=cadd(y1[i],y[i-5]);
}
}
for(i=0;i<5;i++,yyj++)
{
yy[yyj].i=y[i].i;
yy[yyj].r=y[i].r;
}
for(i=0;i<8;i++)
{y1[i].i=y[i].i;
y1[i].r=y[i].r;}
/*for(i=0;i<8;i++)printf(“y$[%d]=%f+%fj,\n”,I,y[i].r,y[i].i);*/
}
plotxy2(RED,show,&yy,18);
}
getch();
closegraph();
}
实验心得:
基本原理:
这两种方法的原理都是当做N点圆卷时,做圆卷的两序列未补零时其两序列长之和小于等于N-1时,两序列的圆卷结果等于线卷结果
原理说明:
无论是重叠相加还是重叠保留,都是将一次圆卷分为几次完成,而
每次圆卷的序列长=输入序列分段后每段的序列长+系统序列长—1+补零至2的n次幂
换而言之,这两种方法每次做圆卷的两个序列是输入序列补零(或补前一段的值)和系统序列补零至上述算式所算出序列长后的序列。
下面分别详细介绍这两种方法。
重叠相加:
这种方法是将输入序列所分的每段尾部补零至所要求的序列长。
系统函数依旧是补零至要求的序列长。
圆卷满足交换律故以系统函数作为乘数。
当系统函数翻转,周期延拓后,由于其主值序列除n=0处的值都已被翻转至负半轴,因此不难想象其在原主值序列的补零处将出现因周期延拓后而产生的除n=0外的系统主值序列。
这里其与重叠舍去法是相同的。
但由于输入序列尾部补零,故其原补零处出现的非零值没有对圆卷造成干扰,相反由于其在最初几个值做圆卷时,没有算上前一段尾部的值(这部分都被补零了)故需要把前一段的尾部(相应的这部分没有算上后端的头部)与当前段相加,再将各段相接即可。
重叠舍去:
这种方法将前面一段的尾部补至本段的头部以保证每段输入序列的序列长为要求的序列长。
系统函数依旧是补零至要求的序列长。
圆卷满足交换律故以系统函数作为乘数。
当系统函数翻转,周期延拓后,由于其主值序列除n=0处的值都已被翻转至负半轴,因此不难想象其在原主值序列的补零处将出现因周期延拓后而产生的除n=0外的系统主值序列。
也就是说,在原补零处的非零值移出相乘范围之前,其值都将是混叠的值(也就是此时圆卷不等于线卷值)。
那么如前面所写,其非零值的数量共为系统函数序列长减一,故其每段圆卷只需将这些值舍去,再将各段拼接即可。
实验二十IIR数字滤波器设计
实验目的:
掌握用双线性变换法设计IIR数字滤波器的方法掌握由低通到高通,带通,带阻的频率变换方法,掌握巴特沃茨滤波器设计的原理
实验介绍如实验指导书
实验内容:
实验截图:
多项式相乘结果:
123-go之前的数据为3级串联的系统分别的系数,之后为3级系统合并的结果
方波时域:
方波过合并系统的结果:
方波分别过3个系统后的结果:
实验程序:
#include
/*****IIRDigitalFilterDesginProcess***/
doubleshow1(doublet,doubley[],intL)
{
if(t>=0&&t<=L)returny[(int)t];
elsereturn0;
}
intorderBtw(int,double,double,double,double,double,...);
intorderBtw(bandType,db1,db2,fs,f1,f2,f3,f4)
intbandType;
doubledb1,db2,fs,f1,f2,f3,f4;
{
doublepi,omga1,omga2,omga3,omga4,p1,p2,p;
intN;
if(bandType<1
||bandType>4
||db1<0
||db2||f1<0
||f2||fs{printf("inputargumenterrorinorderBtw()\n");
getch();
exit
(1);}
pi=M_PI;/*4*atan(1.0);*/
/*********forlow_passandhigh_pass***********/
omga1=tan(pi*f1/fs);
omga2=tan(pi*f2/fs);
if(bandType<3)p=omga2/omga1;
/*********forband_passandband_stop***********/
else{
if(f3||f4||fs{printf("InputargumenterrorinorderBtw()\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=(p1p1:
p2;
}
elseif(bandType==4)
{
p1=-omga2*(omga4-omga1)/(omga2*omga2-omga4*omga1);
p2=-omga3*(omga4-omga1)/(omga3*omga3-omga4*omga1);
p2=-p2;
p=(p1p1:
p2;}
}
p1=pow(10,db1/10)-1;
p2=pow(10,db2/10)-1;
p=0.5*log10(p2/p1)/log10(p);
N=ceil(p);
returnN;
}
intgetBtw(doubleb[][2],intorder)
/*notice:
pleasedeclaredoubled[N],N=filter_order*/
{
intk,M,L;
doublepi;
pi=4*atan(1.0);
M=order/2;L=M;
if(order!
=M*2){b[M][0]=0;b[M][1]=1;L++;}
for(k=0;kb[k][0]=1;b[k][1]=(-2)*cos((2*k+order+1)*pi/(2*order));}
returnL;
}
doublegetc23(double[2][3],
int,int,double,double,double,double,...);
doublegetc23(c,bandType,N,db1,fs,f1,f2,f3,f4)
intbandType,N;
doublec[2][3],db1,fs,f1,f2,f3,f4;
{
doublefc,pi,*b,namda;
pi=4*atan(1.0);
if(bandType==LOWPASS)
{c[0][0]=1;c[0][1]=-c[0][0];c[0][2]=0.0;
c[1][0]=tan(pi*f1/fs);c[1][1]=c[1][0];c[1][2]=0.0;
fc=f1;}
elseif(bandType==HIGHPASS)
{c[0][0]=tan(pi*f2/fs);c[0][1]=c[0][0];c[0][2]=0.0;
c[1][0]=1;c[1][1]=-1;c[1][2]=0.0;
fc=f2;}
elseif(bandType==BANDPASS)
{doubleU,L;
U=tan(pi*f3/fs);L=tan(pi*f2/fs);
c[0][0]=1+U*L;c[0][1]=2*(U*L-1);c[0][2]=1+U*L;
c[1][0]=U-L;c[1][1]=0;c[1][2]=L-U;
fc=f2;}
elseif(bandType==BANDSTOP)
{doubleU,L;
U=tan(pi*f4/fs);L=tan(pi*f1/fs);
c[1][0]=1+U*L;c[1][1]=2*(U*L-1);c[1][2]=1+U*L;
c[0][0]=U-L;c[0][1]=0;c[0][2]=L-U;
fc=f1;}
namda=pow(10,db1/10)-1;
namda=pow(1/namda,0.5/N);
c[1][0]=namda*c[1][0];
c[1][1]=namda*c[1][1];
c[1][2]=namda*c[1][2];
returnfc;
}
voidBtwAFtoDF(doubleH[][2][5],intL,
doubleb[][2],doublec[2][3]
)
/*jiangs=(a+b/z+c/zz)/(d+e/z+f/zz)dairu
h(s)=1/(m*ss+n*s+1)dexinshizhong
dedao:
h(1/z)=[dd+2de/z+(ee+2df)/zz+2ef/zzz+ff/zzzz]
/{m[aa+2ab/z+(bb+2ac)/zz+2bc/zzz+cc/zzzz]
+[dd+2de/z+(ee+2df)/zz+2ef/zzz+ff/zzzz]
+n[ad+(ae+bd)/z+(cd+af+be)/zz+(bf+ce)/zzz+cf/zzzz]}
=[]/
[maa+dd+nad
+(2mab+2de+nae+nbd)/z
+(mbb+2mac+2df+ee+n(cd+af+be))/zz
+(2mbc+2ef+nbf+nce)/zzz
+(mcc+ff+ncf)/zzzz]
****/
{
inti;
for(i=0;i{H[i][0][0]=c[1][0]*c[1][0];
H[i][0][1]=2*c[1][0]*c[1][1];
H[i][0][2]=c[1][1]*c[1][1]+2*c[1][0]*c[1][2];
H[i][0][3]=2*c[1][1]*c[1][2];
H[i][0][4]=c[1][2]*c[1][2];
H[i][1][0]=b[i][0]*c[0][0]*c[0][0]+c[1][0]*c[1][0]+b[i][1]*c[0][0]*c[1][0];
H[i][1][1]=2*b[i][0]*c[0][0]*c[0][1]+2*c[1][0]*c[1][1]+b[i][1]*(c[0][0]*c[1][1]+c[0][1]*c[1][0]);
H[i][1][2]=b[i][0]*(c[0][1]*c[0][1]+2*c[0][0]*c[0][2])+2*c[1][0]*c[1][2]
+c[1][1]*c[1][1]+b[i][1]*(c[0][2]*c[1][0]+c[0][0]*c[1][2]+c[0][1]*c[1][1]);
H[i][1][3]=2*b[i][0]*c[0][1]*c[0][2]+2*c[1][1]*c[1][2]
+b[i][1]*(c[0][1]*c[1][2]+c[0][2]*c[1][1]);
H[i][1][4]=b[i][0]*c[0][2]*c[0][2]+c[1][2]*c[1][2]+b[i][1]*c[0][2]*c[1][2];
}}
doublepzValue(double*a,intN,doublex)
{
inti;
doublepsumI,psumR;
psumR=a[0];psumI=0.0;
for(i=1;ipsumI=psumI-a[i]*sin(x*i);}
return(sqrt(psumR*psumR+psumI*psumI));
}
doubleBtw20lgHz(doublef,double*fs,doubleH[][2][5],int*L)
{
doublewT,sum;
inti;
wT=8*f*atan(1.0)/fs[0];
for(sum=0,i=0;i{sum=sum+20*log10(pzValue(&(H[i][0][0]),5,wT)/pzValue(&(H[i][1][0]),5,wT));}
returnsum;}
/**********************
(1)*************************************/
/***Todohere,accordingtobookp74(20-1),
drawAnalogButtworthFilter'sfigure***/
doubleBtw20lgHs(doublef,double*C,int*N)
{/***programyourfunctionhere****/
}
/****************************************************************/
voidIIR(int,double,double,double,double,double,...);
voidIIR(bandType,db1,db2,fs,f1,f2,f3,f4)
intbandType;
doubledb1,db2,fs,f1,f2,f3,f4;
{intorder,L,i,NK,j,l,n,k;
doublefc,b[10][2],c[2][3],H[10][2][5],d[3];
doubleZ[2][9],Z1[2][9],Fang[25],yy[25];
order=orderBtw(bandType,db1,db2,fs,f1,f2,f3,f4);
for(i=0;i<9;i++)
{
Z[0][i]=0.0;
Z[1][i]=0.0;
Z1[0][i]=0.0;
Z1[1][i]=0.0;
}
for(i=0;i<5;i++)
{
Fang[i]=5.0;
Fang[i+5]=-5.0;
Fang[i+10]=5.0;
Fang[i+15]=0.0;
Fang[i+20]=0.0;
yy[i]=0.0;
yy[i+5]=0.0;
yy[i+10]=0.0;
yy[i+15]=0.0;
yy[i+20]=0.0;
}
L=getBtw(b,order);
/****fc=****/
getc23(c,bandType,order,db1,fs,f1,f2,f3,f4);
BtwAFtoDF(H,L,b,c);
/***************************
(2)***********************************/
/*todohere,pleaseuseprintf()displayDigitalFilter'scoeaffence*/
/*****************************************************************/
for(i=0;i{
for(j=0;j<2;j