# Difference in \`fft\` implementation between stan and numpy/scipy

**URL:** <https://discourse.mc-stan.org/t/difference-in-fft-implementation-between-stan-and-numpy-scipy/28642>\
**Category:** Modeling\
**Created:** [August 24, 2022, 6:06pm UTC](https://discourse.mc-stan.org/t/difference-in-fft-implementation-between-stan-and-numpy-scipy/28642 "2022-08-24T18:06:56Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![tillahoffmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/tillahoffmann/32/14838_2.png) [@tillahoffmann](https://discourse.mc-stan.org/u/tillahoffmann)\
**Post date:** [August 24, 2022, 6:06pm UTC](https://discourse.mc-stan.org/t/difference-in-fft-implementation-between-stan-and-numpy-scipy/28642/1 "2022-08-24T18:06:56Z")

</div>

Thank you for the Fourier implementation in the newest release! Stan [defines the `fft`](https://github.com/stan-dev/docs/blob/1961f88828e406ff31a85dff699b67f686efd7dc/src/functions-reference/complex_matrix_operations.Rmd#L1204-L1213) u\in\mathbb{C}^N of a complex vector v\in\mathbb{C}^N as

u\_n = \sum\_{m \< n} v\_m \cdot \exp\left(\frac{-n \cdot m \cdot 2 \cdot \pi \cdot \sqrt{-1}}{N}\right).

Numpy instead [defines it as](https://numpy.org/doc/stable/reference/routines.fft.html#implementation-details)

u\_n=\sum\_{m=0}^{n-1} v\_m\exp\left(\frac{-2\pi i mn}{N}\right).

These are identical up to the difference in summation (over all indices for numpy and only over “time” indices less than “frequency” indices in Stan). Of course, not everything needs to follow the same conventions as numpy, but I just wanted to understand the summation convention in Stan because it differs from previous definitions I’ve seen, e.g. [on wikipedia](https://en.wikipedia.org/wiki/Fast_Fourier_transform#Definition).

The results also seem to differ when comparing test data.

```stan
// test_fft.stan
data {
    int n;
    vector[n] x;
}

generated quantities {
   complex_vector[n] y = fft(x);
}

```

```python
# test_fft.py
import cmdstanpy
import numpy as np

np.random.seed(0)
model = cmdstanpy.CmdStanModel(stan_file="test_fft.stan")
n = 3
x = np.random.normal(0, 1, n)
fit = model.sample({"n": n, "x": x}, fixed_param=True, iter_warmup=0, iter_sampling=1)
stan_fft, = fit.stan_variable("y")
np_fft = np.fft.fft(x)
for label, fft in [("stan", stan_fft), ("np", np_fft)]:
    print(label)
    print("real", np.round(fft.real, 3))
    print("imag", np.round(fft.imag, 3))
    print()

cmdstanpy.show_versions()

```

```nohighlight
# output
stan
real [3.143 0. 1.075]
imag [0.501 1.075 -0.501]

np
real [3.143 1.075 1.075]
imag [0. 0.501 -0.501]

INSTALLED VERSIONS
---------------------
python: 3.10.0 (default, Jan 15 2022, 13:32:36) [Clang 12.0.5 (clang-1205.0.22.9)]
python-bits: 64
OS: Darwin
OS-release: 21.6.0
machine: arm64
processor: arm
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: ('en_US', 'UTF-8')
cmdstan_folder: /Users/till/.cmdstan/cmdstan-2.30.1
cmdstan: (2, 30)
cmdstanpy: 1.0.4
pandas: 1.4.3
xarray: None
tdqm: None
numpy: 1.22.1
ujson: 5.4.0

```

Do you know why the values differ? As the signal is real, I would’ve expected the zero-frequency term to have zero imaginary part (as numpy does).

---

<div class="post-metadata">

**Author:** ![WardBrian](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wardbrian/32/12078_2.png) [@WardBrian](https://discourse.mc-stan.org/u/WardBrian)\
**Post date:** [August 24, 2022, 6:55pm UTC](https://discourse.mc-stan.org/t/difference-in-fft-implementation-between-stan-and-numpy-scipy/28642/2 "2022-08-24T18:55:10Z")

</div>

Stan wraps the Eigen C++ library for the implementation of FFTs. They have a few notes on different choices they’ve made [here](https://eigen.tuxfamily.org/index.php?title=EigenFFT#Behavior_Differences_vs_other_FFT_libraries).  
I’m not sure if this accounts for the difference you’ve noted or not

---

<div class="post-metadata">

**Author:** ![tillahoffmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/tillahoffmann/32/14838_2.png) [@tillahoffmann](https://discourse.mc-stan.org/u/tillahoffmann)\
**Post date:** [August 24, 2022, 7:16pm UTC](https://discourse.mc-stan.org/t/difference-in-fft-implementation-between-stan-and-numpy-scipy/28642/3 "2022-08-24T19:16:40Z")

</div>

Thank you for the link. They mention scaling which should not affect the results up to a global rescaling. The other difference is how they handle real FFTs where only half the spectrum is required. Neither should lead to the differences observed above, I think.

---

<div class="post-metadata">

**Author:** ![WardBrian](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wardbrian/32/12078_2.png) [@WardBrian](https://discourse.mc-stan.org/u/WardBrian)\
**Post date:** [August 24, 2022, 7:39pm UTC](https://discourse.mc-stan.org/t/difference-in-fft-implementation-between-stan-and-numpy-scipy/28642/4 "2022-08-24T19:39:12Z")

</div>

Oh my, it seems this is a bug not in Stan, but in cmdstanpy:

```stan
generated quantities {
   complex_vector[n] y = fft(x);
   print(get_real(y));
   print(get_imag(y));
}

```

outputs

```nohighlight
Chain [1] [3.14295,1.0746,1.0746]
Chain [1] [0,0.501066,-0.501066]

```

This is not good! Thank you for finding this

---

<div class="post-metadata">

**Author:** ![tillahoffmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/tillahoffmann/32/14838_2.png) [@tillahoffmann](https://discourse.mc-stan.org/u/tillahoffmann)\
**Post date:** [August 24, 2022, 8:08pm UTC](https://discourse.mc-stan.org/t/difference-in-fft-implementation-between-stan-and-numpy-scipy/28642/5 "2022-08-24T20:08:59Z")

</div>

Aha, nice catch! Any idea why cmdstanpy is returning the odd values?

I’ll extract the real and imaginary parts in separate vectors for the time being.

---

<div class="post-metadata">

**Author:** ![WardBrian](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wardbrian/32/12078_2.png) [@WardBrian](https://discourse.mc-stan.org/u/WardBrian)\
**Post date:** [August 24, 2022, 8:22pm UTC](https://discourse.mc-stan.org/t/difference-in-fft-implementation-between-stan-and-numpy-scipy/28642/6 "2022-08-24T20:22:20Z")

</div>

There isn’t a better answer than “we didn’t check our work here” sadly. The CmdStanPy output assumes that complex numbers in the output are stored like an extra final dimension of size 2. But they’re not really, they’re stored strided alongside the other dimensions.

To make this concrete, a vector of 3 elements was being read in as if it was stored as `x.1.real, x.2.real, x.3.real, x.1.imag. x.2.imag, x.3.imag`, but really it is actually written to disk as `x.1.real, x.1.imag, x.2.real, x.2.imag, ...`

We should have tested this, but we didn’t. Our tests only asserted the result had the right shape, not the right order, which I’ll rectify as I fix the bug. My guess is people haven’t used complex outputs very much yet in Stan, so nobody noticed it was broken.

---

<div class="post-metadata">

**Author:** ![tillahoffmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/tillahoffmann/32/14838_2.png) [@tillahoffmann](https://discourse.mc-stan.org/u/tillahoffmann)\
**Post date:** [August 24, 2022, 8:38pm UTC](https://discourse.mc-stan.org/t/difference-in-fft-implementation-between-stan-and-numpy-scipy/28642/7 "2022-08-24T20:38:03Z")

</div>

Thank you for the additional context and the great work with CmdStanPy!

---

<div class="post-metadata">

**Author:** ![WardBrian](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wardbrian/32/12078_2.png) [@WardBrian](https://discourse.mc-stan.org/u/WardBrian)\
**Post date:** [August 24, 2022, 8:49pm UTC](https://discourse.mc-stan.org/t/difference-in-fft-implementation-between-stan-and-numpy-scipy/28642/8 "2022-08-24T20:49:53Z")

</div>

Thank you again for being curious and finding the bug!

I suppose getting the same output now still doesn’t answer your question of why the definition in the docs is nonstandard. For that I’ll ping @Bob_Carpenter, who wrote the docs