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

勘查技术课程设计：信号分析与处理基础

（西南石油大学---资源与环境学院）

对于勘查技术与工程专业的学生来说，《信号分析与处理基础》是一门专业基础课，我是2010级的，

我们是在大三第一学期上的，这门课数学与物理知识要求比较高，不过一开认真仔细学的话，也会学的好的，起码要比那空洞、生奥、蛋疼《弹性波动力学》好学些。

随着课程的结束，《信号分析与处理基础》的课程设计也随之而来，我们是老师布置了 4 个题目，分单号与双号各自做2道，我是单号，做的是滤波与相关。

这次课程设计，注意考验大家的编程能力，目前我们学过得就只有C 语言，可以用Fortran,Matlab 等等，Matlab 可以现学现用，上手快。但是大家也可以挑战下自己的C 语言，提高下自己的编程能力，这是一次很好地机会，真正实用的时刻。我就是用C 语言编的。

其余2题，我也把程序与结果图收集到了这里，以供学弟、学妹们参考之用！我的QQ ：593066480，有什么不懂的，或者好的见教，欢迎来信息交流！

题目如下：

1、滤波

已知原始地震记录x(t),要求：设计滤波器，消除x(t)中10Hz 以下，80Hz 以上的干扰信号。

)()()(t n t s t x +=

)2sin()2sin()(22112

2

t f e

A t f e

A t s t t ππββ--+=

)]2cos()2[sin()(433t f t f A t n ππ+=

建议参数： A 1=1，A 2=0.8，A 3=0.5

f 1=25Hz ，f 2=45Hz ，f 3=5Hz ，f 4=80 取样点数：N=200

抽样间隔：Δ=0.004，?=100

○

1时域滤波： 由)(*)()(t h t x t y = (h(t)为滤波因子)

)sin()cos(2

)(0t t t

t h ωωπ?=

建议参数：第一参数：ff1=20Hz ，ff2=50Hz 第二参数：ff1=25Hz ，ff2=45Hz 抽样点数：M=60 抽样间隔：Δ=0.004 112ff πω=，222ff πω= 2

1

2ωωω-=

?，2

1

20ωωω+=

单独计算：π

ω

?=

2)0(h

要求：画出x(t),h(t),y(t)图形，为了分析方便，也可以画出有效波s(t),干扰波n(t)及其频谱进行分析，如下图:

最后就是答辩，老师问问题，学生回答。主要注意几点就行了 1， 熟悉课本滤波部分知识；

2， 第二参数要比第一参数滤波效果好，因为门第一参数开大了，进来的干扰波也多了，从

第一参数：Fy 、第二参数Fy 图形上可以看出来，干扰波频谱被压小了。 第二参数压制了干扰波，突显了有效波，所有好。

○

2频域滤波 由公式：)()()(k H k X k Y ?=

1， 对x(t)进行补0，28

或者29

总之必须是2的次方，因为要用到FFT 公式与IFFT 公式，进行FFT 变换得到X （k ）；

2， H （k ）的求法：H （k ）也必须与X （k ）点数相同

N

f ?=

?1

， 单门：=??f m 120, =??f m 250

双门：=??f m 120, =??f m 250,3m =N-2m ,4m =N-1m ,

3, )()()(k H k X k Y ?=在用IFFT 反变换得y(t),存在实部与虚部，需要分析与处理 要求：画出H(k)、X(k)、Y(k)图形，并且分析X(k)、Y(k)的区别，还有开单双门的区别与差异，时域与频域滤波谁好、为什么？

3，分析：从单双门实部图形看出，与有效波是完全一样的，但是幅值变大，这体现了滤波突显有效波的特性；从图形看出虚部对结果没有影响；

单双门虚部完全不同，这是由于开单双门效果不同引起的。单门虚部变化大，而且幅值与起伏变化也大，而双门幅值很小，小到可以忽略不计。按理说反变化IFFT后应该只有实部，没有虚部，至于为什么会产生虚部，希望读者自己下去研究下，希望大家相互交流、交流！

第一个注意问题：单门与双门谁好？答案当然是：双门好。因为原始有效波x(t)是实数信号，对应的频谱是偶对称的，单门时，只是让一部分通过了；而双门则全部通过，肯定效果要好！而且实部波形，明显看出双门是干扰部分起伏比单门时要平缓，小得多了。

第二个注意问题：时域与频域哪个好？频域要好，首先从两者波形上可以看出，其次就是频域滤波，只让有效波通过，而之外的完全被滤掉，可谓是真正的理想状态。

第三个注意问题：频域滤波y(t)实部图200点以后，为什么有波形起伏？因为由于时域离散，必将导致频域周期化，这是由于频域周期化的结果。

2、相关

)

2sin()()2sin()(2211t f A t y t f A t x ππ==

建议参数：=1A 1，=1f 30Hz ，抽样点数：M=100 =2A 0.8，=2f 50Hz ，抽样点数：M=150 抽样间隔：Δ=0.004s

要求：画出)()(t y t x 、，)()()(τττ-Rxy Rxy Rxx 、、，)()(ττRyx Ryy 、的图形并分析验证书上自互相关性质的正确性！

3、快速褶积

)

()()()2sin()()2sin()(21222111t x t x t y t f A t x t f A t x *===ππ 建议参数：=1A 1，=1f 25Hz ，抽样点数：M=250 =2A 0.7，=2f 55Hz ，抽样点数：M=200 抽样间隔：Δ=0.004s

4、地震记录的生成和频谱分析

地震记录的生成和频谱分析、信噪比计算以及补0对频谱的影响，给定地震子波的数学表达式：)

2sin()(1

ft A t b t

πα-=

和反射系数序列：=30ξ0.2，=65ξ0.4，=81ξ0.15，=103ξ0.5 =115ξ0.35， =146ξ0.1， =157ξ0.2

产生一个含有随机干扰信号的地震信号：)()()()()()(t n t t b t n t s t x +*=+=ξ （其中n(t)可以用-0.4—+0.4之间的随机数代替） 要求：

1， 制作合成地震记录x(t),并对b(t)和s(t)做频谱分析（用FFT ），其中：b(t)

（N=41,Δ=0.004s ，f=35Hz ，α=100，A 1=1）反射序列和随机干扰N=200； 2， 计算x(t)信噪比，改变n(t)(用-0.8—+0.8之间的随机数代替)的大小再计

算。

信噪比：)()(2

2t n t s SNR ∑

∑= 3， 分别对b(t)后边、前边、中间补0，计算补0前后频谱的变化及补0多少对

频谱的影响。

附带程序：

最麻烦的就是编程序了，开始的时候，是很麻烦，不过只有去啃，就一定会有收获，比如：我开始编褶积程序的时候，整理了几天，上网查，翻图书，后来突然明白，靠公式就可以编出来了。只是FFT需要些功夫,其余都很快！书上介绍的是基2FFT，这个代码网上到处都有，想提升自己的编程能力，就去尝试编基4FFT，以及考虑基nFFT吧，祝大家好运。

1、时域滤波程序代码

#include

#include

#include"conv.cpp"

#include"dft.cpp"

#include"gyh.cpp"

#define T 200//************T表示x(t)点数，R表示n(t)点数

#define R 60

void 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;//******频率取值

double data=0.004;//******************取样点数

double s[T],n[T], x[T],h[R],y[T+R-1];

double si[T]={0},s1[T]={0},s1i[T]={0},Fs[T]={0};

double ni[T]={0},n1[T]={0},n1i[T]={0},Fn[T]={0};

double xi[T]={0},x1[T]={0},x1i[T]={0},Fx[T]={0};

double hi[T]={0},h1[T]={0},h1i[T]={0},Fh[T]={0};

double y1[T+R-1]={0},yb[T+R-1]={0},yb1[T+R-1]={0},Fy[T+R-1]={0};

for(i=0;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)进行频谱分析，用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:\\sh\\time1\\snxFsFnFx1.txt","w");

for(i=0;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)

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;

h[0]=2*dataw/PI;

for(j=1;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

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

fprintf(fp,"%f\t%f\n",y[i],Fy[i]);

fclose(fp);

printf("\nover\n\n");}

2、频域滤波程序代码：

#include

#include

#include"fft.cpp"

#include"ifft.cpp"

#define G 256

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 s[200],n[200];

double x[G]={0},x1[G]={0},h[G]={0},h1[G]={0},y[G]={0},y1[G]={0};

double Fx[G]={0},Fh1[G]={0},Fh2[G]={0},Fy1[G]={0},Fy2[G]={0};

double Y1[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

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

fprintf(fp,"%f\n",Fx[i]);

fclose(fp);

//处理h[t]

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

{

if(i>=int(m1)&&i<=int(m2))

Fh1[i]=1;

else Fh1[i]=0;

}

//计算,导出双门

for(i=0;i

{

if((i>=int(m1)&&i<=int(m2))||(i>=int(m3)&&i<=int(m4))) Fh2[i]=1;

else Fh2[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

{

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

{

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

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

fprintf(fp,"%f\t%f\t%f\t%f\n",Y1[i],Y1i[i],Y2[i],Y2i[i]);

fclose(fp);

printf("\nover\n\n");}

3、相关程序代码

#include

#include

#include"correl.cpp"

#define PI 3.14159265

void main()

{

int i,j;

double A1=1.0,A2=0.8;

double f1=30,f2=50;

double data=0.004,p=100;

double x[100]={0},y[150]={0},xy[249]={0},yx[249]={0};

double xx[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;//导出数据

int i1,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

#include"conv.cpp"

#include"fft.cpp"

#include"uni.cpp"

#define V 64//*********宏定义的参数，可以修改的点数

#define W 256

void main()

{

int i;

FILE *fp;

double A1=1,f=35,a=100;//********************a代表阿尔法α

double data=0.004;//************************* 抽样间隔

double b[V]={0},bi[V]={0},Fb[V]={0},ks[200]={0};

double sz[W]={0},s[200]={0},szi[W]={0},Fsz[W]={0};

double x[W]={0},xi[W]={0},Fx[W]={0},n[200]={0};

double xwp[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[1 56]=0.2;

fp=fopen("D:\\da\\b.txt","w");

for(i=0;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

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

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

fprintf(fp,"%f\n",Fsz[i]);

fclose(fp);

fp=fopen("D:\\da\\Fb.txt","w");

for(i=0;i

fprintf(fp,"%f\n",Fb[i]);

fclose(fp);

fp=fopen("D:\\da\\Fx.txt","w");

for(i=0;i

fprintf(fp,"%f\n",Fx[i]);

fclose(fp);

for(i=0;i

xwp[i]=atan(bi[i]/b[i]);

fp=fopen("D:\\da\\xw.txt","w");

for(i=0;i

fprintf(fp,"%f\n",xwp[i]);

fclose(fp);}

5、信噪比

#include

#include

#include"conv.cpp"

#include"uni.cpp"

#define PI 3.14159265

#define V 64

#define W 256

void main()

{

int i;

double A1=1,f=35,a=100;

double data=0.004,Sr,tps=0,tpn=0;

double b[V]={0},sz[240]={0},s[200]={0},ks[200]={0};

double x[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);}

6、补0影响