| 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
