PDA

View Full Version : problems with FFT / FFTCOMP



hcrisp
11-25-2009, 12:32 PM
A user alerted me to the fact that use of FFT breaks in one case, whereas use of FFTCOMP breaks in another case. Looks like FFT has numerical limitations but FFTCOMP has memory limitations. Neither is fully-functional, and therefore PV-WAVE does not have a single FFT function that works for all cases. Let me illustrate.

Example where FFTCOMP works but FFT does not:


kernel_length = 4001L
h = FLTARR(kernel_length ) + 1.

h(0:2) = 0.
h(3) = 0.0300344
h(4) = 0.412832
h(5) = 0.882798
h(38) = 0.754520
h(39) = 0.254548
h(40) = 2.74063e-005
h(41) = 0.245479
h(42) = 0.745451
h(43) = 0.999973
h(64) = 0.972759
h(65) = 0.595404
h(66) = 0.122645
h(67:80) = 0.
h(81) = 0.240987
h(82) = 0.740877
h(83) = 0.999890
h(3918) = 0.999890
h(3919) = 0.740877
h(3920) = 0.240987
h(3921:3934) = 0.
h(3935) = 0.122645
h(3936) = 0.595404
h(3937) = 0.972759
h(3958) = 0.999973
h(3959) = 0.745451
h(3960) = 0.245479
h(3961) = 2.74063e-005
h(3962) = 0.254548
h(3963) = 0.754520
h(3996) = 0.882798
h(3997) = 0.412832
h(3998) = 0.0300344
h(3999:4000) = 0.

hinv_1 = FFT(h, 1) / kernel_length
hout_1 = FFT(hinv_1, -1)
hinv_2 = FFTCOMP(h, /BACKWARD, /COMPLEX) / kernel_length
hout_2 = FFTCOMP(hinv_2, /COMPLEX) / kernel_length

PLOT, ABS(hout_1)
OPLOT, ABS(hout_2), COLOR=230

; Notice the the FFTCOMP result is numerically correct, while the FFT result is not!


Example where FFT works but FFTCOMP does not:


kernel_length = 12000000L
h = DOUBLE(RANDOMU(s, kernel_length ))
hinv_1 = FFT(h, 1) / kernel_length
hout_1 = FFT(hinv_1, -1)
; Works!

kernel_length = 12000000L
h = DOUBLE(RANDOMU(s, kernel_length ))
hinv_2 = FFTCOMP(h, /BACKWARD, /COMPLEX) / kernel_length
;%%%FFTCOMP: Terminal error: MATH_OUT_OF_MEMORY_1
; Not enough memory, with "n" = 12000000.
;%%%
;%%%Execution halted at FFTCOMP <fftcomp.pro( 134)> (MATHSTAT_140).
;%%%Called from $MAIN$ .


Can anything be done to make both the FFT and/or FFTCOMP functions work for all cases? (I am using WAVE 8.00 on Windows XP.)

rwagner
11-30-2009, 10:10 AM
Hi Hcrisp,
FFT will give correct results using double input datatypes. If you need floats later on, you can cast the returned values down after the call to FFT. You need to use the double precision since your kernel length is a large prime factor (4001).

So change your line:
h = FLTARR(kernel_length ) + 1. to read:
h = DBLARR(kernel_length ) + 1.


You're correct WRT FFTCOMP, I see the out of memory error as well on a system with large amounts of memory. I'll file a change request for this.

-Ryan