I have been trying to write a simple program to perform an fft on a 1D input array using fftw3. Here I am using a seismogram as an input. The output array is, however, coming out to contain only zeroes.
I know that the input is correct as I have tried doing the fft of the same input file in MATLAB as well, which gives correct results. There is no compilation error. I am using f95 to compile this, however, gfortran was also giving pretty much the same results. Here is the code that I wrote:-
program fft
use functions
implicit none
include 'fftw3.f90'
integer nl,row,col
double precision, allocatable :: data(:,:),time(:),amplitude(:)
double complex, allocatable :: out(:)
integer*8 plan
open(1,file='test-seismogram.xy')
nl=nlines(1,'test-seismogram.xy')
allocate(data(nl,2))
allocate(time(nl))
allocate(amplitude(nl))
allocate(out(nl/2+1))
do row = 1,nl
read(1,*,end=101) data(row,1),data(row,2)
amplitude(row)=data(row,2)
end do
101 close(1)
call dfftw_plan_dft_r2c_1d(plan,nl,amplitude,out,FFTW_R2HC,FFTW_PATIENT)
call dfftw_execute_dft_r2c(plan, amplitude, out)
call dfftw_destroy_plan(plan)
do row=1,(nl/2+1)
print *,out(row)
end do
deallocate(data)
deallocate(amplitude)
deallocate(time)
deallocate(out)
end program fft
The nlines() function is a function which is used to calculate the number of lines in a file, and it works correctly. It is defined in the module called functions.
This program pretty much tries to follow the example at http://www.fftw.org/fftw3_doc/Fortran-Examples.html
There might just be a very simple logical error that I am making, but I am seriously unable to figure out what is going wrong here. Any pointers would be very helpful.
This is pretty much how the whole output looks like:-
.
.
.
(0.0000000000000000,0.0000000000000000)
(0.0000000000000000,0.0000000000000000)
(0.0000000000000000,0.0000000000000000)
(0.0000000000000000,0.0000000000000000)
(0.0000000000000000,0.0000000000000000)
.
.
.
My doubt is directly regarding fftw, since there is a tag for fftw on SO, so I hope this question is not off topic
As explained in the comments first by @roygvib and @Ross, the plan subroutines overwrite the input arrays because they try the transform many times with different parameters. I will add some practical use considerations.
You claim you do care about performance. Then there are two possibilities:
FFTW_MEASURE. The planning subroutine is many times slower than actual plan execute subroutine. Use FFTW_ESTIMATE and it will be much faster.FFTW_MEASURE tells FFTW to find an optimized plan by actually computing several FFTs and measuring their execution time. Depending on your machine, this can take some time (often a few seconds). FFTW_MEASURE is the default planning option.
FFTW_ESTIMATE specifies that, instead of actual measurements of different algorithms, a simple heuristic is used to pick a (probably sub-optimal) plan quickly. With this flag, the input/output arrays are not overwritten during planning.
http://www.fftw.org/fftw3_doc/Planner-Flags.html