# Replacing FFT from Numerical Recipes with FFTW

**URL:** <https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508>\
**Category:** Numerics\
**Tags:** question, fftw\
**Created:** [February 18, 2021, 5:34am UTC](https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508 "2021-02-18T05:34:43Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![edwinb-ai](https://avatars.discourse-cdn.com/v4/letter/e/bb73d2/32.png) [@edwinb-ai](https://discourse.julialang.org/u/edwinb-ai)\
**Post date:** [February 18, 2021, 5:34am UTC](https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508/1 "2021-02-18T05:34:43Z")

</div>

Hello there.

I’m tasked to port an old FORTRAN 77 legacy code to Julia. Unfortunately, I’m stuck with the computation of a Fast Fourier Transform (FFT).

The original legacy code uses the `sinft` routine from the [Numerical Recipes](http://www.numerical.recipes/) book (NR). The original version is in FORTRAN 77, but I managed to port it to Julia.  
As it is somewhat big, I posted it in this [GitHub Gist](https://gist.github.com/edwinb-ai/847e948814ac49ec8f7b9530b8ac9e13). I know it gives the correct answer as I have already checked it with the FORTRAN 77 code.

What I’m trying to accomplish is to successfully replace that old subroutine and use `FFTW.jl` instead.

Here is a sample code, which uses the `sinft` function from the gist I linked above

```julia
function main()
    a = [0.0,
        0.47742664734683210,
        0.75851932870459160,
        0.75733104152093889,
        0.53816912369779424,
        0.25985479516207322,
        6.5606350402976049E-002,
        8.4817932682581046E-004,
        1.6555040127766884E-002,
        3.9188171512774567E-002,
        3.4228256134543805E-002,
        1.3420820869486524E-002,
        6.5631416321955491E-004,
        2.7825880087286466E-003,
        9.4219993149478529E-003,
        9.9865908960776392E-003,
        4.6428562500441268E-003,
        3.3327448406347924E-004,
        8.9848071988801400E-004,
        3.6459412006396999E-003,
        4.2443023850689526E-003,
        2.1409386996221803E-003,
        1.8361141730897351E-004,
        3.9074461539924637E-004,
        1.7587883887050227E-003,
        2.1636202084949888E-003,
        1.1499838894481705E-003,
        1.1004675140620131E-004,
        2.0345801653869668E-004,
        9.8374907856643950E-004,
        1.2557106901589935E-003,
        6.9012160386323635E-004
    ]
    n = length(a)
    b = FFTW.r2r(a, FFTW.RODFT00)
    @show b
    c = FFTW.r2r(a, FFTW.RODFT10)
    @show c
    d = FFTW.r2r(a, FFTW.RODFT01)
    @show d
    e = FFTW.r2r(a, FFTW.RODFT11)
    @show e
    sinft!(a, n)
    @show a
end

main()

```

and the output is

```julia
fftw00 = [2.2839824978332457, 3.865808478514381, 4.727281203842373, 4.8731271633513895, 4.430388346987334, 3.52627738343957, 2.3554519677513883, 1.0916035814443859, -0.06044936548148805, -0.9482420224720198, -1.4349198347123648, -1.4922387963947052, -1.1878127966462944, -0.9024039195111244, -0.7544336573013076, -0.6487524977522947, -0.5746854983744873, -0.5109527882913876, -0.4572230988934566, -0.4071294272689463, -0.36317761272948196, -0.32225146396310084, -0.2854515631367159, -0.25067989805628743, -0.21858457535092202, -0.18781719264787639, -0.1587922173412804, -0.13063489068415524, -0.10357909252012387, -0.07705172337385213, -0.05113541033570179, -0.025469727216377414]
fftw10 = [2.0793925679832097, 3.540836745898296, 4.411873130993651, 4.69904918778722, 4.50023727641295, 3.8943771147970456, 3.0155507529728913, 1.9808566633749198, 0.9463285287804625, 0.046019457467148916, -0.5751837487126715, -0.8217911065999712, -0.6986234327172187, -0.5578142627105277, -0.4834741181879493, -0.4335746010153527, -0.4015978869433933, -0.37755252780671394, -0.35778123449159294, -0.34117680841924375, -0.326536173530114, -0.3148024248728375, -0.3041677537130595, -0.2960230382112425, -0.28837505241028905, -0.2828993208737436, -0.2775606516187846, -0.27420195090090593, -0.27077956214449317, -0.26922722005632865, -0.2675076531557324, -0.2676073339655841]
fftw01 = [1.243700805258963, 3.2418116193843214, 4.452525840891483, 4.894734592726747, 4.649826958896702, 3.8623787670102443, 2.720624387729786, 1.4246826482841162, 0.19094388443409294, -0.796962023389162, -1.3819492164469134, -1.5103369064524736, -1.2283983607363878, -0.9173515045448337, -0.7609829467521276, -0.6493873891836436, -0.5742902048686588, -0.5076348011618095, -0.45375468347680226, -0.40149913632711026, -0.35776824660920276, -0.3151620276234529, -0.27871627890030914, -0.2424743837845238, -0.21074645847268259, -0.17859150548271852, -0.1499037585025491, -0.1203668000930191, -0.09359161262773878, -0.0656279913804777, -0.039911622662318735, -0.012688562150492078]
fftw11 = [1.100731615661763, 2.888152469616345, 4.0510814187228, 4.622312215367247, 4.6568191497328115, 4.239205077660919, 3.4831742683814957, 2.506879365762653, 1.4566061868702822, 0.46800391408490905, -0.3038493851830596, -0.7540537756465204, -0.786778367704533, -0.6145829005559422, -0.5181608498286434, -0.4544324822018745, -0.4171609003337509, -0.38806300324536736, -0.3678971086135975, -0.3486464965752973, -0.3339712796773492, -0.3200728029098461, -0.30955018816292434, -0.29966980779798313, -0.29223263070962957, -0.28532961396706014, -0.2802314150004216, -0.27565684140465374, -0.27245794969244197, -0.2698415635994903, -0.2682975643179987, -0.267445642341788]
a = [0.0, 0.9043828218719977, 1.5534403130729006, 1.9790600667238083, 2.185014687850934, 2.207864522352749, 2.065766579627678, 1.7994284395325442, 1.436707753647941, 1.0275526531433679, 0.6135474148730728, 0.25895888711416704, 0.019003685271627706, -0.0361408278912162, -0.016458308980858587, -0.008379067999873493, -0.0015497022834815022, -0.000180369999170496, 0.0007441171767517737, 7.66079446568018e-5, 0.00022392036694572626, -9.59574096880151e-5, 7.835687804802483e-5, -7.522077772007148e-5, 5.0998401781382485e-5, -4.406720381025693e-5, 3.4210803845441395e-5, -2.651523625643755e-5, 2.0307333832680907e-5, -1.4499172994538248e-5, 9.41448920865895e-6, -4.621610797618872e-6]

```

Here, the first four arrays represent each of the four variations of the Fourier sine transforms as defined in the [documentation](http://fftw.org/fftw3_doc/1d-Real_002dodd-DFTs-_0028DSTs_0029.html#g_t1d-Real_002dodd-DFTs-_0028DSTs_0029) for the `FFTW`.  
And `a` is the in-place result of applying the `sinft` function from NR to the original `a` array.

The first thing that is very clear is the first element of the arrays. The `sinft` always leaves the first element as a zero, and then computes everything else. The `FFTW` on the other hand, does not.

Further, every variation of the Fourier sine transform from `FFTW` is giving me different results, as expected.

The question is: how can I reproduce the results from `sinft` using the functions in `FFTW.jl`?  
Any help would be greatly appreciated.

---

<div class="post-metadata">

**Author:** ![edwinb-ai](https://avatars.discourse-cdn.com/v4/letter/e/bb73d2/32.png) [@edwinb-ai](https://discourse.julialang.org/u/edwinb-ai)\
**Post date:** [February 19, 2021, 12:25am UTC](https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508/2 "2021-02-19T00:25:14Z")

</div>

Reading the Numerical Recipes book I came across this definition, after a lot of algebra and computations

F\_k = 2i \sum\_{j=0}^{N-1} f\_j \sin{(\pi j k / N)}

and what I ended up doing was realizing that this means that if I use the `rfft` function I should be able  
to obtain what I need.

In short, I managed to solve it using the following.

- First, pad the original array with a new array twice its size.
- Then, apply the `rfft` function.
- Finally, compute its complex conjugate and extract the imaginary part.

At last, the code I used looks something like this

```julia
function main()
    a = [0.0,
        0.47742664734683210,
        0.75851932870459160,
        0.75733104152093889,
        0.53816912369779424,
        0.25985479516207322,
        6.5606350402976049E-002,
        8.4817932682581046E-004,
        1.6555040127766884E-002,
        3.9188171512774567E-002,
        3.4228256134543805E-002,
        1.3420820869486524E-002,
        6.5631416321955491E-004,
        2.7825880087286466E-003,
        9.4219993149478529E-003,
        9.9865908960776392E-003,
        4.6428562500441268E-003,
        3.3327448406347924E-004,
        8.9848071988801400E-004,
        3.6459412006396999E-003,
        4.2443023850689526E-003,
        2.1409386996221803E-003,
        1.8361141730897351E-004,
        3.9074461539924637E-004,
        1.7587883887050227E-003,
        2.1636202084949888E-003,
        1.1499838894481705E-003,
        1.1004675140620131E-004,
        2.0345801653869668E-004,
        9.8374907856643950E-004,
        1.2557106901589935E-003,
        6.9012160386323635E-004,
    ]
    n = length(a)
    m = 2 * n
    b = zeros(m)
    b[1:n] = a # We pad with zeros
    @show rfft(b) |> conj |> imag
    sinft!(a, n)
    @show a
end

main()

```

and the result is

```julia
(rfft(b) |> conj) |> imag = [-0.0, 0.9043828218719976, 1.5534403130729006, 1.9790600667238079, 2.1850146878509342, 2.207864522352748, 2.0657665796276783, 1.799428439532543, 1.4367077536479407, 1.0275526531433665, 0.613547414873073, 0.2589588871141659, 0.019003685271627713, -0.03614082789121714, -0.01645830898085865, -0.008379067999874062, -0.0015497022834817242, -0.00018036999917095398, 0.0007441171767516019, 7.660794465619292e-5, 0.0002239203669455736, -9.595740968863682e-5, 7.835687804816383e-5, -7.522077772060466e-5, 5.0998401781576774e-5, -4.4067203810568076e-5, 3.421080384571918e-5, -2.651523625655967e-5, 2.0307333832514374e-5, -1.4499172994675291e-5, 9.414489208547927e-6, -4.621610797617137e-6, -0.0]
a = [0.0, 0.9043828218719977, 1.5534403130729006, 1.9790600667238083, 2.185014687850934, 2.207864522352749, 2.065766579627678, 1.7994284395325442, 1.436707753647941, 1.0275526531433679, 0.6135474148730728, 0.25895888711416704, 0.019003685271627706, -0.0361408278912162, -0.016458308980858587, -0.008379067999873493, -0.0015497022834815022, -0.000180369999170496, 0.0007441171767517737, 7.66079446568018e-5, 0.00022392036694572626, -9.59574096880151e-5, 7.835687804802483e-5, -7.522077772007148e-5, 5.0998401781382485e-5, -4.406720381025693e-5, 3.4210803845441395e-5, -2.651523625643755e-5, 2.0307333832680907e-5, -1.4499172994538248e-5, 9.41448920865895e-6, -4.621610797618872e-6]

```

There is a trailing zero in the transformed array using `rfft`, but I can deal with that.

It might not be the best solution, but it’s good enough for me. Although, other solutions are greatly appreciated.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 19, 2021, 4:15am UTC](https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508/3 "2021-02-19T04:15:19Z")

</div>

> [@edwinb-ai](#):
>
> Reading the Numerical Recipes book I came across this definition, after a lot of algebra and computations […]

This is a type-I DST (discrete sine transform), which is equivalent to [`RODFT00`](http://fftw.org/fftw3_doc/1d-Real_002dodd-DFTs-_0028DSTs_0029.html#g_t1d-Real_002dodd-DFTs-_0028DSTs_0029) in FFTW.

However unlike Numerical Recipes, FFTW does not store or compute with the redundant `j=0` and `k=0` elements of the input and output arrays, which are always zero. So, you just omit this when calling FFTW:

```julia
julia> FFTW.r2r(a[2:end], FFTW.RODFT00) / 2
31-element Array{Float64,1}:
  0.9043828218719976
  1.5534403130729009
  1.9790600667238076
  2.1850146878509342
  2.2078645223527475
  2.0657665796276783
  1.7994284395325428
  1.4367077536479407
  1.0275526531433667
  0.6135474148730728
  0.25895888711416587
  0.0190036852716276
 -0.0361408278912172
 -0.01645830898085865
 -0.008379067999874298
 -0.0015497022834815022
 -0.0001803699991712454
  0.000744117176751713
  7.660794465624843e-5
  0.00022392036694568462
 -9.595740968876276e-5
  7.835687804830238e-5
 -7.52207777205971e-5
  5.099840178124371e-5
 -4.406720381067153e-5
  3.421080384558017e-5
 -2.6515236256630104e-5
  2.030733383229233e-5
 -1.4499172994897336e-5
  9.414489208436905e-6
 -4.621610797561626e-6

```

(I also divided by 2 to match the normalization in your `sinft!` function. In practice, you can usually absorb this or any other normalization factor into some other part of your computation.)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 19, 2021, 4:25am UTC](https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508/4 "2021-02-19T04:25:46Z")

</div>

Note that to get the full benefit from FFTW (assuming this computation is performance critical), you should precompute the plan, and possibly preallocate the output:

```julia
julia> using FFTW, BenchmarkTools, LinearAlgebra

julia> b = a[2:end]/2; result = copy(b);

julia> p = FFTW.plan_r2r(result, FFTW.RODFT00, flags=FFTW.MEASURE)
FFTW r2r RODFT00 plan for 31-element array of Float64
(rodft00e-r2hc-pad-31
  (rdft-r2hc-direct-r2c-64 "r2cf_64")
  (rdft-rank0-iter-ci/1-x31))

julia> a2 = a*0;

julia> @btime sinft!($a2, $(length(a)));
  570.503 ns (0 allocations: 0 bytes)

julia> @btime mul!($result, $p, $b); 
  397.590 ns (0 allocations: 0 bytes)

```

Note, by the way, that the _Numerical Recipes_ `sinft` algorithm [suffers from poor accuracy](https://github.com/FFTW/fftw3/blob/17b4252e3da62490d3e206569ebfec552d09afc1/reodft/redft00e-r2hc.c#L24-L30) for long lengths (although this does not seem to be widely known—FFTPACK has the same problem); IIRC, its [rms](https://en.wikipedia.org/wiki/Root_mean_square) floating-point errors grow as O(\sqrt{n}) rather than the O(\sqrt{\log n}) achieved by good FFT algorithms. In consequence, the algorithm that FFTW’s plan `p` is using above is essentially the same as the one you implemented in terms of `rfft`, though for larger sizes it may use a different algorithm.

(That being said, REDFT00 and RODFT00 are two of the less-optimized transforms in FFTW. If you _really_ care about the speed here, there are ways to do better.)

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [February 19, 2021, 2:00pm UTC](https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508/5 "2021-02-19T14:00:01Z")

</div>

I’ve been preallocating and using

```julia
FFTW.plan_r2r!(zstore, FFTW.RODFT00, 1)

```

for several months and thought I was happy. Is there something simple I can do to improve performance.? I’m using the fast sine transform for a Poisson solver.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 19, 2021, 2:18pm UTC](https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508/6 "2021-02-19T14:18:18Z")

</div>

> [@ctkelley](#):
>
> Is there something simple I can do to improve performance.? I’m using the fast sine transform for a Poisson solver.

You could try using a [type-II DST (`RODFT10`)](https://en.wikipedia.org/wiki/Discrete_sine_transform) instead (which correponds to shifting the samples by half a grid point, so it should be easy to adapt for a spectral solver). IIRC, the algorithms in FFTW for a DST-II are a bit more efficient than those for a DST-I.

Also, be careful of the factorization of your sizes. `RODFT00` is most efficient for lengths `n` where `n+1` is highly composite, but for `RODFT11` it should be `n` that is highly composite.

`flags=FFTW.MEASURE` or `flags=FFTW.PATIENT` wouldn’t hurt to try either.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [February 19, 2021, 2:39pm UTC](https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508/7 "2021-02-19T14:39:42Z")

</div>

> [@stevengj](#):
>
> Also, be careful of the factorization of your sizes. `RODFT00` is most efficient for lengths `n` where `n+1` is highly composite, but for `RODFT11` it should be `n` that is highly composite.
> 
> `flags=FFTW.MEASURE` or `flags=FFTW.PATIENT` wouldn’t hurt to try either.

Thanks, @stevengj

I’m stuck with FT00 unless I want to do some massive recoding. n+1 = 2^m for me, which is what it was when I used the fft in matlab for decades. So I’m ok there.

I will look into the flags. I had not thought of that.

---

<div class="post-metadata">

**Author:** ![edwinb-ai](https://avatars.discourse-cdn.com/v4/letter/e/bb73d2/32.png) [@edwinb-ai](https://discourse.julialang.org/u/edwinb-ai)\
**Post date:** [February 19, 2021, 7:40pm UTC](https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508/8 "2021-02-19T19:40:05Z")

</div>

Thanks you so much @stevengj ! Learned a lot from these answers. I’ll definitely try out these solutions as well in my code.

---

<div class="post-metadata">

**Author:** ![haflquantumhalfchaos](https://avatars.discourse-cdn.com/v4/letter/h/85e7bf/32.png) [@haflquantumhalfchaos](https://discourse.julialang.org/u/haflquantumhalfchaos)\
**Post date:** [August 11, 2024, 12:03pm UTC](https://discourse.julialang.org/t/replacing-fft-from-numerical-recipes-with-fftw/55508/9 "2024-08-11T12:03:17Z")

</div>

in f77 fft code the data in function call must be twice the size of original array