SU初学习之程序编写

示例一：

建立一个采样间隔为4ms，200个采样点，主频30Hz的ricker子波。

（一）建立ricker.c文件

首先，在命令行中输入vi  ricker.c，建立一个名为ricker.c的文件，在文件中加入头文件，程序如下：

#include<stdio.h>
#include<math.h>
#include<su.h>
#include<segy.h>
main()
{
    float pi=3.1415926,a,b,t=0.0,lt=0.8,dt=0.004,fp=30,*r1;#lt为总检测时间，dt为每隔dt时间检测一次，fp为主频
    int i,nt;
    nt=(int)(lt/dt);#nt为采样点个数
    FILE *fp1;#创建文件
    r1=alloc1float(nt);#为r1分配空间
    for(i=0;i<nt;i++)
        r1[i]=0.0;#安全起见，先为数组每个元素赋0值
    for(i=0;i<nt;i++)#雷克子波的生成方式
    {
        t=(i-nt/2)*dt;
        a=pi*t*fp;
        b=a*a;
        r1[i]=(1-2*b)*exp(-b);
    }
    fp1=fopen("r1.bin","wb");#打开文件，一个fopen对应一个fclose，wb为打开方式
    for(i=0;i<nt;i++)#向文件中写入数据
         fwrite(&r1[i],sizeof(float),1,fp1);
    fclose(fp1);
}

（二）为ricker.c建立makefile文件

在命令行输入 vi makefile后，再输入

#! /bin/sh
# Make  ricker  programs
INC=${CWPROOT}/include
LIK=${CWPROOT}/lib
LIB=-lsu -lpar -lcwp -lm
ricker:ricker.c
gcc ricker.c -o ricker  -O3 -I${INC} -L${LIK} ${LIB}

保存退出。

注意：makefile与ricker.c在同一个文件夹里。

（三）生成可执行文件

在命令行中输入：make

（四） 运行程序

输入./ricker

使用xgraph<r1.bin n1=200 d1=0.004&

可得到如下图所示的ricker子波。

示例二

生成ricker子波后，对其进行傅里叶变换，计算它的振幅谱与相位谱。

（一）建立一个名为fft.c的文件（其中ricker子波参数与示例一相同），程序如下：

#include "su.h"
#include "segy.h"
#include "header.h"
#include <time.h>
void ricker(float *r1,int nt)
{
	float pi=3.1415926,a,b,t=0.0,lt=0.8,dt=0.004;
	int i;
	float fp=30;
	for(i=0;i<nt;i++)
		r1[i]=0.0;
	for(i=0;i<nt;i++)
	{
		t=(i-nt/2)*dt;
		a=pi*t*fp;
		b=a*a;
		r1[i]=(1-2*b)*exp(-b);
	}
}
main()
{
	int nt,ntfft,nw,i,iw;
	float dt=0.004,lt=0.8;
	float *r1,*amp,*pha,*a,*b,*r2f;
	complex *r1f;
	FILE *fp1,*fp2,*fp3,*fp4;
	nt=(int)(lt/dt);
	r1=alloc1float(nt);
	for(i=0;i<nt;i++)
		r1[i]=0.0;
		ricker(r1,nt);
		fp1=fopen("r1.bin","wb");
	for(i=0;i<nt;i++)
    	fwrite(&r1[i],sizeof(float),1,fp1);
	fclose(fp1);
	ntfft=npfar(nt);
	nw=(int)(ntfft/2+1);
	r1f=alloc1complex(nw);
	amp=alloc1float(nw);
	pha=alloc1float(nw);
	a=alloc1float(nw);
	b=alloc1float(nw);
	for(iw=0;iw<nw;iw++)
	{
    	r1f[iw]=cmplx(0.0,0.0);
    	amp[iw]=0.0;
    	pha[iw]=0.0;
    	a[iw]=0.0;
    	b[iw]=0.0;
	}

	pfarc(-1,ntfft,r1,r1f);
	for(iw=0;iw<nw;iw++)
	{
    	a[iw]=r1f[iw].r;
    	b[iw]=r1f[iw].i;
	}
	for(iw=0;iw<nw;iw++)
	{	
    	amp[iw]=sqrt(a[iw]*a[iw]+b[iw]*b[iw]);
    	pha[iw]=atan(b[iw]/a[iw]);
	}
	fp2=fopen("amp.bin","wb");
	for(i=0;i<nw;i++)
    	fwrite(&amp[i],sizeof(float),1,fp2);
	fclose(fp2);
	fp3=fopen("pha.bin","wb");
	for(i=0;i<nw;i++)
    	fwrite(&pha[i],sizeof(float),1,fp3);
	fclose(fp3);
}

 (二) 为fft.c建立makefile文件

在命令行输入 vi  makefile后，再输入

#! /bin/sh
INC=${CWPROOT}/include
LIK=${CWPROOT}/lib
LIB=-lsu -lpar -lcwp -lm
fft:fft.c
gcc fft.c -o fft  -O3 -I${INC} -L${LIK} ${LIB}

保存退出。

注意：makefile与fft.c在同一个文件夹里。

（三）生成可执行文件

在命令行中输入：make

（四） 运行程序

输入./fft

使用

xgraph<amp.bin n1=200 d1=0.004&

xgraph<pha.bin n1=200 d1=0.004&

可得到如下振幅谱与相位谱：