[ncl-talk] some clarifications on fft2df function

[Dennis Shea]

Buongiorno Guido,

I think what you are doing is correct.

Unfortunately,  I think I wrote the "pretty limited" documentation ... :-(
Partly, it is "pretty limited" because I don't use 'fft2d'.
-------------------------------------------------------------------

"I believe that fft2df is returning the real and imaginary part of the
Fourier transform"

Yes.

---

I believe NCL is using:
   https://www2.cisl.ucar.edu/resources/legacy/fft5
   https://www2.cisl.ucar.edu/resources/legacy/fft5/documentation

NCL's fft2{f/b} are interfaces to

   - RFFT2I
   <https://www2.cisl.ucar.edu/resources/legacy/fft5/documentation#rfft2i.html>
   2D real initialization   <=== users do not see this. The interface call this
   - RFFT2B
   <https://www2.cisl.ucar.edu/resources/legacy/fft5/documentation#rfft2b.html>
   2D real backward     <===
   - RFFT2F
   <https://www2.cisl.ucar.edu/resources/legacy/fft5/documentation#rfft2f.html>
   2D real forward         <===

------
The following interface software clearly shows that it is doing nothing
more than calling the fortran subroutines and explicitly extracting the
real and imaginary coefficients. (This latter step is necessary because NCL
does not support complex variables.)

-----
[SNIP from NCL's  C-Fortran interface.]

/*
 * Call the Fortran routines.
 */
  ier = 0;
  NGCALLF(drfft2i,DRFFT2I)(&il, &im, wsave, &ilwsave, &ier);
  NGCALLF(drfft2f,DRFFT2F)(&ildim, &il, &im, tmp_r, wsave, &ilwsave, work,
&ilwork,
                             &ier);
  if(ier) {
    NhlPError(NhlFATAL,NhlEUNKNOWN,"fft2df: ier = %d", ier);
    return(NhlFATAL);
  }
/*
 * Copy tmp_r back to the appropriate locations in coef.
 */
  if(type_coef == NCL_float) {
    for(i = 0; i < m; i++ ) {
      for(j = 0; j < l21; j++ ) {
        ic0 = i*l21 + j;
        ic1 = ml21 + ic0;
        ir0 = i*ldim + 2*j;
        ir1 = ir0 + 1;
        ((float*)coef)[ic0] = (float)tmp_r[ir0];
        ((float*)coef)[ic1] = (float)tmp_r[ir1];
      }
    }
  }
  else {
    for(i = 0; i < m; i++ ) {
      for(j = 0; j < l21; j++ ) {
        ic0 = i*l21 + j;
        ic1 = ml21 + ic0;
        ir0 = i*ldim + 2*j;
        ir1 = ir0 + 1;
        ((double*)coef)[ic0] = tmp_r[ir0];
        ((double*)coef)[ic1] = tmp_r[ir1];
      }
    }
  }

[SNIP]
+++++++++++++++++++++++
+++++++++++++++++++++++

A pic I got from 'somehere'. Sorry, I don't know the source.

http://www.cgd.ucar.edu/~shea/fft2d_image_ncl

Hope all this verbal meandering is useful!

Cheers
D

On Wed, Mar 8, 2017 at 8:24 AM, Guido Cioni <guidocioni at gmail.com> wrote:

> Hi all,
> I'm trying to perform a spectral analysis of surface fluxes (mainly
> sensible heat) in order to identify spatial scales of surface heterogeneity
> [Baidya Roy et al., 2003]. In order to test the method I'm first using
> idealized simulations where the heterogeneity is prescribed with one dry
> and one wet patch aligned in the x-dimension with very different fluxes.
> What I'd expect to find in this case is a peak at the scale of ~100km,
> which is the size of the patch.
>
> I've been trying to use the fft2df function, but the documentation is
> pretty limited, so that's why I'm asking here. The idea is to go from
> sensible heat flux, shfl_s(time, lat, lon), to the Fourier transform
> F(lhfl_s)(time, k_lat, k_lon) and then somehow collapse the 2-D DFT into a
> one-dimensional by averaging. Since the idealized domain is limited but
> with periodic boundary conditions I don't think that tapering would be
> necessary.
>
> So I first apply
>
>   coef:=fft2df(shfl_s(:,:))
>   coefmod:=sqrt(coef(0,:,:)^2+coef(1,:,:)^2)
>   spectrum=dim_avg_n_Wrap(coefmod, 0)
>
> I believe that fft2df is returning the real and imaginary part of the
> Fourier transform, so that I have to compute the absolute value and then an
> average to get rid of the y-coordinate, which is useless given the
> heterogeneity gradient in my case is everywhere parallel to the x-axis. Is
> the procedure correct?
> I get a spectrum which is something like that (don't mind the x-axis...),
> but I would have expected something different...
>
>
>
> Guido Cioni
> http://guidocioni.altervista.org
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170308/d9e914b4/attachment.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screen Shot 2017-03-08 at 16.23.22.png
Type: image/png
Size: 51828 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170308/d9e914b4/attachment.png