Re: [ferret_users] power spectrum 1d integrate to GMT

Dear Russel Fiedler,

Thank you for the replied.
Yes, I've tried the external FFT functions, and the results is like the attached figures.
Actually I need the graph in log scale for y-axes (not in meter as in figure attached).
Following is my script by applying FFT functions.

!==start

          cancel win/all

          cancel mode logo

          use "/home/dessy/Documents/verification.cdf"

          set win 1

          plot/nolab/line=1/dash/l=1897:2640/thick=2 simulation

          plot/nolab/overlay/line=1/thick=2 verification

          !==

          set win 2

          let FFT_simulation=simulation

          let FFT_simulationfft=FFTA(FFT_simulation[l=1897:2640])

          set variable/title="Amplitude Spectrum" FFT_simulationfft

          plot/line=1/dash FFT_simulationfft

          let FFT_verification=verification

          let FFT_verificationfft=FFTA(FFT_verification[l=1:744])

          set variable/title="Amplitude Spectrum" FFT_verificationfft

          plot/overlay/line=1/thick=3/nolab FFT_verificationfft

          !==

          set win 3

          LET FFT_nf = `FFT_simulationfft,return=lend`

          LET FFT_nyquist = 0.5

          LET FFT_freq1 = FFT_nyquist/ FFT_nf

          DEFINE AXIS/T=`FFT_freq1`:`FFT_nyquist`:`FFT_freq1` faxis

          DEFINE GRID/T=faxis gfftfreq

          LET a = T[g=gfftfreq]

          Let per = 1./a

          !==

          LET FFT_nf2 = `FFT_verificationfft,return=lend`

          LET FFT_nyquist2 = 0.5

          LET FFT_freq2 = FFT_nyquist2/ FFT_nf2

          DEFINE AXIS/T=`FFT_freq2`:`FFT_nyquist2`:`FFT_freq2` faxis2

          DEFINE GRID/T=faxis2 gfftfreq2

          LET b = T[g=gfftfreq2]

          Let per2 = 1./b

          !==

          PLOT/VS/dash/line=1/HLIMITS=0:40:2/thick=3/nolab/TITLE="Amplitude
          Spectrum" per[l=1:`FFT_nf`], FFT_simulationfft

          PLOT/OVERLAY/VS/LINE=1/thick=1/HLIMITS=0:40:2/nolab
          per2[l=1:`FFT_nf2`], FFT_verificationfft

          !PPL XFOR (I2)

          !PPL XLAB Period, hours

          !PPL YLAB 

          !PPL PLOT

          !zoom

          set win 4

          PLOT/VS/dash/line=1/nolab/HLIMITS=10:28:2/thick=3/TITLE="Amplitude
          Spectrum" per[l=1:`FFT_nf`], FFT_simulationfft

          PLOT/OVERLAY/VS/nolab/LINE=1/thick=1/HLIMITS=10:28:2
          per2[l=1:`FFT_nf2`], FFT_verificationfft

          !end of script

I've failed search in the previous Ferret Email Archives which similarly to my case.

I am looking forward to hearing any further suggestion about this.

Thank you very much.

Best wishes,

dessy

From: Russ Fiedler <russell.fiedler@xxxxxxxx>
To: dessy berlianty <dessyberlianty@xxxxxxxxx>
Cc: Ferret Users <ferret_users@xxxxxxxx>
Sent: Tuesday, September 4, 2012 4:07 PM
Subject: Re: [ferret_users] power spectrum 1d integrate to GMT

Hi dessy,

Have you tried the external FFT functions that come with Ferret?

--------------

FFT_AMP(A)
Computes fft amplitude spectra, normalized by 1/N, WITH CORRECTED
FREQ AXIS
A: Variable with regular time axis. Specify time explicitly e.g.
fft_amp(var[l=1,96])
FFTA_SAMPLE(A)
Computes fft amplitude spectra, normalized by 1/N
A: Variable with regular time axis. Specify time explicitly e.g.
ffta_sample(var[l=1,120])
FFTA(A)
Computes fft amplitude spectra, normalized by 1/N
A: Variable with regular time axis. Specify time explicitly e.g.
ffta(var[l=1,120])
FFT_PHAS(A)
Computes fft phase, WITH CORRECTED FREQ AXIS
A: Variable with regular time axis. Specify time explicitly e.g.
fft_phas(v[l=1,96])
FFTP(A)
Computes fft phase
--------------

Cheers,

Russ

On Tue, 2012-09-04 at 15:13 +1000, dessy berlianty wrote:
> Dear Ferret Users,
>
> I need some suggestions to solve my problem.
> In order to make log-log scale of power spectrum 1d (based on time
> series), I have been tried gmt_power.jnl as suggestion by Steve Hankin
> in this link.
> After installation of GMT, I loss my way to continue the execution of
> gmt_power.jnl because of the following errors:
>
>
> yes? set data monthly_navy_winds
> yes? go gmt_power uwnd[x=180,y=0] uwnd x180y0 64
> LISTing to file x180y0.tseries
> >>> Using GMT command: spectrum1d x180y0.tseries -Nx180y0 -S64
> -D30.4375
> sh: 1: spectrum1d: not found
> **TMAP ERR: non-existent or not on line
> x180y0.xpower
> SET DATA/EZ/var=freq x180y0.xpower
> Command file, command group, or REPEAT execution aborted
> yes?
>
>
>
> Thank you in advance for your time, and highly appreciate for any
> helps.
>
>
>
> Best wishes,
>
>
> dessy
>