| 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
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
