# Why are the return-values of fft and fftfreq so terrible?

## The situation

When I calculate the fft of a signal, I want to plot it almost every time. However, there is a number of operations that I want to do each time, as the output from `fft` and `fftfreq` is honestly terrible.

So what is so bad about them? Several things.

## Problems with fft and fftfreq

First of all, the amplitude needs to be corrected for the length of the signal.

Secondly, there is no kwarg to get a onesided spectrum, which one often wants. And if you start trying to select the indices yourself, you will see the following horror:

``````julia> fftfreq(5)
5-element AbstractFFTs.Frequencies{Float64}:
0.0
0.2
0.4
-0.4
-0.2
``````

Why are the indices stored in this way? It makes no sense, and results in ugly lines if you plot with a line between each point, as they are not ordered. Luckily, there is the `fftshift` function to the rescue, which sorts the frequencies properly. Unfortunatly, it does not save us from having a signal of even length:

``````julia> fftfreq(6)
6-element AbstractFFTs.Frequencies{Float64}:
0.0
0.16666666666666666
0.3333333333333333
-0.5
-0.3333333333333333
-0.16666666666666666
``````

The largest frequency is negative - nonsense! If we plot a onesided spectrum, we almost always use the positive frequencies, so why do the negative ones get an extra valid datapoint?

To get around all of this, I defined the following function to do this processing for me:

``````
"""
fftFAP(s::AbstractArray; fs=1, onesided=true) = freqs, amps, phases

Convenience-function for computing the frequencies,
amplitudes and phases of the FFT of a signal, correcting
the amplitudes for length, sorting the frequencies
(if twosided), and including the maximal frequency
(if onsided, normally negative.)
"""
function fftFAP(s::AbstractArray; fs=1, onesided=true)
N = length(s)
frex = fftfreq(N, fs)

s_fft = fft(s)
amps = abs.(s_fft) ./ N
phases = angle.(s_fft)

if onesided
amps[begin+1:end] *= 2  #amplitude-correct for onesided
if N |> isodd  # 0 in middle, even number of elements to each side

frex, amps, phases = getindex.((frex, amps, phases), Ref(mask))
else  # for even length s, frex has the form  [0, 0.25, -0.5, -0.25]
frex, amps, phases = circshift.((frex, amps, phases), Ref(-1))  # frex = [0.25, -0.5, -0.25, 0]
mask = range(N, step=-1, length=NĂˇ2+1)  # +1 as we include a negative freq, the largest freq
frex, amps, phases = getindex.((frex, amps, phases), Ref(mask))
frex = abs.(frex)  # It is more conventional to plot positive freqs if onesided.
end
else
frex, amps, phases = fftshift.((frex, amps, phases))  #! Allocating twice by doing this
end

frex, amps, phases
end
``````

My issue is that I believe that most users want all of this processing, most of the time. Which begs the question:

### Why is there no function in the FFT packages that does this processing?

I much prefer the `fft` function from `FFTW.jl` to be consistent with the definition of the FFT. Any subsequent data massaging is better left to the users.

If you are calculating the one sided spectra of real signals, then Iâd suggest using `rfft` instead of `fft`.

2 Likes

I usually use DSP.jl for that kind of things, makes things much easier :

Agreed, consistency is a must. I am arguing that the data massaging function should ship with the fft function, with the package. And the dosctring should refer normal users to the assisted version.

That is a great point, and something i just discovered myself. But the AbstractFFTs doctrings should then refer users to use periodogram rather than FFTW, right? It just seems like it takes too much fiddling/looking to make/find such good solutions, and that they should be advertised from the start. So it is largely a documentation problem I guess.

A periodogram is not an FFT (which computes a discrete Fourier transform / DFT). Youâre making pretty strong assumptions of what a ânormalâ FFT user is if you are assuming that most people want a periodogram rather than a DFT.

e.g. some of us want to solve PDEs, or compute Chebyshev polynomials, or perform convolutions, or do other tasks that rest on a DFT as a foundation. An FFT packageâs only job is to provide that foundation.

17 Likes

Fair point. I have been quite biased in my basis for a ânormal userâ.

My personal context is as an engineering student, of which there are a few. My best guess is that all engineers learn about FFTâs at some point. A key part is hands-on calculation and visualisation. So can we agree that at minimum, a significant portion of users will a significant part of the time have similar issues with the return-valjes from âfftâ and âfftfreqâ as me?

Also, my personal feeling that one of the best things about programming is the ability to automate tedious and repeated tasks for any number of people and for a long time. I just feel like not providing an improved version of the function I defined in the OP is such a missed opportunity.

And if the FFTW package does not provide such a function, who should?

This is exactly what I need, except for the scaling by length. Thanks!

1 Like

I realize that I am talking to one of the creators of FFTW, and I want to be clear:
I am grateful for the great work put into making an amazing library, faster than anything else. I am also sure that there are good algorithmic reasons why the return-values are ordered as they are.

I am however simply left feeling like the most used FFT library for julia could be more âjulianâ in itâs API, e.g. by dispatching on the element type to use rfft automatically. In this case, the output could be onesided by default. Also, a method could be added with a second positional argument `fs`, which would return the fft values and frequencies as a named tuple, which would make `fft(s, fs)` equivalent to `fft(s), fftfreq(length(s), fs)` - nothing different in terms of function, but some simple sugar for the end user. Again, the work I am critizizing is fantastic work. Please read it in the spirit of the greed of julia (from the why we created julia blogpost).

I know I should simply âvote with my feetâ and make such a âmore julianâ (by my estimate) library myself, building on top of the amazing tooling provided by people such as yourself. I might get around to it at some point. If I did it today, it would be called âEasyFFTsâ - please lit me know if there are better potential names.

For what itâs worth, I had a go at some point. But instead of âvalues and frequencies as a named tupleâ they live together in a fancy array:

``````using AxisKeys, FFTW
wrapdims([1,im,-1,-im, 1,im,-1,-im] * 100, time=1:8) |> fft
fftshift(ans)
wrapdims([1,-1,1,-1] * 100, time=0:0.1:0.3) |> rfft
``````

`rfft` is a different transformation than `fft`. Think of `fft(x)` as `F*x` where `F` is a matrix. `rfft` is a different matrix (a submatrix of `F`). Because of that, I donât think it would be appropriate to silently change `fft(x)` to compute `rfft(x)` depending on the type of `x`.

I agree that FFTs can be tricky to use, but thatâs mostly because people have a hard time understanding the mathematical operation it represents. (For example, its relationship to the continuous Fourier transform is a perennial source of confusionâI canât tell you how many times someone has complained that the FFT of a Gaussian pulse is not exactly a Gaussian, or been confused about how to take a derivative.) This kind of confusion seems to have been endemic to FFTs ever since they were popularized in the 1960s. For example, the Numerical Recipes authors write:

Fourier routines: persistent reports of bugs, but NOT ONE has proved to be correct. Usually turns out that the bug-reporter hasnât read the book and doesnât understand the ordering of frequency components

8 Likes

Right. We are coming from very different places here, where you prioritize the appropriateness whereas I am primarily concerned with ease-of-use. I simply want to hide all the complexity into convenience-functions, making the uninformed user never learn the details, but making it easier to use. It is incredibly important to have your perspective in the discussion, whereas my perspective is a ânice-to-haveâ rather than a âmust-haveâ. Which is why I feel like something to the effect of EasyFFTs, with a name implying that appropriateness of change of return values is not of top priority, could be a good solution.

Youâre more than welcome to attempt this. However, I suspect that as soon as people get past the âletâs plot our first picture of the spectrumâ stage they will find that they have to deal with what an FFT actually computes, and youâll find that your package doesnât help so much unless it targets a very narrow range of applications.

Even if you target a narrow application, youâll then run into a lot of complications from the application itself. For example, if you want to target only plotting periodograms, youâll very quickly run into complications arising from windowing / spectral leakage, averaging of noise, etcetera. (Which would still be useful as a package, though DSP.jl already has some of this!)

3 Likes

This seems to go in your direction too :

@jonathanBieler FourierTools.jl lists the following âmain featuresâ:

• sinc interpolation allows to up and downsample a (bandlimited) signal
• FFT based convolutions
• array/image rotation
• array/image shifting (including noteworthy subpixel shifts)
• array/image shearing
• several tools like `ffts`, `ft`, `fftshift_view` etc. allowing simpler use with Fourier transforms

None of this has anything to do with deviating from the actual DFT to make life simpler for end users, which is my goal. So as far as I can tell, the effors are completely orthogonal.

I have finally gotten my `EasyFFTs.jl` enough into shape to register it:

I would be very happy if someone would test it by adding the github repo manually

``````using Pkg