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:
