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

```
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

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