0%

如何调用fftw进行Fast Fourier Transform

  我们常常要看信号的振幅谱来进行分析,那傅里叶变换就必不可少。如果水平不错你可以试着自己写。当然有很多已经写好的包,非常方便,例如这里要讲到的fftw[1]

  他的用法其实人家主页说明文档讲的很清楚。我这里就记录一下怎么用之来对sac文件进行读取并计算fft。程序是用fortran写的,C的话可以参考fftw说明文档。
首先把主程序贴出来:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
program main
use globe_data // 全局变量
use sacio // sac 头变量
implicit none
integer i,nerr,itest
character (180) :: sacfile,tmp
type(sac_head) :: sachead
complex(8),allocatable,dimension(:) :: s1,s2
if (iargc().ne.2)stop 'Usage: fft sacfile '
call getarg(1,sacfile) // 读参数到sacfile,即sac 文件
call read_sachead(trim(sacfile),sachead,nerr) //读sac文件头变量
npts=sachead%npts
dt=sachead%delta
nn=2
npow=1
do while(nn.le.npts)
nn=nn*2
npow=npow+1
enddo
nk=nn/2+1
halfn=20
dom=dble(1.0/nn/dt)
allocate(sig(nn,3))
call read_sac(trim(sacfile), sig(:,1),sachead,nerr) //读sac文件数据到sig(:,1)
allocate(s1(nn),s2(nn))
s1=czero
s1(1:npts)=dcmplx(dble(sig(1:npts,1)),0.0d0)
call dfftw_plan_dft_1d(plan,nn,s1,s2,FFTW_FORWARD, FFTW_ESTIMATE)
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)
sachead%npts=nk
sig(:,1)=real(dreal(s2))
call write_sac(trim(sacfile)//'_fft',sig(:,1),sachead,nerr) ! problem with nerr=-1
deallocate(sig,s1,s2)
end program

这里sig(:,:)是个二维数组,其实用一维的就够了哈。
sacio.f90 是一个module,定义了sac文件的头,并含有sac读写程序。需要的话给我发信息,或者邮件^_^。
globe_data.f90也是一个module,定义了全局变量:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
module globe_data
integer,parameter :: nmax=2000000,nstmax=1000
real(4),dimension(4,2):: fre
real(4),allocatable,dimension(:,:):: sig,sigo
real(4),allocatable,dimension(:) :: sigt
real :: dt
real(8) :: dom
integer :: halfn
integer :: ncom,comb,npts
integer :: nn,npow,nk,nf
complex(8),allocatable,dimension(:,:):: seisout
complex(8),parameter:: czero=(0.0d0,0.0d0)
integer,parameter :: FFTW_ESTIMATE=64,FFTW_MEASURE=1
integer,parameter :: FFTW_FORWARD=-1,FFTW_BACKWARD=1
integer(8) :: plan,plan1,plan2,plan3
end module

编译需要一个makefile:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
FC=gfortran
FFLAG=-lfftw3 -I /usr/include -L /usr/lib64 -fbounds-check
objects=call_fft.o sacio.o globe_data.o
all:sacio.mod globe_data.mod call_fft
.f.o:
$(FC) $(FFLAG) $< -c
%.o:%.f90
$(FC) $(FFLAG) $< -c
sacio.mod:sacio.f90
$(FC) $< -c
globe_data.mod:globe_data.f90
$(FC) $< -c
call_fft:$(objects)
$(FC) $^ -o $@ $(FFLAG) -lm
clean:
-rm *.o *.mod

注意这里要加上-lfftw3表示调用fftw,-I给定fftw头的路径,-L给定fftw的lib。如果是按照默认路径安装的fftw,那一般都不用指定-I和-L,因为默认路径一般已经包含在编译环境里了。
好了大功告成了,但不保证拷贝下来就能运行通过哦。有问题发邮件吧,我们再交流^_^


  1. 1.http://www.fftw.org/
有疑问或建议请发邮件到junxie01@gmail.com。