課程設計：信號分析與處理C語言編程

勘查技術課程設計：信號分析與處理基礎（西南石油大學---資源與環(huán)境學院）對于勘查技術與工程專業(yè)的學生來說，《信號分析與處理基礎》是一門專業(yè)基礎課，我是2010級的，我們是在大三第一學期上的，這門課數學與物理知識要求比較高，不過一開認真仔細學的話，也會學的好的，起碼要比那空洞、生奧、蛋疼《彈性波動力學》好學些。隨著課程的結束，《信號分析與處理基礎》的課程設計也隨之而來，我們是老師布置了4個題目，分單號與雙號各自做2道，我是單號，做的是濾波與相關。這次課程設計，注意考驗大家的編程能力，目前我們學過得就只有C語言，可以用Fortran,Matlab等等，Matlab可以現學現用，上手快。但是大家也可以挑戰(zhàn)下自己的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)的區(qū)別，還有開單雙門的區(qū)別與差異，時域與頻域濾波誰好、為什么？分析：從單雙門實部圖形看出，與有效波是完全一樣的，但是幅值變大，這體現了濾波突顯有效波的特性；從圖形看出虛部對結果沒有影響；單雙門虛部完全不同，這是由于開單雙門效果不同引起的。單門虛部變化大，而且幅值與起伏變化也大，而雙門幅值很小，小到可以忽略不計。按理說反變化IFFT后應該只有實部，沒有虛部，至于為什么會產生虛部，希望讀者自己下去研究下，希望大家相互交流、交流！第一個注意問題：單門與雙門誰好？答案當然是：雙門好。因為原始有效波x(t)是實數信號，對應的頻譜是偶對稱的，單門時，只是讓一部分通過了；而雙門則全部通過，肯定效果要好！而且實部波形，明顯看出雙門是干擾部分起伏比單門時要平緩，小得多了。第二個注意問題：時域與頻域哪個好？頻域要好，首先從兩者波形上可以看出，其次就是頻域濾波，只讓有效波通過，而之外的完全被濾掉，可謂是真正的理想狀態(tài)。第三個注意問題：頻域濾波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);//統(tǒng)一導出 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ǒng)(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