1D convolution of real data and real response


kurti
10-26-2006, 07:19 PM
Hi all,

I 'just' need to convolve a real signal (length n) with a response (length n). The main program is F90.

In Matlab, the following produces the exact results (compared to the clumsy summation technique):

fft_d = fft(data)
fft_r = fft(response)
conv = fft_d .* fft_r
ic = real(ifft(conv))

Put simply, I want to do this in Fortran(90).

[1] NR's convlv does not work.
Somehow, the code gets problems with declarations (I use DOUBLE PRECISION in my f90 program). Somewhere along the computation, there's always a NaN in the convlv machinery.
I also tried zero-padding back and forth--nothing works.

[2] When I use four1 instead of Matlab's fft (and using REAL in the main program), results are getting better. But, without zero-padding, I get somewhat a doubled response (symmetric), with zero padding, the result looks QUITE good, except a systematic misfit between the summation and the convolution--which I have no clue how to get rid of. It's too big and nonlinear (but symmetric) to be tolerated and scaled, respectively.

[3] I tried the same with the fft.f90 from Alan Miller's page, but the result is quite the same as in the four1 case, except no hassle with declarations. (I also converted the convlv pack with to_f90.f90, but this didn't work out, too...)

I have the feeling it has something to do with the symmetry of the forward FFT's? Since Matlab's ifft checks for it, and gives back a real vector of length n, it may very well be that this is what I have to do before I use any fft routine to do the inverse fft.

Does anybody has an idea how this should look like? I don't get the answer by studying NR's convlv.

I appreciate

Kurti