normalization in FFTW3

Question

I am trying to use FFTw3 in Visual Studio using the Intel Fortran compiler. I am confused about normalization. I have setup a r2r spectral to grid plan for sinx siny components.

 plan_spectral_to_grid_sin_sin = fftwf_plan_r2r_2d(Ny, Nx, dummy, dummy, FFTW_RODFT00, FFTW_RODFT00, FFTW_MEASURE)
 plan_grid_to_spectral_sin_sin = fftwf_plan_r2r_2d(Ny, Nx, dummy, dummy, FFTW_RODFT00, FFTW_RODFT00, FFTW_MEASURE)

(Note the reversal of the x, y dimensions since this is a C library and FORTRAN is column-major)

I understand (maybe this is incorrect) that the normalization factor is

sin_sin_norm = (2.0 * REAL(currentDims%IXMAX + 1) * 2.0 * REAL(currentDims%IYMAX + 1))

So if I divide my spectral components A(k,l)sinkxsinly by this value and do the transform

call fftwf_execute_r2r(plan_spectral_to_grid_sin_sin, horzArr%G1, horzArr%G1)

I should be able to do the transform back

call fftwf_execute_r2r(plan_grid_to_spectral_sin_sin, horzArr%G1, horzArr%G1)

and get the original spectral components WITHOUT normalization. But that is not what I am seeing.

For reference

horzArr%G1 is a 2-d array of dimension currentDims%IXMAX, currentDims%IYMAX

• Martin Brown: I can't recall which normalisation convention FFTW3 uses for that r2r transform. If you take it on a round trip without trying to normalise it at all what factor does it get multiplied by? Fence post errors are all too easy in C/FORTRAN hybrid interface code.
• Vladimir F Героям слава: What are you seeing? Is your result simply different by a multiplicative factor? Can you compute the factor? The normalisation in FFTW always requires the logical transform size.
• Vladimir F Героям слава: Note that I have no idea what currentDims%IXMAX is and how it is related to Nx. Try to show a minimal reproducible example.
• Brian Reinhold: @VladimirFГероямслава I guess I should have put horzArr%G1 is a 2-d array of dimension currentDims%IXMAX, currentDims%IYMAX. I was careless.
• Brian Reinhold: @MartinBrown Good idea and that's exactly how I am trying to approach it. I am doing a spectral meteorological model and there one typically does spectral to grid, compute nonlinear terms, and then grid to spectral. Since we work in spectral space the 'inverse' transform is the grid to spectral which is opposite of the documentation definition. That confused me for a while. Given that, I have lined up some code that does spectral to grid and then grid to spectral one after the other. That way I can find out what the normalization is with the help of the debugger.

Answer 1

Following the comment made by Martin Brown using the debugger and some modest code changes I found that taking

sin_sin_norm = (2.0 * REAL(currentDims%IXMAX + 1) * 2.0 * REAL(currentDims%IYMAX + 1))

I could set a spectral component horzArr%G1(2,2) = 1.0 / sin_sin_norm and do

call fftwf_execute_r2r(plan_spectral_to_grid_sin_sin, horzArr%G1, horzArr%G1)

followed by

call fftwf_execute_r2r(plan_grid_to_spectral_sin_sin, horzArr%G1, horzArr%G1)

and I would get my spectral component back. The AI I was using stated otherwise; normalization ahead of time was not sufficient. AI has been a great help (especially with syntaxes as you go from language to language and library to library) but one still has to be careful!