Re: spectral energy from 2D FFT

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.shtml

"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

      M = 2
      L = 6 ; fastest varying dimension
      x = new ( (/M,L/), "float")
      x(0,0) = 0.968
      x(0,1) = 0.067
      x(0,2) = 0.478
      x(0,3) = 0.910
      x(0,4) = 0.352
      x(0,5) = 0.933

      x(1,0) = 0.654
      x(1,1) = 0.021
      x(1,2) = 0.512
      x(1,3) = 0.202
      x(1,4) = 0.940
      x(1,5) = 0.204

      N = M*L
                                     ; [re/im] [M] [L/2+1]
      coef = fft2df (x) ; [2] x [2] x [4]
      xNew = fft2db(coef)
      print("---")
      print("x="+x+" xNew="+xNew)

      epsx = 0.00001 ; make sure transforms work
      if (any(xNew.lt.(x-epsx) .or. xNew.gt.(x+epsx))) then
          print("fft2df <==> fft2db: something wrong")
          exit
      end if

      print("---")
      print("coef: "+coef)

      do m=0,M-1
         print("---> m="+m)
         print("coef(real): "+coef(0,m,:)+" coef(imag)="+coef(1,m,:))
      end do

      avgx = avg(x) ; overall average
      print("---")
      print("avg(x)="+avgx+" coef(0,0,0)="+coef(0,0,0))

      varxs = variance(x) ; overall sample (unbiased) variance
      varxp = (variance(x)*(N-1))/N ; overall population variance

      coef2 = coef(:,:,1:)^2 ; *ignore* the (:,:,0) coef
      varc2 = sum(coef2) ; variance
      print("---")
      print("varc2="+varc2+" varxp="+varxp+" varxs="+varxs)

             ; real imag
      varcri = sum(coef(0,:,1:)^2 + coef(1,:,1:)^2)
      print("---")
      print("varcri="+varcri)

      print("---")
      print("varc2/sqrt(2)="+varc2/sqrt(2))

      print("---")
      print("varc2/sqrt(pi)="+varc2/sqrt(3.14159))

      print("---")
      print("varc2/sqrt(pi/2)="+varc2/sqrt(3.14159/2))

      print("---")
      print("varc2/(pi/2)="+varc2/(3.14159/2))

      print("---")
      print("varc2/(sqrt(pi)/2)="+varc2/(sqrt(3.14159)/2))

                                       ; partition by component
      avgM = dim_avg(x) ; avg for each M
      varM = (dim_variance(x)*(L-1))/L ; population variance for each M
      print("---")
      print("avgM="+avgM+" varM="+varM)

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Feb 12 10:54:02 2010