Re: cfftf vs ezfftf

Thanks, Dennis.

So you're doing 'sum' all.

But if you do a test for a 'sine-wave', you'll see the difference:
Should we also need to do 'sum' here, i.e., sqrt((cf(0,1)^2+cf(1,1)^2)
+ (cf(0,23)^2+cf(1,23)^2)) for wave-amplitude for wave#1? But this does
not appears to be the way used in the NCL space-time spectra analysis.
Did I miss some here about waves resolved to eastward and westward
propagating.

Any comments for this?

Thanks.

Houjun

====
Here is the sine-wave test script and output: 'cfftf' got the amplitude
of about '1' [0.9754332] for wave 1 only if "N/2" is used for
normalization [but then the 'ave' is 20, which should be 10].

Houjun

========
begin

N = 24

x = new ((/N/), "float")

pi = 4.0*atan(1.0)
print(pi)

do i = 0, N-1
x(i) = 10. + sin(2.*pi*i/(N-1))
end do

print(x)

; cfftf

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) +" "+sprintf("%9.3f", cf(0,:))+"
"+sprintf("%9.3f", cf(1,:)) )

; ezfftf

cf_ezfftf = ezfftf (x)

print(cf_ezfftf)

; amplitude

amplitude = new ((/N/), "float")
amplitude = sqrt((cf(0,:)^2 + cf(1,:)^2))

print(amplitude)

amplitude_ezfftf = new ((/N/2/), "float")

amplitude_ezfftf = sqrt((cf_ezfftf(0,:)^2 + cf_ezfftf(1,:)^2))

print(amplitude_ezfftf)

end
=======

v1n3 338% ncl cfftf_vs_ezfftf-for-sine-wave.ncl
  Copyright (C) 1995-2010 - All Rights Reserved
  University Corporation for Atmospheric Research
  NCAR Command Language Version 5.2.0
  The use of this software is governed by a License Agreement.
  See http://www.ncl.ucar.edu/ for more details.

Variable: pi
Type: float
Total Size: 4 bytes
             1 values
Number of Dimensions: 1
Dimensions and sizes: [1]
Coordinates:
(0) 3.141593

Variable: x
Type: float
Total Size: 96 bytes
             24 values
Number of Dimensions: 1
Dimensions and sizes: [24]
Coordinates:
Number Of Attributes: 1
   _FillValue : -999
(0) 10
(1) 10.2698
(2) 10.51958
(3) 10.73084
(4) 10.88789
(5) 10.97908
(6) 10.99767
(7) 10.94226
(8) 10.81697
(9) 10.63109
(10) 10.3984
(11) 10.13617
(12) 9.863833
(13) 9.601599
(14) 9.368912
(15) 9.18303
(16) 9.057739
(17) 9.002331
(18) 9.020916
(19) 9.112115
(20) 9.269164
(21) 9.480416
(22) 9.730204
(23) 10
(0) N=24 df=0.0416667 avg(x)=10

Variable: cf
Type: float
Total Size: 192 bytes
             48 values
Number of Dimensions: 2
Dimensions and sizes: [2] x [24]
Coordinates:
Number Of Attributes: 2
   npts : 24
   frq : <ARRAY of 24 elements>
(0) 0.00000 20.000 0.000
(1) 0.04167 0.127 -0.967
(2) 0.08333 -0.016 0.058
(3) 0.12500 -0.013 0.031
(4) 0.16667 -0.012 0.021
(5) 0.20833 -0.012 0.015
(6) 0.25000 -0.012 0.012
(7) 0.29167 -0.012 0.009
(8) 0.33333 -0.012 0.007
(9) 0.37500 -0.011 0.005
(10) 0.41667 -0.011 0.003
(11) 0.45833 -0.011 0.002
(12) 0.50000 -0.011 0.000
(13) -0.45833 -0.011 -0.002
(14) -0.41667 -0.011 -0.003
(15) -0.37500 -0.011 -0.005
(16) -0.33333 -0.012 -0.007
(17) -0.29167 -0.012 -0.009
(18) -0.25000 -0.012 -0.012
(19) -0.20833 -0.012 -0.015
(20) -0.16667 -0.012 -0.021
(21) -0.12500 -0.013 -0.031
(22) -0.08333 -0.016 -0.058
(23) -0.04167 0.127 0.967

Variable: cf_ezfftf
Type: float
Total Size: 96 bytes
             24 values
Number of Dimensions: 2
Dimensions and sizes: [2] x [12]
Coordinates:
Number Of Attributes: 2
   npts : 24
   xbar : 10
(0,0) 0.1273196
(0,1) -0.01554397
(0,2) -0.01287115
(0,3) -0.0121421
(0,4) -0.01183357
(0,5) -0.01167448
(0,6) -0.01158284
(0,7) -0.01152635
(0,8) -0.01149109
(0,9) -0.0114695
(0,10) -0.01145764
(0,11) -0.005726894
(1,0) 0.9670882
(1,1) -0.05801087
(1,2) -0.03107371
(1,3) -0.02103074
(1,4) -0.0154218
(1,5) -0.01167448
(1,6) -0.008887824
(1,7) -0.006654739
(1,8) -0.004759768
(1,9) -0.003073242
(1,10) -0.001508427
(1,11) 0

Variable: amplitude
Type: float
Total Size: 96 bytes
             24 values
Number of Dimensions: 1
Dimensions and sizes: [24]
Coordinates:
Number Of Attributes: 1
   _FillValue : -999
(0) 20
(1) 0.9754332
(2) 0.06005727
(3) 0.03363394
(4) 0.0242842
(5) 0.01943876
(6) 0.01651021
(7) 0.01459985
(8) 0.01330948
(9) 0.01243787
(10) 0.01187409
(11) 0.01155651
(12) 0.01145379
(13) 0.01155651
(14) 0.01187409
(15) 0.01243787
(16) 0.01330948
(17) 0.01459985
(18) 0.01651021
(19) 0.01943876
(20) 0.0242842
(21) 0.03363394
(22) 0.06005727
(23) 0.9754332

Variable: amplitude_ezfftf
Type: float
Total Size: 48 bytes
             12 values
Number of Dimensions: 1
Dimensions and sizes: [12]
Coordinates:
Number Of Attributes: 1
   _FillValue : -999
(0) 0.9754332
(1) 0.06005727
(2) 0.03363394
(3) 0.0242842
(4) 0.01943876
(5) 0.01651021
(6) 0.01459985
(7) 0.01330948
(8) 0.01243787
(9) 0.0118741
(10) 0.01155651
(11) 0.005726894
v1n3 339%

On 2/7/2012 6:16 PM, Dennis Shea wrote:
> 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 19:38:37 2012