课程设计信号分析与处理C语言编程.docx
《课程设计信号分析与处理C语言编程.docx》由会员分享，可在线阅读，更多相关《课程设计信号分析与处理C语言编程.docx（20页珍藏版）》请在冰豆网上搜索。
课程设计信号分析与处理C语言编程
勘查技术课程设计：
信号分析与处理基础
（西南石油大学---资源与环境学院）
对于勘查技术与工程专业的学生来说，《信号分析与处理基础》是一门专业基础课，我是2010级的，我们是在大三第一学期上的，这门课数学与物理知识要求比较高，不过一开认真仔细学的话，也会学的好的，起码要比那空洞、生奥、蛋疼《弹性波动力学》好学些。
随着课程的结束，《信号分析与处理基础》的课程设计也随之而来，我们是老师布置了4个题目，分单号与双号各自做2道，我是单号，做的是滤波与相关。
这次课程设计，注意考验大家的编程能力，目前我们学过得就只有C语言，可以用Fortran,Matlab等等，Matlab可以现学现用，上手快。
但是大家也可以挑战下自己的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=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的次方，因为要用到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，分析：
从单双门实部图形看出，与有效波是完全一样的，但是幅值变大，这体现了滤波突显有效波的特性；从图形看出虚部对结果没有影响；
单双门虚部完全不同，这是由于开单双门效果不同引起的。
单门虚部变化大，而且幅值与起伏变化也大，而双门幅值很小，小到可以忽略不计。
按理说反变化IFFT后应该只有实部，没有虚部，至于为什么会产生虚部，希望读者自己下去研究下，希望大家相互交流、交流！
第一个注意问题：
单门与双门谁好？
答案当然是：
双门好。
因为原始有效波x（t）是实数信号，对应的频谱是偶对称的，单门时，只是让一部分通过了；而双门则全部通过，肯定效果要好！
而且实部波形，明显看出双门是干扰部分起伏比单门时要平缓，小得多了。
第二个注意问题：
时域与频域哪个好？
频域要好，首先从两者波形上可以看出，其次就是频域滤波，只让有效波通过，而之外的完全被滤掉，可谓是真正的理想状态。
第三个注意问题：
频域滤波y（t）实部图200点以后，为什么有波形起伏？
因为由于时域离散，必将导致频域周期化，这是由于频域周期化的结果。
2、相关
建议参数：
1，30Hz，抽样点数：
M=100
0.8，50Hz，抽样点数：
M=150
抽样间隔：
Δ=0.004s
要求：
画出，，的图形并分析验证书上自互相关性质的正确性！
3、快速褶积
建议参数：
1，25Hz，抽样点数：
M=250
0.7，55Hz，抽样点数：
M=200
抽样间隔：
Δ=0.004s
4、地震记录的生成和频谱分析
地震记录的生成和频谱分析、信噪比计算以及补0对频谱的影响，给定地震子波的数学表达式：
和反射系数序列：
0.2，0.4，0.15，0.5
0.35，0.1，0.2
产生一个含有随机干扰信号的地震信号：
（其中n（t）可以用-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需要些功夫,其余都很快！
书上介绍的是基2FFT，这个代码网上到处都有，想提升自己的编程能力，就去尝试编基4FFT，以及考虑基nFFT吧，祝大家好运。
1、时域滤波程序代码
#include
#include
#include"conv.cpp"
#include"dft.cpp"
#include"gyh.cpp"
#defineT200//************T表示x（t）点数，R表示n（t）点数
#defineR60
voidmain（）
{
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{
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;ifprintf（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{
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;jfprintf（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;ifprintf（fp,"%f\t%f\n",y[i],Fy[i]）;
fclose（fp）;
printf（"\nover\n\n"）;}
2、频域滤波程序代码：
#include
#include
#include"fft.cpp"
#include"ifft.cpp"
#defineG256
voidmain（）
{
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;ifprintf（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;ifprintf（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{
if（i>=int（m1）&&i<=int（m2））
Fh1[i]=1;
elseFh1[i]=0;
}
//计算,导出双门
for（i=0;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{
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;ifprintf（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;ifprintf（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"
#definePI3.14159265
voidmain（）
{
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
#include"conv.cpp"
#include"fft.cpp"
#include"uni.cpp"
#defineV64//*********宏定义的参数，可以修改的点数
#defineW256
voidmain（）
{
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;ifprintf（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;ifprintf（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;ifprintf（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;ifprintf（fp,"%f\n",Fsz[i]）;
fclose（fp）;
fp=fopen（"D:
\\da\\Fb.txt","w"）;
for（i=0;ifprintf（fp,"%f\n",Fb[i]）;
fclose（fp）;
fp=fopen（"D:
\\da\\Fx.txt","w"）;
for（i=0;ifprintf（fp,"%f\n",Fx[i]）;
fclose（fp）;
for（i=0;ixwp[i]=atan（bi[i]/b[i]）;
fp=fopen（"D:
\\da\\xw.txt","w"）;
for（i=0;ifprintf（fp,"%f\n",xwp[i]）;
fclose（fp）;}
5、信噪比
#include
#include
#include"conv.cpp"
#include"uni.cpp"
#definePI3.14159265
#defineV64
#defineW256
voidmain（）
{
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）;
升级会员 