| wavenum, frequency FFT plot [message #66173] |
Thu, 23 April 2009 13:57  |
kishore1818
Messages: 24
Registered: June 2007
|
Junior Member |
|
|
Hi,
I have time vs latitudinal values and I would like to estimate the
power spectral density 2-D do the 2-D FFT to get the frequency VS wave
number .
Example : data sets are (30,8); 30=days values and 8 latitudinal avg.
data
I am planning to wave number, frequency and power contour plot.
x-axis zonal wave number range is -7 to 7
Y-axis frequency range is -1. To 1.
Does anybody have a good example program on 2D FFT in IDL program.
Thanking you,
Kishore
|
|
|
|
| Re: wavenum, frequency FFT plot [message #66221 is a reply to message #66173] |
Wed, 29 April 2009 10:49  |
kishore1818
Messages: 24
Registered: June 2007
|
Junior Member |
|
|
On Apr 29, 9:38 am, "R.G. Stockwell" <noemai...@please.com> wrote:
> <kishore1...@gmail.com> wrote in message
>
> news:6157cb9b-a549-4b42-8845-798916fbda9b@c18g2000prh.google groups.com...
> On Apr 23, 3:03 pm, "R.G. Stockwell" <noemai...@please.com> wrote:
> Hi Bob,
> Thanks for your e-mail message.
> Actually I have 365 days daily avg values, every 10 deg latitudinal
> averages.
> ie., (365,18). I need to calculate the psd for 30days data sets for
> all latitudinal values.
> But I am planning to plot contour: wvae no. (x axis), frequency (y-
> axis) and z is power.
> could you give a example program then it is very easy to understand to
> me and also I am beginner in IDL.
>
> Thanking you,
> Kishore
>
> Here you go:
>
> ndays = 365
> nlats = 18
>
> ; make fake data with DC component
> data = randomn(seed,ndays,nlats) + 0.2
>
> spe =fft(data) ; 2D fft
>
> ; make the frequency and wavenumber arrays
> freqs = findgen(ndays)/ndays
> ; put in negative freqs
> freqs[ndays/2+1:*] = -1.0 + freqs[ndays/2+1:*]
>
> wavenumbers = findgen(nlats)
> ; put in negative wavenumbers
> wavenumbers[nlats/2+1:*] = -nlats + wavenumbers[nlats/2+1:*]
>
> ;FFT has postive freqs first
> ; (0, 1/NT, 2/NT up to nyquist, then most negative
> ; and then counting towards zero)
>
> ;if you want, you can switch it to the "normal way"
> ; of viewing the spectr, with 0,0 at the middle.
> ; so shift the spectra.
> spe = shift(spe,ndays/2, nlats/2-1) ; note ndays is odd, nlats is even
> ;shift the freqs
> freqs = shift(freqs,ndays/2)
> ;shift the wavenumbers
> wavenumbers = shift(wavenumbers,nlats/2-1)
>
> ; ALSO you want the wavenumber on the x axis, so transpose
> spe = transpose(spe)
>
> ; make levels for contour plot
> nlevels = 40
> levels =
> findgen(nlevels+2)/nlevels*(max(abs(spe))-min(abs(spe)))+min (abs(spe),/nan)
> contour,abs(spe),wavenumbers,freqs,levels=levels,/fill,$
> xtitle='Meridional Wavenumber',ytitle='Frequency (1/day)'
>
> end
Thanks for providing the detailed code.
Kishore
|
|
|
|
| Re: wavenum, frequency FFT plot [message #66224 is a reply to message #66173] |
Wed, 29 April 2009 09:38 |
R.G. Stockwell
Messages: 363
Registered: July 1999
|
Senior Member |
|
|
<kishore1818@gmail.com> wrote in message
news:6157cb9b-a549-4b42-8845-798916fbda9b@c18g2000prh.google groups.com...
On Apr 23, 3:03 pm, "R.G. Stockwell" <noemai...@please.com> wrote:
Hi Bob,
Thanks for your e-mail message.
Actually I have 365 days daily avg values, every 10 deg latitudinal
averages.
ie., (365,18). I need to calculate the psd for 30days data sets for
all latitudinal values.
But I am planning to plot contour: wvae no. (x axis), frequency (y-
axis) and z is power.
could you give a example program then it is very easy to understand to
me and also I am beginner in IDL.
Thanking you,
Kishore
Here you go:
ndays = 365
nlats = 18
; make fake data with DC component
data = randomn(seed,ndays,nlats) + 0.2
spe =fft(data) ; 2D fft
; make the frequency and wavenumber arrays
freqs = findgen(ndays)/ndays
; put in negative freqs
freqs[ndays/2+1:*] = -1.0 + freqs[ndays/2+1:*]
wavenumbers = findgen(nlats)
; put in negative wavenumbers
wavenumbers[nlats/2+1:*] = -nlats + wavenumbers[nlats/2+1:*]
;FFT has postive freqs first
; (0, 1/NT, 2/NT up to nyquist, then most negative
; and then counting towards zero)
;if you want, you can switch it to the "normal way"
; of viewing the spectr, with 0,0 at the middle.
; so shift the spectra.
spe = shift(spe,ndays/2, nlats/2-1) ; note ndays is odd, nlats is even
;shift the freqs
freqs = shift(freqs,ndays/2)
;shift the wavenumbers
wavenumbers = shift(wavenumbers,nlats/2-1)
; ALSO you want the wavenumber on the x axis, so transpose
spe = transpose(spe)
; make levels for contour plot
nlevels = 40
levels =
findgen(nlevels+2)/nlevels*(max(abs(spe))-min(abs(spe)))+min (abs(spe),/nan)
contour,abs(spe),wavenumbers,freqs,levels=levels,/fill,$
xtitle='Meridional Wavenumber',ytitle='Frequency (1/day)'
end
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
