信号分析实用程序.docx
《信号分析实用程序.docx》由会员分享，可在线阅读，更多相关《信号分析实用程序.docx（27页珍藏版）》请在冰豆网上搜索。
信号分析实用程序
《信号分析》实用程序
C&FORTRAN
来自张军华老师的博客：
注：
大部分程序是用C语言编写的，少部分是用FORTRAN编写的，可以考虑用C++或者Java改写。
实际地震资料处理中语言不是重点。
对于地震信号处理方向感兴趣的同学可以关注一下张老师。
另外对于测井信号处理方向感兴趣的也有一定用处，测井实验用Matlab的居多，数学专业对于物探或测井都是挺好的发展方向。
张军华，男，党员，1965年4月生，博士，教授，浙江绍兴人，中国石油学会、SEG会员，校、院教代会代表，省和教育部科技评审专家。
1987年物探专业本科毕业留校任教，现在地球物理系从事石油物探的教学与研究工作。
主要研究成果和方向：
1）并行计算和地震资料并行处理方法；2）薄层的地球物理特征研究；3）相干体、复合属性的提取与应用；4）自适应拉东变换和多次波去噪；5）静校正方法研究；6）小波变换和多尺度处理、解释方法等。
目录
目录1
1．FFT子函数（fft_sub.c）1
2．DFT子函数（dft_sub.c）4
3．频谱泄露现象及窗函数的作用——实验五之二5
1、程序5
4．窗函数的描述9
1）布莱克曼窗9
5．频率域数字滤波程序——实验六11
1、程序11
2、参数文件13
6．频谱分析程序14
7．重采样程序15
8．自相关与互相关程序16
1）自相关程序16
2）互相关程序17
9．随机数的产生程序18
10．希尔伯特变换的物理意义18
11．Fourier级数、离散谱与振幅谱、相位谱的概念22
12．DFT和FFT子程序及二者速度比较——实验五之一24
13．多频信号的混叠与分频解释的体会26
14．卷积的意义和合成地震记录的制作27
15．信号的描述29
1、雷克子波的生成29
2、指数衰减函数的描述29
3、傅立叶核函数（Sinc函数）30
1．FFT子函数（fft_sub.c）
voidfft（floatsr[],floatsx[],intm0,intinv）
{ inti,j,lm,li,k,lmx,np,lix,mm2;
doublet1,t2,c,s,cv,st,ct;
if（m0<0）
return;
lmx=1;
for（i=1;i<=m0;++i）
lmx+=lmx; //form2**m0
cv=2.0*PI/（double）lmx;
ct=cos（cv）; st=-inv*sin（cv）;
np=lmx;mm2=m0-2;
/*fftbutterflynumeration*/
for（i=1;i<=mm2;++i）
{
lix=lmx;lmx/=2;
c=ct;s=st;
for（li=0;li {
j=li;k=j+lmx;
t1=sr[j]-sr[k];t2=sx[j]-sx[k];
sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=t1;sx[k]=t2;
++j;++k;
t1=sr[j]-sr[k];t2=sx[j]-sx[k];
sr[j]+=sr[k];sx[j]+=sx[k];
sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;
}
for（lm=2;lm {
cv=c;c=ct*c-st*s;s=st*cv+ct*s;
for（li=0;li {
j=li+lm;k=lmx+j;
t1=sr[j]-sr[k];t2=sx[j]-sx[k];
sr[j]+=sr[k];sx[j]+=sx[k];
sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;
}
}
cv=ct;ct=2.0*ct*ct-1.0;st=2.0*st*cv;
}
/*4pointsDFT*/
if（m0>=2）
for（li=0;li {
j=li;k=j+2;
t1=sr[j]-sr[k];t2=sx[j]-sx[k];
sr[j]+=sr[k];
sx[j]+=sx[k];
sr[k]=t1;sx[k]=t2;
++j;++k;
t1=sr[j]-sr[k];t2=sx[j]-sx[k];
sr[j]+=sr[k];sx[j]+=sx[k];
sr[k]=inv*t2;sx[k]=-inv*t1;
}
/*2pointsDFT*/
for（li=0;li {
j=li;k=j+1;
t1=sr[j]-sr[k];t2=sx[j]-sx[k];
sr[j]+=sr[k];sx[j]+=sx[k];
sr[k]=t1;sx[k]=t2;
}
/*sortaccordingtobitreversal*/
lmx=np/2;j=0;
for（i=1;i {
k=lmx;
while（k<=j）
{
j-=k;k/=2;
}
j+=k;
if（i {
t1=sr[j];sr[j]=sr[i];sr[i]=t1;
t1=sx[j];sx[j]=sx[i];sx[i]=t1;
}
}
/* ifInverseFFT,multiply1.0/np */
if（inv!
=-1）
return;
t1=1.0/np;
for（i=0;i {
sr[i]*=t1;sx[i]*=t1;
}}
2．DFT子函数（dft_sub.c）
voiddft（xr,xi,flag）
floatxr[N],xi[N];
intflag;
{floatXR[N],XI[N];
intk,n;
floatsum1,sum2,cita;
for（k=0;k<=N-1;k++）
{sum1=0.0;
sum2=0.0;
{for（n=0;n<=N-1;n++）
{ cita=2.0*3.1415926/N*n*k;
sum1=sum1+xr[n]*cos（cita）+flag*xi[n]*sin（cita）;
sum2=sum2-flag*xr[n]*sin（cita）+xi[n]*cos（cita）;
}
XR[k]=sum1;
XI[k]=sum2;
}
}
if（flag==1）
for（k=0;k<=N-1;k++）
{
xr[k]=XR[k];
xi[k]=XI[k];
}
else
for（k=0;k<=N-1;k++）
{
xr[k]=XR[k]/N;
xi[k]=XI[k]/N;
}
3．频谱泄露现象及窗函数的作用——实验五之二
1、程序
#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#definePI3.1415926
/* sr,si：
双精度型一维数组,输入（输出）信号的实部和虚部*/
/* m0:
2的次方数,2**m0=nfft*/
/* inv=1forwardtransform;inv=-1inversetransform */
voidfft（floatsr[],floatsi[],intm0,intinv）
{
inti,j,lm,li,k,lmx,np,lix,mm2;
doublet1,t2,c,s,cv,st,ct;
if（m0<0）
return;
lmx=1;
for（i=1;i<=m0;++i）
lmx+=lmx; //form2**m0
cv=2.0*PI/（double）lmx;
ct=cos（cv）; st=-inv*sin（cv）;
np=lmx;mm2=m0-2;
/*fftbutterflynumeration*/
for（i=1;i<=mm2;++i）
{
lix=lmx;lmx/=2;
c=ct;s=st;
for（li=0;li {
j=li;k=j+lmx;
t1=sr[j]-sr[k];t2=si[j]-si[k];
sr[j]+=sr[k];si[j]+=si[k];sr[k]=t1;si[k]=t2;
++j;++k;
t1=sr[j]-sr[k];t2=si[j]-si[k];
sr[j]+=sr[k];si[j]+=si[k];
sr[k]=c*t1-s*t2;si[k]=s*t1+c*t2;
}
for（lm=2;lm {
cv=c;c=ct*c-st*s;s=st*cv+ct*s;
for（li=0;li {
j=li+lm;k=lmx+j;
t1=sr[j]-sr[k];t2=si[j]-si[k];
sr[j]+=sr[k];si[j]+=si[k];
sr[k]=c*t1-s*t2;si[k]=s*t1+c*t2;
}
}
cv=ct;ct=2.0*ct*ct-1.0;st=2.0*st*cv;
}
/*4pointsDFT*/
if（m0>=2）
for（li=0;li {
j=li;k=j+2;
t1=sr[j]-sr[k];t2=si[j]-si[k];
sr[j]+=sr[k];
si[j]+=si[k];
sr[k]=t1;si[k]=t2;
++j;++k;
t1=sr[j]-sr[k];t2=si[j]-si[k];
sr[j]+=sr[k];si[j]+=si[k];
sr[k]=inv*t2;si[k]=-inv*t1;
}
/*2pointsDFT*/
for（li=0;li {
j=li;k=j+1;
t1=sr[j]-sr[k];t2=si[j]-si[k];
sr[j]+=sr[k];si[j]+=si[k];
sr[k]=t1;si[k]=t2;
}
/*sortaccordingtobitreversal*/
lmx=np/2;j=0;
for（i=1;i {
k=lmx;
while（k<=j）
{
j-=k;k/=2;
}
j+=k;
if（i {
t1=sr[j];sr[j]=sr[i];sr[i]=t1;
t1=si[j];si[j]=si[i];si[i]=t1;
}
}
/* ifInverseFFT,multiply1.0/np */
if（inv!
=-1）
return;
t1=1.0/np;
for（i=0;i {
sr[i]*=t1;si[i]*=t1;
}
}
voidmain（）
{voidfft（）;
float*xr,*xi;
float*w;
inti,np,nfft,k,flag;
floatt,dt,df,f,f1,f2,f3;
FILE*fpar,*fp1,*fp2;
charfil1[80],fil2[80];
fpar=fopen（"filter_window_par.txt","r"）;
fscanf（fpar,"%s",fil1）;
fscanf（fpar,"%s",fil2）;
printf（"%s\n",fil1）;
printf（"%s\n",fil2）;
fscanf（fpar,"%d",&np）;
printf（"np=%d\n",np）;
fscanf（fpar,"%f",&dt）;
printf（"dt=%8.3fms\n",dt）;
fscanf（fpar,"%f%f%f",&f1,&f2,&f3）;
printf（"f1=%8.3ff2=%8.3ff3=%8.3f\n",f1,f2,f3）;
fscanf（fpar,"%d",&flag）;
if（flag==1）printf（"Blackmanwindowedfunction\n"）;
elseprintf（"Rectangularwindowedfunction\n"）;
dt=dt/1000.0;
//calculatefftpoint
k=log（np）/log
（2）;
if（np>pow（2,k））k=k+1;
nfft=pow（2,k）;
df=1.0/（nfft*dt）;
printf（"nfft=%d k=%d\n",nfft,k）;
printf（"dt=%8.4f df=%8.4f\n",dt,df）;
//allocatememory
xr=（float*）calloc（nfft,sizeof（float））;
xi=（float*）calloc（nfft,sizeof（float））;
w=（float*）malloc（np*sizeof（float））;
//ifflag=1generateaBlackmanwindowedfunction
//elsegeneratearectangularwindowedfunction
if（flag==1）
for（i=0;i {t=2.0*PI*i/（np-1）;
w[i]=0.42-0.5*cos（t）+0.08*cos（2*t）; //BlackmanWindowFunction
}
else
for（i=0;i {
w[i]=1.0; //RectangularWindowFunction
}
fp1=fopen（fil1,"w"）;
//inputdata
for（i=0;i {t=i*dt;
xr[i]=sin（2*PI*f1*t）+sin（2*PI*f2*t）+sin（2*PI*f3*t）;
xr[i]=xr[i]*w[i];
fprintf（fp1,"%10.4f%12.4f\n",t,xr[i]）;
}
fclose（fp1）;
//fillzero
if（np {for（i=np;i xr[i]=0;
}
for（i=0;i xi[i]=0.0;
//calculatefft
fft（xr,xi,k,1）;
//outputspectrum
fp2=fopen（fil2,"w"）;
for（i=0;i {f=i*df;
fprintf（fp2,"%10.4f%12.4f\n",f,sqrt（xr[i]*xr[i]+xi[i]*xi[i]））;
}
fclose（fp2）;
free（xr）;free（xi）;free（w）;
}
2、参数文件
xt5_1.txt
xt5_amp_1.txt
200 //点数可增加或减少
4.0
10.035.0 70.0
1 //1—是Blackman窗，改成其他值就是矩形函数。
4．窗函数的描述
1）布莱克曼窗
C BLACKMANWINDOW
PARAMETER（N=201）
REALW（N）
PI=3.1415926
OPEN（1,FILE='BLACKMAN.DAT',FORM='FORMATTED',STATUS='UNKNOWN'）
DO10I=1,N
T=2*PI*（I-1）/（N-1）
10 W（I）=0.42-0.5*COS（T）+0.08*COS（2*T）
DO20I=1,N
20 WRITE（1,111）I-1,W（I）
111 FORMAT（I5,F8.4）
CLOSE
（1）
END
2）海宁窗
C HANNINGWINDOW
PARAMETER（N=201）
REALW（N）
PI=3.1415926
OPEN（1,FILE='HANNING.DAT',FORM='FORMATTED',STATUS='UNKNOWN'）
DO10I=1,N
10 W（I）=0.5*（1-COS（2*PI*（I-1）/（N-1）））
DO20I=1,N
20 WRITE（1,111）I-1,W（I）
111 FORMAT（I5,F8.4）
CLOSE
（1）
END
3）汉明窗
C HAMMINGWINDOW
PARAMETER（N=201）
REALW（N）
PI=3.1415926
OPEN（1,FILE='HAMMING.DAT',FORM='FORMATTED',STATUS='UNKNOWN'）
DO10I=1,N
10 W（I）=0.54-0.46*COS（2*PI*（I-1）/（N-1））
DO20I=1,N
20 WRITE（1,111）I-1,W（I）
111 FORMAT（I5,F8.4）
CLOSE
（1）
END
5．频率域数字滤波程序——实验六
1、程序
#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#definePI3.1415926
#include"fft_sub.c"
voidmain（）
{voidfft（）;
float*xr,*xi;
float*H;
inti,np,nfft,k;
floatt,dt,df,f,z,fc1,fc2;
FILE*fpar,*fp1,*fp2,*fp3,*fp4,*fp5;
charfil1[80],fil2[80],fil3[80],fil4[80],fil5[80];
fpar=fopen（"filter.par","r"）;
fscanf（fpar,"%s",fil1）;
fscanf（fpar,"%s",fil2）;
fscanf（fpar,"%s",fil3）;
fscanf（fpar,"%s",fil4）;
fscanf（fpar,"%s",fil5）;
printf（"%s\n",fil1）;
printf（"%s\n",fil2）;
printf（"%s\n",fil3）;
printf（"%s\n",fil4）;
printf（"%s\n",fil5）;
fscanf（fpar,"%d",&np）;
printf（"np=%d\n",np）;
fscanf（fpar,"%f",&dt）;
printf（"dt=%8.3fms\n",dt）;
dt=（float）（dt/1000.0）;
fscanf（fpar,"%f%f",&fc1,&fc2）;
printf（"fc1=%5.1f fc2=%5.1f\n",fc1,fc2）;
fclose（fpar）;
//calculatefftpoint
k=（int）（log（np）/log
（2））;
if（np>pow（2,k））k=k+1;
nfft=（int）pow（2,k）;
df=（float）（1.0/（nfft*dt））;
printf（"nfft=%d k=%d\n",nfft,k）;
printf（"dt=%8.4f df=%8.4f\n",dt,df）;
//allocatememory
xr=（float*）calloc（nfft,sizeof（float））;
xi=（float*）calloc（nfft,sizeof（float））;
H=（float*）calloc（nfft,sizeof（float））;
//generateafilter
fp1=fopen（fil1,"w"）;
for（i=0;i<=nfft/2;i++）
{f=i*df;
if（f<=fc2&&f>=fc1）H[i]=1.0;
elseH[i]=0.0;
fprintf（fp1,"%4d%10.4f%12.4f\n",i,f,H[i]）;
}
for（i=nfft/2+1;i {f=i*df;
H[i]=H[nfft-i];
fprintf（fp1,"%4d%10.4f%12.4f\n",i,f,H[i]）;
}
fclose（fp1）;
printf（"TheFilterisok!
\n"）;
fp2=fopen（fil2,"r"）;
//inputdata
for（i=0;i {fscanf（fp2,"%f%f",&t,&z）;
xr[i]=z;
}
fclose（fp2）;
printf（"TheDataisok!
\n"）;
//fillzero
if（np {for（i=np;i xr[i]=0;
}
for（i=0;i xi[i]=0.0;
//calc
升级会员 