I recently had reason to revisit the FFT math and realized that I've never written an FFT library really, so here goes:
let
signal = (0..1023).mapIt(
complex64(sin(TAU * 0.1 * float64(it))))
# abs gets us back into real space
# false for forward FFT, true for inverse (TODO flip it? separate names?)
frequencies = fft(signal, false).mapIt(abs(it))
The library implements a couple of well-known algorithms for powers-of-2, primes and combinations thereof, thus maintaining decent O(N log N) speed.
https://github.com/arnetheduck/nim-fftr
Enjoy!
This is awesome! Do you have any performance measurements?
I assume it is CPU only?
Do you have any performance measurements?
They're in the readme along with similar benches for fftw
CPU only
Yes, though it would be fun to add GPU as well and have it auto-select
Yes, though it would be fun to add GPU as well and have it auto-select
That would be pretty cool, let us know if you do it! :-)
BTW, there is a small issue in your example. It is missing an import of std/sequtils (to be able to use mapIt).
Also, do you plan to add this library to the nimble package list? I had to install it using the github URL. It might also be a good idea to add a reference to this this to the scinim package list.
Now 2x faster ;)
small issue in your example.
fixed, thanks
friendly license
"Add MIT license" button clicked to appease those that find this .. important.
# false for forward FFT, true for inverse (TODO flip it? separate names?)
i.m.h.o separate names: fft, ifft
Just playing with fftr I noticed two things,
When going back from the frequency domain to signal, I have to divide the real part of the result by the length of the signal. Is it the convention that I do that, or should the library?
It seems that depending on the length of the signal, the ifft result is flipped sign wise. When is this to be expected?
import fftr, std/[math, sequtils]
let
length = 456 #64
signal = (0..length).mapIt(complex64(sin(TAU * 0.1 * float64(it))))
frequencies = fft(signal, false)
backToSignal = fft(frequencies, true).mapIt(it.re / length.float64)
echo signal
echo "\n"
echo frequencies
echo "\n"
echo backToSignal
I'm not sure, but I would expect that division by N should included in the outcome, without your intervention. I would look at the "competitor" implementation.
Y = fft(X)
X = iftt(Y)
When going back from the frequency domain to signal, I have to divide the real part of the result by the length of the signal. Is it the convention that I do that, or should the library?
Looks like every library does its own thing here - ie scale, not scale, offer options to scale - for now, fftr does no scaling (since scaling might be needed in either direction and sometimes by the square root, if done in both directions) - this is probably a reasonable starting point though maybe confusing for ... python and matlab users. hm.. see also https://stackoverflow.com/questions/62881003/fftw-backwards-transform-is-multiplied-by-n
It seems that depending on the length of the signal, the ifft result is flipped sign wise.
Bug, should be fixed now
signal = (0..length)
btw, you have an off-by-one in normalisation due to the unfortunate design in nim that intervals are closed.. the length of this signal is 457 elements.
publish
eh I guess I'll get to it eventually but feel free to PR it yourself in the meantime ;)