Re: cfftf vs ezfftf

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Tue Feb 07 2012 - 18:16:13 MST

Hi Houjon

I don't think there is an ambiguity. cfft returns a two-sided spectrum.
The appropriate normalization is (1/N) ... not (2/N).

The attached script provides a simple test. It compares the
classic biased variance ( SUM(x^2)/N ) with the variances
computed via cfftf and ezfftf.

        N=24 xVar=193.207

cfftf: cf2sum =193.207
ezfftf: cfez_var=193.207

There are different ways ffts return the data.

One brief discussion is at:

http://www.wavemetrics.com/products/igorpro/dataanalysis/signalprocessing/powerspectra.htm

On 02/02/2012 03:30 PM, Houjun Wang wrote:
> Hi,
>
> It appears that there is an ambiguity with the normalization in 'cfftf':
>
> To obtain a normalized transform the output must be divided by "N", the
> size of the rightmost dimension.
>
> This works fine to obtain the 'average'; but when computing the power
> spectrum, the magnitude would be small by a factor of 'two' [for
> amplitude, and 2^2=4 for power spectrum], when compared with that from
> 'ezfftf'. See the example below - in the example, I used [N/2], instead
> of [N], just to match the power spectra, but then the average needs to
> be divided by 'two'.
>
> Please let us know whether this is correct [in this case, there will
> also be a 2^2 factors in the examples of the "space-time spectra" in
> Wheeler-Kiladis Space-Time Spectra that use 'cfftf'].
>
> Thanks.
>
> Houjun
>
> =======
> begin
>
> x = (/ 1002, 1017, 1018, 1020, 1018, 1027, \
> 1028, 1030, 1012, 1012, 982, 1012, \
> 1001, 996, 995, 1011, 1027, 1025, \
> 1030, 1016, 996, 1006, 1002, 982 /)
>
> N = dimsizes(x)
>
> df = 1.0/N
> xAvg = avg(x)
>
> print ("N="+N+" df="+df+" avg(x)="+xAvg)
>
> cf = cfftf (x, 0.0, 0) ; imaginary part set to 0.0
>
> printVarSummary(cf)
>
> cf = cf/N*2
>
> print(sprintf("%9.5f",cf@frq
> <https://web3.acd.ucar.edu/horde/imp/message.php?mailbox=INBOX.Sent&index=81#>)
> +" "+sprintf("%9.3f", cf(0,:))+" "+sprintf("%9.3f", cf(1,:)) )
>
> cf_ezfftf = ezfftf (x)
>
> print(cf_ezfftf)
>
> power = new ((/N/), "float")
> power = sqrt((cf(0,:)^2 + cf(1,:)^2))
>
> print(power)
>
> power_ezfftf = new ((/N/2/), "float")
>
> power_ezfftf = sqrt((cf_ezfftf(0,:)^2 + cf_ezfftf(1,:)^2))
>
> print(power_ezfftf)
>
> end
> =========
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk

Received on Tue Feb 7 18:16:21 2012

This archive was generated by hypermail 2.1.8 : Thu Feb 09 2012 - 13:33:26 MST