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

1、课程设计信号分析与处理C语言编程勘查技术课程设计：信号分析与处理基础（西南石油大学-资源与环境学院） 对于勘查技术与工程专业的学生来说，信号分析与处理基础是一门专业基础课，我是2010级的，我们是在大三第一学期上的，这门课数学与物理知识要求比较高，不过一开认真仔细学的话，也会学的好的，起码要比那空洞、生奥、蛋疼弹性波动力学好学些。随着课程的结束，信号分析与处理基础的课程设计也随之而来，我们是老师布置了4 个题目，分单号与双号各自做2道，我是单号，做的是滤波与相关。这次课程设计，注意考验大家的编程能力，目前我们学过得就只有C语言，可以用Fortran,Matlab等等，Matlab可以现学现用，

2、上手快。但是大家也可以挑战下自己的C语言，提高下自己的编程能力，这是一次很好地机会，真正实用的时刻。我就是用C语言编的。其余2题，我也把程序与结果图收集到了这里，以供学弟、学妹们参考之用！我的QQ：593066480，有什么不懂的，或者好的见教，欢迎来信息交流！题目如下：1、滤波已知原始地震记录x(t),要求：设计滤波器，消除x(t)中10Hz以下，80Hz以上的干扰信号。建议参数： A1=1，A2=0.8，A3=0.5 f1=25Hz，f2=45Hz，f3=5Hz，f4=80 取样点数：N=200抽样间隔：=0.004，=100时域滤波：由 (h(t)为滤波因子)建议参数：第一参数：ff1=

3、20Hz，ff2=50Hz 第二参数：ff1=25Hz，ff2=45Hz 抽样点数：M=60 抽样间隔：=0.004 ， ，单独计算： 要求：画出x(t),h(t),y(t)图形，为了分析方便，也可以画出有效波s(t),干扰波n(t)及其频谱进行分析，如下图: 最后就是答辩，老师问问题，学生回答。主要注意几点就行了1，熟悉课本滤波部分知识；2，第二参数要比第一参数滤波效果好，因为门第一参数开大了，进来的干扰波也多了，从第一参数：Fy、第二参数Fy图形上可以看出来，干扰波频谱被压小了。第二参数压制了干扰波，突显了有效波，所有好。频域滤波由公式：1，对x(t)进行补0，28或者29总之必须是2的次

4、方，因为要用到FFT公式与IFFT公式，进行FFT变换得到X（k）；2，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)的区别，还有开单双门的区别与差异，时域与频域滤波谁好、为什么？3，分析：从单双门实部图形看出，与有效波是完全一样的，但是幅值变大，这体现了滤波突显有效波的特性；从图形看出虚部对结果没有影响；单双门虚部完全不同，这是由于开单双门效果不同引起的。单门虚部变化大，而且幅值与起伏变化也大，而双门

5、幅值很小，小到可以忽略不计。按理说反变化IFFT后应该只有实部，没有虚部，至于为什么会产生虚部，希望读者自己下去研究下，希望大家相互交流、交流！第一个注意问题：单门与双门谁好？答案当然是：双门好。因为原始有效波x(t)是实数信号，对应的频谱是偶对称的，单门时，只是让一部分通过了；而双门则全部通过，肯定效果要好！而且实部波形，明显看出双门是干扰部分起伏比单门时要平缓，小得多了。第二个注意问题：时域与频域哪个好？频域要好，首先从两者波形上可以看出，其次就是频域滤波，只让有效波通过，而之外的完全被滤掉，可谓是真正的理想状态。第三个注意问题：频域滤波y(t)实部图200点以后，为什么有波形起伏？因为由

6、于时域离散，必将导致频域周期化，这是由于频域周期化的结果。2、相关建议参数：1，30Hz，抽样点数：M=100 0.8，50Hz，抽样点数：M=150 抽样间隔：=0.004s要求：画出，的图形并分析验证书上自互相关性质的正确性！3、快速褶积建议参数：1，25Hz，抽样点数：M=250 0.7，55Hz，抽样点数：M=200 抽样间隔：=0.004s4、地震记录的生成和频谱分析地震记录的生成和频谱分析、信噪比计算以及补0对频谱的影响，给定地震子波的数学表达式：和反射系数序列：0.2，0.4，0.15，0.5 0.35， 0.1， 0.2产生一个含有随机干扰信号的地震信号：（其中n(t)可以用-

7、0.4+0.4之间的随机数代替）要求：1，制作合成地震记录x(t),并对b(t)和s(t)做频谱分析（用FFT），其中：b(t)（N=41,=0.004s，f=35Hz，=100，A1=1）反射序列和随机干扰N=200；2，计算x(t)信噪比，改变n(t)(用-0.8+0.8之间的随机数代替)的大小再计算。信噪比：3，分别对b(t)后边、前边、中间补0，计算补0前后频谱的变化及补0多少对频谱的影响。附带程序：最麻烦的就是编程序了，开始的时候，是很麻烦，不过只有去啃，就一定会有收获，比如：我开始编褶积程序的时候，整理了几天，上网查，翻图书，后来突然明白，靠公式就可以编出来了。只是FFT需要些功夫

8、,其余都很快！书上介绍的是基2FFT，这个代码网上到处都有，想提升自己的编程能力，就去尝试编基4FFT，以及考虑基nFFT吧，祝大家好运。1、时域滤波程序代码#include#include#includeconv.cpp#includedft.cpp#includegyh.cpp#define T 200/*T表示x(t)点数，R表示n(t)点数#define R 60void main() int i,j,B=100;/*B表示贝塔的值 FILE *fp; double A1=1.0,A2=0.8,A3=0.5;/*产生x(t) double f1=25,f2=45,f3=5,f4=80;

9、/*频率取值 double data=0.004;/*取样点数 double sT,nT, xT,hR,yT+R-1; double siT=0,s1T=0,s1iT=0,FsT=0; double niT=0,n1T=0,n1iT=0,FnT=0; double xiT=0,x1T=0,x1iT=0,FxT=0; double hiT=0,h1T=0,h1iT=0,FhT=0; double y1T+R-1=0,ybT+R-1=0,yb1T+R-1=0,FyT+R-1=0; for(i=0;iT;i+) si=A1*exp(-B*pow(i*data),2) * sin(2*PI*f1*i*

10、data) + A2*exp(-B*pow(i*data),2) * sin(2*PI*f2*i*data); ni=A3*(sin(2*PI*f3*i*data)+cos(2*PI*f4*i*data); xi=si+ni; /对s(t)进行频谱分析，用DFT dft(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:shtime1snxFsFnFx1.txt,w);for(i=0;iT;i+) fprintf(fp,%ft%ft%ft%ft%ft%fn,si,ni,x

11、i,Fsi,Fni,Fxi);fclose(fp);/产生h(t) float ff1=20,ff2=50;/*可以修改的h(t)参数 double w1,w2,dataw,w0; w1=2*PI*ff1,w2=2*PI*ff2; dataw=(w2-w1)/2,w0=(w2+w1)/2; h0=2*dataw/PI; for(j=1;jR;j+) hj=(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:shtime1hFh.txt,w); for(

12、j=0;jR;j+) fprintf(fp,%ft%fn,hj,Fhj); 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:shtime1yFy1.txt,w);for(i=0;iT+R-1;i+) fprintf(fp,%ft%fn,yi,Fyi);fclose(fp); printf(novernn);2、频域滤波程序代码：#include#include#includefft.cpp#includeifft.cpp#define G 256

13、void main() int i,B=100; FILE *fp;/定义指针文件 double A1=1.0,A2=0.8,A3=0.5;/产生x(t) double f1=25,f2=45,f3=5,f4=80;/参数 double data=0.004;/抽样间隔 double s200,n200; double xG=0,x1G=0,hG=0,h1G=0,yG=0,y1G=0; double FxG=0,Fh1G=0,Fh2G=0,Fy1G=0,Fy2G=0; double Y1G=0,Y1iG=0,Y2G=0,Y2iG=0; for(i=0;i200;i+) si=A1*exp(-B

14、*pow(i*data),2) * sin(2*PI*f1*i*data) + A2*exp(-B*pow(i*data),2) * sin(2*PI*f2*i*data); ni=A3*(sin(2*PI*f3*i*data)+cos(2*PI*f4*i*data); xi=si+ni; /导出xt fp=fopen(D:shpyx.txt,w); for(i=0;iG;i+) fprintf(fp,%fn,xi); fclose(fp); /xt频谱计算fft(x,x1,Fx,G,1); /导出X(k) fp=fopen(D:shpyFx.txt,w); for(i=0;iG;i+) fp

15、rintf(fp,%fn,Fxi); fclose(fp); /处理ht double dataf; double m1,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=int(m1)&i=int(m2) Fh1i=1; else Fh1i=0; /计算,导出双门 for(i=0;i=int(m1)&i=int(m3)&i=int(m4) Fh2i=1; else Fh2i=0; fp=fopen(D:shpyFh12.txt,w); for(i=0;i

16、256;i+) fprintf(fp,%ft%fn,Fh1i,Fh2i); fclose(fp); for(i=0;iG;i+)/单门滤波 Y1i=xi*Fh1i; Y1ii=x1i*Fh1i; Fy1i=sqrt(pow(Y1i,2)+pow(Y1ii,2); for(i=0;iG;i+)/开双门 Y2i=xi*Fh2i; Y2ii=x1i*Fh2i; Fy2i=sqrt(pow(Y2i,2)+pow(Y2ii,2); fp=fopen(D:shpyFy12.txt,w);/导出单双门 for(i=0;iG;i+) fprintf(fp,%ft%fn,Fy1i,Fy2i); fclose(f

17、p);/Yk作反变换IFFT,并且导出实部，虚步 ifft(Y1,Y1i,G,-1); ifft(Y2,Y2i,G,-1); fp=fopen(D:shpyyt12.txt,w); for(i=0;iG;i+) fprintf(fp,%ft%ft%ft%fn,Y1i,Y1ii,Y2i,Y2ii); fclose(fp);printf(novernn);3、相关程序代码#include#include#includecorrel.cpp#define PI 3.14159265void main() int i,j; double A1=1.0,A2=0.8; double f1=30,f2=5

18、0; double data=0.004,p=100; double x100=0,y150=0,xy249=0,yx249=0; double xx199,yy299,xyft249; for(i=0;i100;i+) xi=A1*sin(2*PI*f1*i*data); for(j=0;j150;j+)/计算y(t) yj=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,

19、100,yx); FILE *fp,*fp1,*fp2,*fp3,*fp4,*fp5,*fp6;/导出数据 int i1,j1,k,l,m,n; fp1=fopen(D:shdatax.txt,w); fp2=fopen(D:shdatay.txt,w); fp3=fopen(D:shdataxy.txt,w); fp4=fopen(D:shdataxx.txt,w); fp5=fopen(D:shdatayx.txt,w); fp6=fopen(D:shdatayy.txt,w); for(i=0;i249;i+) xyfti=yxi; fp=fopen(D:shdataxyt.txt,w)

20、; for(i=0;i249;i+) fprintf(fp,%fn,xyfti); fclose(fp); for(i1=0;i1100;i1+) fprintf(fp1,%fn,xi1); fclose(fp1); for(j1=0;j1150;j1+) fprintf(fp2,%fn,yj1); fclose(fp2); for(k=0;k249;k+) fprintf(fp3,%fn,xyk); fclose(fp3); for(l=0;l199;l+) fprintf(fp4,%fn,xxl); fclose(fp4); for(m=0;m249;m+) fprintf(fp5,%fn

21、,yxm); fclose(fp5); for(n=0;n299;n+) fprintf(fp6,%fn,yyn); fclose(fp6);4、地震子波及频谱分析#include#includeconv.cpp#includefft.cpp#includeuni.cpp#define V 64/*宏定义的参数，可以修改的点数#define W 256void main() int i; FILE *fp; double A1=1,f=35,a=100;/*a代表阿尔法 double data=0.004;/* 抽样间隔 double bV=0,biV=0,FbV=0,ks200=0; dou

22、ble szW=0,s200=0,sziW=0,FszW=0; double xW=0,xiW=0,FxW=0,n200=0; double xwpV; for(i=0;i41;i+) bi=A1*exp(-a*i*data)*sin(2*PI*f*i*data);/*产生(t) ks29=0.2,ks64=0.4,ks80=0.15,ks102=0.5,ks114=0.35,ks145=0.1,ks156=0.2; fp=fopen(D:dab.txt,w); for(i=0;iV;i+) fprintf(fp,%fn,bi); fclose(fp); fp=fopen(D:daks.txt

23、,w); for(i=0;i200;i+) fprintf(fp,%fn,ksi); fclose(fp); con(b,ks,sz,41,200); fp=fopen(D:das.txt,w); for(i=0;iW;i+) fprintf(fp,%fn,szi); fclose(fp); for(i=0;i200;i+)/*甩掉后40个样点 si=szi; uni(-0.4,0.4,200,n);/*产生随机数 fp=fopen(D:dan.txt,w); for(i=0;i200;i+) fprintf(fp,%fn,ni); fclose(fp); for(i=0;i200;i+) x

24、i=si+ni; fp=fopen(D:dax.txt,w); for(i=0;iW;i+) fprintf(fp,%fn,xi); fclose(fp);fft(sz,szi,Fsz,W,1);/*fft(b,bi,Fb,V,1);/*用fft对bt,xt频谱分析fft(x,xi,Fx,W,1); fp=fopen(D:daFs.txt,w); for(i=0;iW;i+) fprintf(fp,%fn,Fszi); fclose(fp); fp=fopen(D:daFb.txt,w); for(i=0;iV;i+) fprintf(fp,%fn,Fbi); fclose(fp); fp=f

25、open(D:daFx.txt,w); for(i=0;iW;i+) fprintf(fp,%fn,Fxi); fclose(fp); for(i=0;iV;i+) xwpi=atan(bii/bi); fp=fopen(D:daxw.txt,w); for(i=0;iV;i+) fprintf(fp,%fn,xwpi); fclose(fp);5、信噪比#include#include#includeconv.cpp#includeuni.cpp#define PI 3.14159265#define V 64#define W 256void main() int i; double A1

26、=1,f=35,a=100; double data=0.004,Sr,tps=0,tpn=0; double bV=0,sz240=0,s200=0,ks200=0; double xW=0,n200=0; for(i=0;i41;i+) bi=A1*exp(-a*i*data)*sin(2*PI*f*i*data); ks29=0.2,ks64=0.4,ks80=0.15; ks102=0.5,ks114=0.35,ks145=0.1,ks156=0.2; con(b,ks,sz,41,200); for(i=0;i200;i+)/甩掉后40个样点 si=szi; uni(-0.4,0.4,200,n);/随机数 for(i=0;i200;i+) xi=si+ni; for(i=0;i200;i+)/信噪比计算 tps+=pow(si,2); tpn+=pow(ni,2);