Re: spectral energy from 2D FFT

From: Bjoern Maronga <maronga_at_nyahnyahspammersnyahnyah>
Date: Wed Feb 17 2010 - 02:02:54 MST

Hi,

thanks for the reply.

I was playing aroung with fft2df for a while for calculating spectra. From the
FFT I got a (N x N/2+1) shaped array. Then I calculated the wave numbers by k
= sqrt(i^2 + j^2), where i=1/N to Nyquist frequency and j=1/N to Nyquist
frequency. The result is a wavenumber for each value in the array. After that
I defined a 1d array with wavenumbers k1d=1/N to Nyquist frequency.

And now comes the part which might be called into question. I calculated the
spectral energy for each k1d by averaging over all values with this
wavenumber in the 2D array, or more precisely I averaged over "similar"
k-values ( where k >= k1d(x) and k < k1d(x+1) ).

However, the results look quite good in my case, but I could not manage to
find the same variance as in the original data. I normalized the energy with
the data's variance, so that the area under the curve finally equals the
total variance.

But still I am wondering why fft2df has a (N x N/2+1) shape. If I'm right, it
should be (N/2+1 x N/2+1), because from N/2+2 onwards, all wave numbers will
be larger than Nyquist frequency (= 2*sampling rate) and therefore have to be
cut off.

Corrections or supporting comments on the thoughts above are of course
appreciated.

> NCL uses:
>
> fft2df:
> http://www.cisl.ucar.edu/css/software/fftpack5/rfft2f.html
>
> fft2db:
> http://www.cisl.ucar.edu/css/software/fftpack5/rfft2b.html
>
> ========
>
> I wrote the attached test script.
> The function pair fft2df <==> fft2db work in the sense
> that the original array can be reconstructed from the derived coefficients.
>
> However, I also can not figure out how to get the variance from
> these coefficients. The coefficients are normalized. Sometimes
> a sqrt(2), sqrt(2/pi), sqrt(pi), etc or some such coefficient
> can be used but I just can not seem to determine anything that works.
>
> Attached is a simple test script.
>
> =========
> Ideas anyone???
>
>
>
> =======
> Note:
> http://www.ncl.ucar.edu/Document/Functions/Contributed/wave_number_spc.shtm
>l
>
> "wave_number_spc" requires a global grid.
> It uses spherical harmonics to (hopefully) get the wave number spectrum.
>
> You can see the code by:
>
> %> less contributed.ncl
>
> Then search for wave_number_spc
>
> ====
>
> An example that uses fft2df
>
> http://www.ncl.ucar.edu/Applications/twoPointCor.shtml
>
> =====
>
> Good luck
>
> Bjoern Maronga wrote:
> > I have problems carrying out proper spectra from a discrete
> > two-dimensional field A (e.g. vertical velocity). Unfortunately, NCL does
> > not support the function ezfft2f yet. So I was trying to use fft2df for a
> > 2D-FFT. Following literature (e.g. Stull), the variance of my field A
> > should be equal to the sum of square of the norm of the Fourier transform
> > F_A:
> >
> > variance(A) = sum from i=1 to N-1 of | F_A(i) |^2
> > where | F_A(i) |^2 = (F_Arealpart(i)^2 + F_Aimag.part(i)^2
> >
> > for a 1D series. I tried to use this to check whether the FFT is working
> > properly. But I am confused. The equation should adapted to 2D be:
> >
> > variance(A) = sum from i=1 to N-1 sum from j=1 to N-1 of | F_A(i,j) |^2
> > where | F_A(i,j) |^2 = (F_Arealpart(i,j)^2 + F_Aimag.part(i,j)^2
> >
> > But the results are quite different (0.265 and 0.48). The function fft2df
> > transforms the array A(M x N) to an array of size (M x N/2+1). Now I'm
> > not sure, how to sum up to get correct variances. I tried summing up the
> > full array as well as only (N/2+1,N/2+1) and cared for the mean values in
> > i=0 and j=0.
> >
> > Furthermore I want to calculate the variance (energy) spectrum from the
> > Fourier transform F_A and plot variance against wave number, in
> > consideration of the grid spacing of 100m in the original data.
> >
> > I haven't found any useful information about how to handle the FFT result
> > of a 2D field and how to carry out the desired spectrum.
> >
> > Has anybody a solution or hints?
> >
> > Best regards,
> > Bjoern Maronga
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Feb 17 02:03:08 2010

This archive was generated by hypermail 2.1.8 : Tue Feb 23 2010 - 08:26:41 MST