Re: ezfftf

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Fri Dec 02 2011 - 09:21:21 MST

Actually, after another test, it depends on if the
the series length is even or odd. This determine
how the last coefficient should be handled.

  if (N%2 .eq. 0) then
      var_zfft2 = 0.5*sum(z(0:nfreq-2)) + z(nfreq-1) ; even
  else
      var_zfft2 = 0.5*sum(z) ; odd
  end if

On 12/2/11 8:33 AM, Dennis Shea wrote:
> 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
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Dec 2 09:21:30 2011

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