# Re: ezfftf

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Fri Dec 02 2011 - 08:33:34 MST

ezfftf is a direct call to the FFTPACK routine at:
http://www.netlib.org/fftpack/ezfftf.f
or, with better documentation
ftp://ftp.ucar.edu/dsl/lib/fftpack/ezfftf.f

===
As with any FFT, the issue is with how the coefficients are normalized.
For ezfftf:

var_zfft2 = 0.5*sum(z(0:nfreq-2)) + z(nfreq-1)

This is a *biased* estimate of the variance. [ie: 1/N and not 1/(N-1) ]
===
I was a bit surprised. I thought the variance would be
var_zfft2 = 0.5*sum(z)
but the core documentation indicates differently.

Attached is a test script.

On 12/2/11 2:01 AM, Wolfgang Langhans wrote:
> Hi,
>
> I am a bit confused about the outcome of the ezfftf Fourier Forward
> analysis. Should not the sum over all squared norms equal the total
> biased variance of the data series? How can I understand the different
> results obtained in the little example below? How comes that the sum
> amounts almost twice the actual variance?
>
> Thanks for your help in understanding this!
> Wolfgang
>
> ;======================================================
> low = -5.0
> high = 5.0
> con = (high - low) / 32766.0 ; 32766.0 forces a 0.0 to 1.0 range
> Z = new(100,float)
> do i=0,dimsizes(Z)-1
> Z(i) = low + con * rand()
> end do
>
> zmean = avg(Z)
> zvarbiased = 1./dimsizes(Z) * sum( (Z-zmean)^2 ) ;total biased
> variance of original data
>
> zfft = ezfftf(Z)
> z = zfft(0,:)^2 + zfft(1,:)^2
> nfreq = dimsizes(z)
> totvar = sum(z(0:nfreq-1)) ; sum over all squared norms of the
> complex Fourier transforms
>
> print("Variance: "+ zvarbiased + " Sum over spectrum: "+totvar )
> ;======================================================
> _______________________________________________
> 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 Fri Dec 2 08:33:42 2011

This archive was generated by hypermail 2.1.8 : Fri Dec 02 2011 - 16:10:18 MST