課程設(shè)計(jì)：信號分析與處理C語言編程

勘查技術(shù)課程設(shè)計(jì)：信號分析與處理基礎(chǔ)（西南石油大學(xué)---資源與環(huán)境學(xué)院）對于勘查技術(shù)與工程專業(yè)的學(xué)生來說，《信號分析與處理基礎(chǔ)》是一門專業(yè)基礎(chǔ)課，我是2010級的，我們是在大三第一學(xué)期上的，這門課數(shù)學(xué)與物理知識要求比較高，不過一開認(rèn)真仔細(xì)學(xué)的話，也會學(xué)的好的，起碼要比那空洞、生奧、蛋疼《彈性波動力學(xué)》好學(xué)些。隨著課程的結(jié)束，《信號分析與處理基礎(chǔ)》的課程設(shè)計(jì)也隨之而來，我們是老師布置了4個題目，分單號與雙號各自做2道，我是單號，做的是濾波與相關(guān)。這次課程設(shè)計(jì)，注意考驗(yàn)大家的編程能力，目前我們學(xué)過得就只有C語言，可以用Fortran,Matlab等等，Matlab可以現(xiàn)學(xué)現(xiàn)用，上手快。但是大家也可以挑戰(zhàn)下自己的C語言，提高下自己的編程能力，這是一次很好地機(jī)會，真正實(shí)用的時刻。我就是用C語言編的。其余2題，我也把程序與結(jié)果圖收集到了這里，以供學(xué)弟、學(xué)妹們參考之用！我的QQ：593066480，有什么不懂的，或者好的見教，歡迎來信息交流！題目如下：濾波已知原始地震記錄x(t),要求：設(shè)計(jì)濾波器，消除x(t)中10Hz以下，80Hz以上的干擾信號。建議參數(shù)：A1=1，A2=0.8，A3=0.5f1=25Hz，f2=45Hz，f3=5Hz，f4=80取樣點(diǎn)數(shù)：N=200抽樣間隔：Δ=0.004，?=100eq\o\ac(○,1)時域?yàn)V波：由(h(t)為濾波因子)建議參數(shù)：第一參數(shù)：ff1=20Hz，ff2=50Hz第二參數(shù)：ff1=25Hz，ff2=45Hz抽樣點(diǎn)數(shù)：M=60抽樣間隔：Δ=0.004，，單獨(dú)計(jì)算：要求：畫出x(t),h(t),y(t)圖形，為了分析方便，也可以畫出有效波s(t),干擾波n(t)及其頻譜進(jìn)行分析，如下圖:最后就是答辯，老師問問題，學(xué)生回答。主要注意幾點(diǎn)就行了熟悉課本濾波部分知識；第二參數(shù)要比第一參數(shù)濾波效果好，因?yàn)殚T第一參數(shù)開大了，進(jìn)來的干擾波也多了，從第一參數(shù)：Fy、第二參數(shù)Fy圖形上可以看出來，干擾波頻譜被壓小了。第二參數(shù)壓制了干擾波，突顯了有效波，所有好。eq\o\ac(○,2)頻域?yàn)V波由公式：對x(t)進(jìn)行補(bǔ)0，28或者29總之必須是2的次方，因?yàn)橐玫紽FT公式與IFFT公式，進(jìn)行FFT變換得到X（k）；H（k）的求法：H（k）也必須與X（k）點(diǎn)數(shù)相同，單門：20,50雙門：20,50,=N-,=N-,3,在用IFFT反變換得y(t),存在實(shí)部與虛部，需要分析與處理要求：畫出H(k)、X(k)、Y(k)圖形，并且分析X(k)、Y(k)的區(qū)別，還有開單雙門的區(qū)別與差異，時域與頻域?yàn)V波誰好、為什么？分析：從單雙門實(shí)部圖形看出，與有效波是完全一樣的，但是幅值變大，這體現(xiàn)了濾波突顯有效波的特性；從圖形看出虛部對結(jié)果沒有影響；單雙門虛部完全不同，這是由于開單雙門效果不同引起的。單門虛部變化大，而且幅值與起伏變化也大，而雙門幅值很小，小到可以忽略不計(jì)。按理說反變化IFFT后應(yīng)該只有實(shí)部，沒有虛部，至于為什么會產(chǎn)生虛部，希望讀者自己下去研究下，希望大家相互交流、交流！第一個注意問題：單門與雙門誰好？答案當(dāng)然是：雙門好。因?yàn)樵加行Р▁(t)是實(shí)數(shù)信號，對應(yīng)的頻譜是偶對稱的，單門時，只是讓一部分通過了；而雙門則全部通過，肯定效果要好！而且實(shí)部波形，明顯看出雙門是干擾部分起伏比單門時要平緩，小得多了。第二個注意問題：時域與頻域哪個好？頻域要好，首先從兩者波形上可以看出，其次就是頻域?yàn)V波，只讓有效波通過，而之外的完全被濾掉，可謂是真正的理想狀態(tài)。第三個注意問題：頻域?yàn)V波y(t)實(shí)部圖200點(diǎn)以后，為什么有波形起伏？因?yàn)橛捎跁r域離散，必將導(dǎo)致頻域周期化，這是由于頻域周期化的結(jié)果。相關(guān)建議參數(shù)：1，30Hz，抽樣點(diǎn)數(shù)：M=1000.8，50Hz，抽樣點(diǎn)數(shù)：M=150抽樣間隔：Δ=0.004s要求：畫出，，的圖形并分析驗(yàn)證書上自互相關(guān)性質(zhì)的正確性！快速褶積建議參數(shù)：1，25Hz，抽樣點(diǎn)數(shù)：M=2500.7，55Hz，抽樣點(diǎn)數(shù)：M=200抽樣間隔：Δ=0.004s地震記錄的生成和頻譜分析地震記錄的生成和頻譜分析、信噪比計(jì)算以及補(bǔ)0對頻譜的影響，給定地震子波的數(shù)學(xué)表達(dá)式：和反射系數(shù)序列：0.2，0.4，0.15，0.50.35，0.1，0.2產(chǎn)生一個含有隨機(jī)干擾信號的地震信號：（其中n(t)可以用-0.4—+0.4之間的隨機(jī)數(shù)代替）要求：制作合成地震記錄x(t),并對b(t)和s(t)做頻譜分析（用FFT），其中：b(t)（N=41,Δ=0.004s，f=35Hz，α=100，A1=1）反射序列和隨機(jī)干擾N=200；計(jì)算x(t)信噪比，改變n(t)(用-0.8—+0.8之間的隨機(jī)數(shù)代替)的大小再計(jì)算。信噪比：分別對b(t)后邊、前邊、中間補(bǔ)0，計(jì)算補(bǔ)0前后頻譜的變化及補(bǔ)0多少對頻譜的影響。附帶程序：最麻煩的就是編程序了，開始的時候，是很麻煩，不過只有去啃，就一定會有收獲，比如：我開始編褶積程序的時候，整理了幾天，上網(wǎng)查，翻圖書，后來突然明白，靠公式就可以編出來了。只是FFT需要些功夫,其余都很快！書上介紹的是基2FFT，這個代碼網(wǎng)上到處都有，想提升自己的編程能力，就去嘗試編基4FFT，以及考慮基nFFT吧，祝大家好運(yùn)。1、時域?yàn)V波程序代碼#include<stdio.h>#include<math.h>#include"conv.cpp"#include"dft.cpp"#include"gyh.cpp"#defineT200//************T表示x(t)點(diǎn)數(shù)，R表示n(t)點(diǎn)數(shù)#defineR60voidmain(){ inti,j,B=100;//*********************B表示貝塔的值 FILE*fp; doubleA1=1.0,A2=0.8,A3=0.5;//********產(chǎn)生x(t) doublef1=25,f2=45,f3=5,f4=80;//******頻率取值 doubledata=0.004;//******************取樣點(diǎn)數(shù) 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)進(jìn)行頻譜分析，用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)一導(dǎo)出 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);//產(chǎn)生h(t) floatff1=20,ff2=50;//*****************可以修改的h(t)參數(shù) 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);//導(dǎo)出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");}頻域?yàn)V波程序代碼：#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;//產(chǎn)生x(t) doublef1=25,f2=45,f3=5,f4=80;//參數(shù) 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]; } //導(dǎo)出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]頻譜計(jì)算fft(x,x1,Fx,G,1);//導(dǎo)出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;//計(jì)算,導(dǎo)出單門for(i=0;i<G;i++){ if(i>=int(m1)&&i<=int(m2)) Fh1[i]=1; elseFh1[i]=0;} //計(jì)算,導(dǎo)出雙門 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");//導(dǎo)出單雙門for(i=0;i<G;i++) fprintf(fp,"%f\t%f\n",Fy1[i],Fy2[i]); fclose(fp);//Y[k]作反變換IFFT,并且導(dǎo)出實(shí)部，虛步 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");}相關(guā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++)//計(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;//導(dǎo)出數(shù)據(jù) 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//*********宏定義的參數(shù)，可以修改的點(diǎn)數(shù)#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);//*****產(chǎn)生ξ(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個樣點(diǎn) s[i]=sz[i]; uni(-0.4,0.4,200,n);//*****************產(chǎn)生隨機(jī)數(shù) 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個樣點(diǎn) s[i]=sz[i]; uni(-0.4,0.4,200,n);//隨機(jī)數(shù) for(i=0;i<200;i++) x[i]=s[i]+n[i]; for(i=0;i<200;i++)//信噪比計(jì)算 { tps+=pow(s[i],2); tpn+=pow(n[i],2); } Sr=tps/tpn; printf("SNR1=%f\n",Sr);}補(bǔ)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);}《最后》所用到的子程序相關(guān)：voidcorrel(doublex[],doubley[],intN,intM,doublez[]){ inti,j,l,ts=0;//ts代表相關(guān)序列的第一個數(shù)的序列號 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");}褶積：//傳入數(shù)組1,2以及存儲數(shù)組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子函數(shù)，傳遞實(shí)部，虛部，存儲sx數(shù)組，振幅數(shù)組。點(diǎn)數(shù)，正反#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子函數(shù)，傳遞實(shí)部，虛部，存儲振幅數(shù)組，點(diǎn)數(shù)，正反#include<math.h>#definePI3.14159265voidfft(doublex[],doubley[],doubleFn[],intN,intsign){ inti,j,r,m1,m2,m3,m4,k1,k2,l,k; doubleu