Re: IDL FFT (spec -> interferogram) [message #30031 is a reply to message #30029] |
Thu, 04 April 2002 08:32   |
Randall Skelton
Messages: 169 Registered: October 2000
|
Senior Member |
|
|
It is always the case... I struggled for hours with the IDL FFT routine and
the numerical recipes cook book and now that I have written the 'please
help' request to the newsgroup, I think I've found it.
Read in the data:
IDL> spec = dcomplexarr(2048)
IDL> read_cmplx, 'spec2048.dat', spec
IDL> igm = dcomplexarr(2048)
IDL> read_cmplx, 'igm2048.dat', igm
Fourier Transform and unwrap the phase
IDL> for i=0, n_elements(spec)-1, 2 do spec[i] = -spec[i]
IDL> idl_igm = fft(spec)
IDL> for i=0, n_elements(idl_igm)-1, 2 do idl_igm[i] = -idl_igm[i]
IDL> plot, idl_igm
(NB: the for loops can be replaced by a where statement and index that
gabs an even number sequence but the above is a little more clear)
There is no need to fold the spectrum about the Nyquist frequency and
therefore the above does not require a division by 2 amplitude correction.
Cheers,
Randall
On Thu, 4 Apr 2002, Randall Skelton wrote:
> Hi all,
>
> Having read through all of the FFT posts that google groups keeps, I am no
> closer to understanding why I am unable to transform a spectrum into an
> interferogram using IDL. All of the data files, procedures, and pictures
> of this are at http://tulip.atm.ox.ac.uk/~rhskelto/fft-help/
>
> Given two files:
>
> 1) 'spec.dat' contains 512 points of complex spectral data
>
> 2) 'igm.dat' contains 512 points of complex interferogram data that was
> derived from 'spec.dat' using a prime factor FFT written in C. This is
> the correct interferogram as far as I am concerned. The plot
>
> Read in the data:
>
> IDL> spec = dcomplexarr(512)
> IDL> read_cmplx, 'spec.dat', spec
> IDL> igm = dcomplexarr(512)
> IDL> read_cmplx, 'igm.dat', igm
>
> Plot the expected result:
>
> IDL> plot, igm
> IDL> write_jpeg,'igm.jpg',tvrd()
>
> Do the Fourier Transform in IDL (based on Paul van Delst's examples):
>
> IDL> spec2 = temporary( [ spec, reverse( spec[ 1: n_elements(spec) - 2 ] ) ] )
> IDL> idl_igm = fft(temporary(spec2), /double, /inverse)
> IDL> idl_igm = shift(idl_igm, -1 * (n_elements(spec)-1))
>
> Plot the IDL result:
>
> IDL> plot, idl_igm
> IDL> write_jpeg, 'idl_igm.jpg', tvrd()
>
>
> The result 'idl_igm' contains twice the number of points (minus 2)
> because of the required reflection about the Nyquist frequency.
> Moreover, the result appears to be modulated (almost like a frequency
> chirp)? I recall having a similar problem with a 2pi phase-wrapping in
> MathCad a number of years ago that gave similar results but I cannot
> remember how to fix it. I also cannot seem to reproduce the AIRS
> interferograms shown on Paul's site...
>
> My question is, how do I get the desired result (i.e. 'igm.jpg') in IDL?
>
> Cheers,
> Randall
>
> IDL Version 5.3, Linux RH 7.1
>
>
|
|
|