勘查技术课程设计:信号分析与处理基础(西南石油大学---资源与环境学院)对于勘查技术与工程专业的学生来说,《信号分析与处理基础》是一门专业基础课,我是2010级的,我们是在大三第一学期上的,这门课数学与物理知识要求比较高,不过一开认真仔细学的话,也会学的好的,起码要比那空洞、生奥、蛋疼《弹性波动力学》好学些。随着课程的结束,《信号分析与处理基础》的课程设计也随之而来,我们是老师布置了4个题目,分单号与双号各自做2道,我是单号,做的是滤波与相关。这次课程设计,注意考验大家的编程能力,目前我们学过得就只有C语言,可以用Fortran,Matlab等等,Matlab可以现学现用,上手快。但是大家也可以挑战下自己的C语言,提高下自己的编程能力,这是一次很好地机会,真正实用的时刻。我就是用C语言编的。其余2题,我也把程序与结果图收集到了这里,以供学弟、学妹们参考之用!题目如下:1、滤波已知原始地震记录x(t),要求:设计滤波器,消除x(t)中10Hz以下,80Hz以上的干扰信号。建议参数:A1=1,A2=,A3=f1=25Hz,f2=45Hz,f3=5Hz,f4=80取样点数:N=200抽样间隔:Δ=,?=100○1时域滤波:由)(*)()(thtxty(h(t)为滤波因子)建议参数:第一参数:ff1=20Hz,ff2=50Hz第二参数:ff1=25Hz,ff2=45Hz抽样点数:M=60抽样间隔:Δ=112ff,222ff212,2120单独计算:2)0(h要求:画出x(t),h(t),y(t)图形,为了分析方便,也可以画出有效波s(t),干扰波n(t)及其频谱进行分析,如下图:最后就是答辩,老师问问题,学生回答。主要注意几点就行了1,熟悉课本滤波部分知识;2,第二参数要比第一参数滤波效果好,因为门第一参数开大了,进来的干扰波也多了,从第一参数:Fy、第二参数Fy图形上可以看出来,干扰波频谱被压小了。第二参数压制了干扰波,突显了有效波,所有好。○2频域滤波由公式:)()()(kHkXkY1,对x(t)进行补0,28或者29总之必须是2的次方,因为要用到FFT公式与IFFT公式,进行FFT变换得到X(k);2,H(k)的求法:H(k)也必须与X(k)点数相同Nf1,单门:fm120,fm250双门:fm120,fm250,3m=N-2m,4m=N-1m,3,)()()(kHkXkY在用IFFT反变换得y(t),存在实部与虚部,需要分析与处理要求:画出H(k)、X(k)、Y(k)图形,并且分析X(k)、Y(k)的区别,还有开单双门的区别与差异,时域与频域滤波谁好、为什么?3,分析:从单双门实部图形看出,与有效波是完全一样的,但是幅值变大,这体现了滤波突显有效波的特性;从图形看出虚部对结果没有影响;单双门虚部完全不同,这是由于开单双门效果不同引起的。单门虚部变化大,而且幅值与起伏变化也大,而双门幅值很小,小到可以忽略不计。按理说反变化IFFT后应该只有实部,没有虚部,至于为什么会产生虚部,希望读者自己下去研究下,希望大家相互交流、交流!第一个注意问题:单门与双门谁好?答案当然是:双门好。因为原始有效波x(t)是实数信号,对应的频谱是偶对称的,单门时,只是让一部分通过了;而双门则全部通过,肯定效果要好!而且实部波形,明显看出双门是干扰部分起伏比单门时要平缓,小得多了。第二个注意问题:时域与频域哪个好?频域要好,首先从两者波形上可以看出,其次就是频域滤波,只让有效波通过,而之外的完全被滤掉,可谓是真正的理想状态。第三个注意问题:频域滤波y(t)实部图200点以后,为什么有波形起伏?因为由于时域离散,必将导致频域周期化,这是由于频域周期化的结果。2、相关建议参数:1A1,1f30Hz,抽样点数:M=1002A,2f50Hz,抽样点数:M=150抽样间隔:Δ=要求:画出)()(tytx、,)()()(RxyRxyRxx、、,)()(RyxRyy、的图形并分析验证书上自互相关性质的正确性!3、快速褶积建议参数:1A1,1f25Hz,抽样点数:M=2502A,2f55Hz,抽样点数:M=200抽样间隔:Δ=4、地震记录的生成和频谱分析地震记录的生成和频谱分析、信噪比计算以及补0对频谱的影响,给定地震子波的数学表达式:)2sin()(1ftAtbt和反射系数序列:30,65,81,103115,146,157产生一个含有随机干扰信号的地震信号:)()()()()()(tnttbtntstx(其中n(t)可以用—+之间的随机数代替)要求:1,制作合成地震记录x(t),并对b(t)和s(t)做频谱分析(用FFT),其中:b(t)(N=41,Δ=,f=35Hz,α=100,A1=1)反射序列和随机干扰N=200;2,计算x(t)信噪比,改变n(t)(用—+之间的随机数代替)的大小再计算。信噪比:)()(22tntsSNR3,分别对b(t)后边、前边、中间补0,计算补0前后频谱的变化及补0多少对频谱的影响。附带程序:最麻烦的就是编程序了,开始的时候,是很麻烦,不过只有去啃,就一定会有收获,比如:我开始编褶积程序的时候,整理了几天,上网查,翻图书,后来突然明白,靠公式就可以编出来了。只是FFT需要些功夫,其余都很快!书上介绍的是基2FFT,这个代码网上到处都有,想提升自己的编程能力,就去尝试编基4FFT,以及考虑基nFFT吧,祝大家好运。1、时域滤波程序代码#include#include#include#include#include#defineT200//************T表示x(t)点数,R表示n(t)点数#defineR60voidmain(){inti,j,B=100;//*********************B表示贝塔的值FILE*fp;doubleA1=,A2=,A3=;//********产生x(t)doublef1=25,f2=45,f3=5,f4=80;//******频率取值doubledata=;//******************取样点数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;iT;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\\,w);for(i=0;iT;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;jR;j++){h[j]=(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\\,w);for(j=0;jR;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\\,w);for(i=0;iT+R-1;i++)fprintf(fp,%f\t%f\n,y[i],Fy[i]);fclose(fp);printf(\nover\n\n);}2、频域滤波程序代码:#include#include#include#include#defineG256voidmain(){inti,B=100;FILE*fp;//定义指针文件doubleA1=,A2=,A3=;//产生x(t)doublef1=25,f2=45,f3=5,f4=80;//参数doubledata=;//抽样间隔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;i200;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\\,w);for(i=0;iG;i++)fprintf(fp,%f\n,x[i]);fclose(fp);//x[t]频谱计算fft(x,x1,Fx,G,1);//导出X(k)fp=fopen(D:\\sh\\py\\,w);for(i=0;iG;i++)fprintf(fp,%f\n,Fx[i]);fclose(fp);//处理h[t]doubledataf;doublem1,m2,m3,m4;dataf=(data*G);m1=dataf,m2=dataf;m3=G-m2,m4=G-m1;//计算,导出单门for(i=0;iG;i++){if(i=int(m1)&&i=int(m2))Fh1[i]=1;elseFh1[i]=0;}//计算,导出双门for(i=0;iG;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\\,w);for(i=0;i256;i++)fprintf(fp,%f\t%f\n,Fh1[i],Fh2[i]);fclose(fp);for(i=0;iG;i++)//单门滤波{Y1[i]=x[i]*Fh1[i];Y1i[i]=x1[i]*Fh1[i];Fy1[i]=sqrt(pow(Y1[i],2