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
            mask = 1:frex.n_nonnegative
            
            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 :

https://docs.juliadsp.org/stable/periodograms/#DSP.Periodograms.periodogram-Union{Tuple{AbstractVector{T}},%20Tuple{T}}%20where%20T<:Number

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.

15 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

6 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 :