课程设计：信号分析与处理C语言编程

勘查技术课程设计：信号分析与处理基础（西南石油大学---资源与环境学院）对于勘查技术与工程专业的学生来说，《信号分析与处理基础》是一门专业基础课，我是2010级的，我们是在大三第一学期上的，这门课数学与物理知识要求比较高，不过一开认真仔细学的话，也会学的好的，起码要比那空洞、生奥、蛋疼《弹性波动力学》好学些。随着课程的结束，《信号分析与处理基础》的课程设计也随之而来，我们是老师布置了4个题目，分单号与双号各自做2道，我是单号，做的是滤波与相关。这次课程设计，注意考验大家的编程能力，目前我们学过得就只有C语言，可以用Fortran,Matlab等等，Matlab可以现学现用，上手快。但是大家也可以挑战下自己的C语言，提高下自己的编程能力，这是一次很好地机会，真正实用的时刻。我就是用C语言编的。其余2题，我也把程序与结果图收集到了这里，以供学弟、学妹们参考之用！我的QQ：593066480，有什么不懂的，或者好的见教，欢迎来信息交流！题目如下：滤波已知原始地震记录x(t),要求：设计滤波器，消除x(t)中10Hz以下，80Hz以上的干扰信号。建议参数：A1=1，A2=0.8，A3=0.5f1=25Hz，f2=45Hz，f3=5Hz，f4=80取样点数：N=200抽样间隔：Δ=0.004，ß=100eq\o\ac(○,1)时域滤波：由(h(t)为滤波因子)建议参数：第一参数：ff1=20Hz，ff2=50Hz第二参数：ff1=25Hz，ff2=45Hz抽样点数：M=60抽样间隔：Δ=0.004，，单独计算：要求：画出x(t),h(t),y(t)图形，为了分析方便，也可以画出有效波s(t),干扰波n(t)及其频谱进行分析，如下图:最后就是答辩，老师问问题，学生回答。主要注意几点就行了熟悉课本滤波部分知识；第二参数要比第一参数滤波效果好，因为门第一参数开大了，进来的干扰波也多了，从第一参数：Fy、第二参数Fy图形上可以看出来，干扰波频谱被压小了。第二参数压制了干扰波，突显了有效波，所有好。eq\o\ac(○,2)频域滤波由公式：对x(t)进行补0，28或者29总之必须是2的次方，因为要用到FFT公式与IFFT公式，进行FFT变换得到X（k）；H（k）的求法：H（k）也必须与X（k）点数相同，单门：20,50双门：20,50,=N-,=N-,3,在用IFFT反变换得y(t),存在实部与虚部，需要分析与处理要求：画出H(k)、X(k)、Y(k)图形，并且分析X(k)、Y(k)的区别，还有开单双门的区别与差异，时域与频域滤波谁好、为什么？分析：从单双门实部图形看出，与有效波是完全一样的，但是幅值变大，这体现了滤波突显有效波的特性；从图形看出虚部对结果没有影响；单双门虚部完全不同，这是由于开单双门效果不同引起的。单门虚部变化大，而且幅值与起伏变化也大，而双门幅值很小，小到可以忽略不计。按理说反变化IFFT后应该只有实部，没有虚部，至于为什么会产生虚部，希望读者自己下去研究下，希望大家相互交流、交流！第一个注意问题：单门与双门谁好？答案当然是：双门好。因为原始有效波x(t)是实数信号，对应的频谱是偶对称的，单门时，只是让一部分通过了；而双门则全部通过，肯定效果要好！而且实部波形，明显看出双门是干扰部分起伏比单门时要平缓，小得多了。第二个注意问题：时域与频域哪个好？频域要好，首先从两者波形上可以看出，其次就是频域滤波，只让有效波通过，而之外的完全被滤掉，可谓是真正的理想状态。第三个注意问题：频域滤波y(t)实部图200点以后，为什么有波形起伏？因为由于时域离散，必将导致频域周期化，这是由于频域周期化的结果。相关建议参数：1，30Hz，抽样点数：M=1000.8，50Hz，抽样点数：M=150抽样间隔：Δ=0.004s要求：画出，，的图形并分析验证书上自互相关性质的正确性！快速褶积建议参数：1，25Hz，抽样点数：M=2500.7，55Hz，抽样点数：M=200抽样间隔：Δ=0.004s地震记录的生成和频谱分析地震记录的生成和频谱分析、信噪比计算以及补0对频谱的影响，给定地震子波的数学表达式：和反射系数序列：0.2，0.4，0.15，0.50.35，0.1，0.2产生一个含有随机干扰信号的地震信号：（其中n(t)可以用-0.4—+0.4之间的随机数代替）要求：制作合成地震记录x(t),并对b(t)和s(t)做频谱分析（用FFT），其中：b(t)（N=41,Δ=0.004s，f=35Hz，α=100，A1=1）反射序列和随机干扰N=200；计算x(t)信噪比，改变n(t)(用-0.8—+0.8之间的随机数代替)的大小再计算。信噪比：分别对b(t)后边、前边、中间补0，计算补0前后频谱的变化及补0多少对频谱的影响。附带程序：最麻烦的就是编程序了，开始的时候，是很麻烦，不过只有去啃，就一定会有收获，比如：我开始编褶积程序的时候，整理了几天，上网查，翻图书，后来突然明白，靠公式就可以编出来了。只是FFT需要些功夫,其余都很快！书上介绍的是基2FFT，这个代码网上到处都有，想提升自己的编程能力，就去尝试编基4FFT，以及考虑基nFFT吧，祝大家好运。1、时域滤波程序代码#include<stdio.h>#include<math.h>#include"conv.cpp"#include"dft.cpp"#include"gyh.cpp"#defineT200//************T表示x(t)点数，R表示n(t)点数#defineR60voidmain(){ inti,j,B=100;//*********************B表示贝塔的值 FILE*fp; doubleA1=1.0,A2=0.8,A3=0.5;//********产生x(t) doublef1=25,f2=45,f3=5,f4=80;//******频率取值 doubledata=0.004;//******************取样点数 doubles[T],n[T],x[T],h[R],y[T+R-1]; doublesi[T]={0},s1[T]={0},s1i[T]={0},Fs[T]={0}; doubleni[T]={0},n1[T]={0},n1i[T]={0},Fn[T]={0}; doublexi[T]={0},x1[T]={0},x1i[T]={0},Fx[T]={0}; doublehi[T]={0},h1[T]={0},h1i[T]={0},Fh[T]={0}; doubley1[T+R-1]={0},yb[T+R-1]={0},yb1[T+R-1]={0},Fy[T+R-1]={0}; for(i=0;i<T;i++) { s[i]=A1*exp(-B*pow((i*data),2))*sin(2*PI*f1*i*data)+A2*exp(-B*pow((i*data),2))*sin(2*PI*f2*i*data); n[i]=A3*(sin(2*PI*f3*i*data)+cos(2*PI*f4*i*data)); x[i]=s[i]+n[i]; }//对s(t)进行频谱分析，用DFTdft(s,si,s1,s1i,Fs,T,1); dft(n,ni,n1,n1i,Fn,T,1); dft(x,ni,x1,x1i,Fx,T,1);//统一导出 fp=fopen("D:\\sh\\time1\\snxFsFnFx1.txt","w");for(i=0;i<T;i++) fprintf(fp,"%f\t%f\t%f\t%f\t%f\t%f\n",s[i],n[i],x[i],Fs[i],Fn[i],Fx[i]);fclose(fp);//产生h(t) floatff1=20,ff2=50;//*****************可以修改的h(t)参数 doublew1,w2,dataw,w0; w1=2*PI*ff1,w2=2*PI*ff2; dataw=(w2-w1)/2,w0=(w2+w1)/2; h[0]=2*dataw/PI; for(j=1;j<R;j++) { h[j]=(2.0/(PI*j*data))*cos(w0*j*data)*sin(dataw*j*data); } dft(h,hi,h1,h1i,Fh,R,1); gyh(Fh,R); fp=fopen("D:\\sh\\time1\\hFh.txt","w"); for(j=0;j<R;j++) fprintf(fp,"%f\t%f\n",h[j],Fh[j]); fclose(fp);//褶积滤波得到y(t)con(x,h,y,T,R);//对y(t)作傅里叶变换dft(y,y1,yb,yb1,Fy,T+R-1,1);//导出y及频谱Y(k)fp=fopen("D:\\sh\\time1\\yFy1.txt","w");for(i=0;i<T+R-1;i++) fprintf(fp,"%f\t%f\n",y[i],Fy[i]);fclose(fp);printf("\nover\n\n");}频域滤波程序代码：#include<stdio.h>#include<math.h>#include"fft.cpp"#include"ifft.cpp"#defineG256voidmain(){ inti,B=100; FILE*fp;//定义指针文件 doubleA1=1.0,A2=0.8,A3=0.5;//产生x(t) doublef1=25,f2=45,f3=5,f4=80;//参数 doubledata=0.004;//抽样间隔 doubles[200],n[200]; doublex[G]={0},x1[G]={0},h[G]={0},h1[G]={0},y[G]={0},y1[G]={0}; doubleFx[G]={0},Fh1[G]={0},Fh2[G]={0},Fy1[G]={0},Fy2[G]={0}; doubleY1[G]={0},Y1i[G]={0},Y2[G]={0},Y2i[G]={0}; for(i=0;i<200;i++) { s[i]=A1*exp(-B*pow((i*data),2))*sin(2*PI*f1*i*data)+A2*exp(-B*pow((i*data),2))*sin(2*PI*f2*i*data); n[i]=A3*(sin(2*PI*f3*i*data)+cos(2*PI*f4*i*data)); x[i]=s[i]+n[i]; } //导出x[t] fp=fopen("D:\\sh\\py\\x.txt","w"); for(i=0;i<G;i++) fprintf(fp,"%f\n",x[i]); fclose(fp); //x[t]频谱计算fft(x,x1,Fx,G,1);//导出X(k) fp=fopen("D:\\sh\\py\\Fx.txt","w"); for(i=0;i<G;i++) fprintf(fp,"%f\n",Fx[i]); fclose(fp); //处理h[t]doubledataf;doublem1,m2,m3,m4;dataf=1.0/(data*G);m1=20.0/dataf,m2=50.0/dataf;m3=G-m2,m4=G-m1;//计算,导出单门for(i=0;i<G;i++){ if(i>=int(m1)&&i<=int(m2)) Fh1[i]=1; elseFh1[i]=0;} //计算,导出双门 for(i=0;i<G;i++) { if((i>=int(m1)&&i<=int(m2))||(i>=int(m3)&&i<=int(m4))) Fh2[i]=1; elseFh2[i]=0; } fp=fopen("D:\\sh\\py\\Fh12.txt","w");for(i=0;i<256;i++) fprintf(fp,"%f\t%f\n",Fh1[i],Fh2[i]); fclose(fp); for(i=0;i<G;i++)//单门滤波 { Y1[i]=x[i]*Fh1[i]; Y1i[i]=x1[i]*Fh1[i]; Fy1[i]=sqrt(pow(Y1[i],2)+pow(Y1i[i],2)); } for(i=0;i<G;i++)//开双门 { Y2[i]=x[i]*Fh2[i]; Y2i[i]=x1[i]*Fh2[i]; Fy2[i]=sqrt(pow(Y2[i],2)+pow(Y2i[i],2));} fp=fopen("D:\\sh\\py\\Fy12.txt","w");//导出单双门for(i=0;i<G;i++) fprintf(fp,"%f\t%f\n",Fy1[i],Fy2[i]); fclose(fp);//Y[k]作反变换IFFT,并且导出实部，虚步 ifft(Y1,Y1i,G,-1); ifft(Y2,Y2i,G,-1); fp=fopen("D:\\sh\\py\\yt12.txt","w");for(i=0;i<G;i++) fprintf(fp,"%f\t%f\t%f\t%f\n",Y1[i],Y1i[i],Y2[i],Y2i[i]); fclose(fp);printf("\nover\n\n");}相关程序代码#include<stdio.h>#include<math.h>#include"correl.cpp"#definePI3.14159265voidmain(){ inti,j; doubleA1=1.0,A2=0.8; doublef1=30,f2=50; doubledata=0.004,p=100; doublex[100]={0},y[150]={0},xy[249]={0},yx[249]={0}; doublexx[199],yy[299],xyft[249]; for(i=0;i<100;i++) x[i]=A1*sin(2*PI*f1*i*data); for(j=0;j<150;j++)//计算y(t) y[j]=A2*exp(-p*j*data*j*data)*sin(2*PI*f2*j*data); correl(x,y,100,150,xy); correl(x,x,100,100,xx); correl(y,y,150,150,yy); correl(y,x,150,100,yx); FILE*fp,*fp1,*fp2,*fp3,*fp4,*fp5,*fp6;//导出数据 inti1,j1,k,l,m,n; fp1=fopen("D:\\sh\\data\\x.txt","w"); fp2=fopen("D:\\sh\\data\\y.txt","w"); fp3=fopen("D:\\sh\\data\\xy.txt","w"); fp4=fopen("D:\\sh\\data\\xx.txt","w"); fp5=fopen("D:\\sh\\data\\yx.txt","w"); fp6=fopen("D:\\sh\\data\\yy.txt","w"); for(i=0;i<249;i++) xyft[i]=yx[i]; fp=fopen("D:\\sh\\data\\xyt.txt","w"); for(i=0;i<249;i++) fprintf(fp,"%f\n",xyft[i]); fclose(fp); for(i1=0;i1<100;i1++) fprintf(fp1,"%f\n",x[i1]); fclose(fp1); for(j1=0;j1<150;j1++) fprintf(fp2,"%f\n",y[j1]); fclose(fp2); for(k=0;k<249;k++) fprintf(fp3,"%f\n",xy[k]); fclose(fp3); for(l=0;l<199;l++) fprintf(fp4,"%f\n",xx[l]); fclose(fp4); for(m=0;m<249;m++) fprintf(fp5,"%f\n",yx[m]); fclose(fp5); for(n=0;n<299;n++) fprintf(fp6,"%f\n",yy[n]); fclose(fp6);}4、地震子波及频谱分析#include<stdio.h>#include"conv.cpp"#include"fft.cpp"#include"uni.cpp"#defineV64//*********宏定义的参数，可以修改的点数#defineW256voidmain(){ inti; FILE*fp; doubleA1=1,f=35,a=100;//********************a代表阿尔法α doubledata=0.004;//*************************抽样间隔 doubleb[V]={0},bi[V]={0},Fb[V]={0},ks[200]={0}; doublesz[W]={0},s[200]={0},szi[W]={0},Fsz[W]={0}; doublex[W]={0},xi[W]={0},Fx[W]={0},n[200]={0}; doublexwp[V]; for(i=0;i<41;i++) b[i]=A1*exp(-a*i*data)*sin(2*PI*f*i*data);//*****产生ξ(t) ks[29]=0.2,ks[64]=0.4,ks[80]=0.15,ks[102]=0.5,ks[114]=0.35,ks[145]=0.1,ks[156]=0.2; fp=fopen("D:\\da\\b.txt","w"); for(i=0;i<V;i++) fprintf(fp,"%f\n",b[i]); fclose(fp); fp=fopen("D:\\da\\ks.txt","w"); for(i=0;i<200;i++) fprintf(fp,"%f\n",ks[i]); fclose(fp); con(b,ks,sz,41,200); fp=fopen("D:\\da\\s.txt","w"); for(i=0;i<W;i++) fprintf(fp,"%f\n",sz[i]); fclose(fp); for(i=0;i<200;i++)//*******************甩掉后40个样点 s[i]=sz[i]; uni(-0.4,0.4,200,n);//*****************产生随机数 fp=fopen("D:\\da\\n.txt","w"); for(i=0;i<200;i++) fprintf(fp,"%f\n",n[i]); fclose(fp); for(i=0;i<200;i++) x[i]=s[i]+n[i]; fp=fopen("D:\\da\\x.txt","w"); for(i=0;i<W;i++) fprintf(fp,"%f\n",x[i]); fclose(fp);fft(sz,szi,Fsz,W,1);//**********fft(b,bi,Fb,V,1);//************************用fft对b[t],x[t]频谱分析fft(x,xi,Fx,W,1);fp=fopen("D:\\da\\Fs.txt","w"); for(i=0;i<W;i++) fprintf(fp,"%f\n",Fsz[i]); fclose(fp);fp=fopen("D:\\da\\Fb.txt","w"); for(i=0;i<V;i++) fprintf(fp,"%f\n",Fb[i]); fclose(fp); fp=fopen("D:\\da\\Fx.txt","w"); for(i=0;i<W;i++) fprintf(fp,"%f\n",Fx[i]); fclose(fp); for(i=0;i<V;i++) xwp[i]=atan(bi[i]/b[i]); fp=fopen("D:\\da\\xw.txt","w"); for(i=0;i<V;i++) fprintf(fp,"%f\n",xwp[i]); fclose(fp);}5、信噪比#include<stdio.h>#include<math.h>#include"conv.cpp"#include"uni.cpp"#definePI3.14159265#defineV64#defineW256voidmain(){ inti; doubleA1=1,f=35,a=100; doubledata=0.004,Sr,tps=0,tpn=0; doubleb[V]={0},sz[240]={0},s[200]={0},ks[200]={0}; doublex[W]={0},n[200]={0}; for(i=0;i<41;i++) b[i]=A1*exp(-a*i*data)*sin(2*PI*f*i*data); ks[29]=0.2,ks[64]=0.4,ks[80]=0.15; ks[102]=0.5,ks[114]=0.35,ks[145]=0.1,ks[156]=0.2; con(b,ks,sz,41,200); for(i=0;i<200;i++)//甩掉后40个样点 s[i]=sz[i]; uni(-0.4,0.4,200,n);//随机数 for(i=0;i<200;i++) x[i]=s[i]+n[i]; for(i=0;i<200;i++)//信噪比计算 { tps+=pow(s[i],2); tpn+=pow(n[i],2); } Sr=tps/tpn; printf("SNR1=%f\n",Sr);}补0影响#include<stdio.h>#include"fft.cpp"#defineV64voidmain(){ inti; FILE*fp; doubleA1=1,f=35,a=100; doubledata=0.004; doubleb[V]={0},bi[V]={0},Fb[V]={0}; for(i=0;i<41;i++) b[i]=A1*exp(-a*i*data)*sin(2*PI*f*i*data); fp=fopen("D:\\daa\\b1.txt","w"); for(i=0;i<V;i++) fprintf(fp,"%f\n",b[i]); fclose(fp);fft(b,bi,Fb,V,1);//**************用fft对b[t]频谱分析fp=fopen("D:\\daa\\Fb1.txt","w"); for(i=0;i<V;i++) fprintf(fp,"%f\n",Fb[i]); fclose(fp);}《最后》所用到的子程序相关：voidcorrel(doublex[],doubley[],intN,intM,doublez[]){ inti,j,l,ts=0;//ts代表相关序列的第一个数的序列号 doublet=0; l=1-M; for(i=-M+1;i<N;i++) { for(j=0;j<N;j++) if(j-i>=0&&(j-i)<M) t+=x[j]*y[j-i]; z[i+M-1]=t; t=0; } printf("\n**********\n"); printf("ts=%d\n",l); printf("**********\n");}褶积：//传入数组1,2以及存储数组3；1,2的长度voidcon(doublea[],doubleb[],doublec[],intM,intL){ inti,j,N; N=M+L-1; for(i=0;i<N;i++) { doubletp=0.0; for(j=0;j<M;j++) { if((i-j)>=0&&(i-j)<L) tp+=a[j]*b[i-j]; } c[i]=tp,tp=0.0; }}DFT：//dft子函数，传递实部，虚部，存储sx数组，振幅数组。点数，正反#include<math.h>#definePI3.14159265voiddft(doublex[],doubley[],doublex1[],doubley1[],doubleFn[],intN,intsign){ inti,j; doubleQ,D,C,c,s; Q=2*PI/N; for(i=0;i<N;i++) { for(j=0;j<N;j++) { D=Q*i*j; c=cos(D); s=sin(D)*sign; x1[i]+=c*x[j]+s*y[j]; y1[i]+=c*y[j]-s*x[j]; } } if(sign==-1) { C=1.0/N; for(i=0;i<N;i++) { x1[i]=C*x1[i],y1[i]=C*y1[i]; }}for(i=0;i<N;i++)Fn[i]=sqrt(pow(x1[i],2)+pow(y1[i],2));}FFT://fft子函数，传递实部，虚部，存储振幅数组，点数，正反#include<math.h>#definePI3.14159265voidfft(doublex[],doubley[],doubleFn[],intN,intsign){ inti,j,r,m1,m2,m3,m4,k1,k2,l,k; doubleu