Google Highway based FOSS FFT library

[zsliu98]

Recently I have decided to swtich my FFT/SIMD library to Google Highway in order to obtain the dark power of raw/semi-raw SIMD intrinsics. However, there is no FFT library built with Google Highway. Therefore, I decide to write one myself, mainly following ideas from OTFFT.

If you want to use it, see GitHub - ZL-Audio/zldsp_fft: A FFT library based on Google Highway · GitHub , it is under Apache-2.0 license. For now, it supports double/float C2C/R2C power-of-two FFT, targeting SSE2/SSE4/AVX2/NEON. AVX512/SVE are not supported as I do not have machines to develop/test

If you want to benchmark or contribute (especially regarding performance on x86-64), see GitHub - ZL-Audio/zldsp_fft_develop: FFT implementation and analysis · GitHub , it is also under Apache-2.0 license (EXCEPT FOR other FFT libraries in this repo!).

___

Now, the question comes down to the performance. Here are the throughput of real-to-complex transform (higher the better).

Disclaimer: I have tried my best to make the benchmark fair. However, there might be mistakes regarding the set-up of other FFT libraries. There might also be fluctuations in hardware, though they should get covered by the benchmark framework.

Feel free to suggest any ideas perhaps it is the best I can get with Stockham?

[zsliu98]

Some optimizations for medium/high orders And update KFR library to the latest release. Well done, KFR!

Benchmark - High-performance FFT Library zldsp::fft

[ttg]

For JUCE users, it’ll be useful if you’ll be able to make the JUCE team add it as a backend to existing juce::FFT. That way it’s a drop-in that can be easily tested within JUCE apps/plug-ins.

Another thing I was wondering is:

To activate SIMD, you must set the correct compilation flags. Highway will automatically detect the target architecture from these flags.

Especially in Intel/AMD. From my understanding of Highway, their own magic is that they have a binary part compiled with the right flag for each SIMD and they can choose it at runtime (unless explicitly set to compile time). This feels counterintuitive to the Highway design or am I missing something? (while ideal and not always possible, our goal is to have best performance on each device, so if I have AVX256 it’ll dynamically go there while same binary would still be able to branch to SSE3 on older device without crashing).

[zsliu98]

Thanks for responding.

@ttg

For JUCE users, it’ll be useful if you’ll be able to make the JUCE team add it as a backend to existing juce::FFT. That way it’s a drop-in that can be easily tested within JUCE apps/plug-ins.

I do try to design the API to make it similar to juce::dsp::FFT (for example, it has an optimized version of get RFFT squared magnitude response). However, I am not very sure how to integrate it, and I do think it needs to be approved by JUCE team first. Some things that might need to be considered:

• it needs to add Highway to CMake and some proper flags
• juce::dsp::FFT current APIs do not support SoA output while SoA is much simpler for almost any spectrum processing

@ttg

Especially in Intel/AMD. From my understanding of Highway, their own magic is that they have a binary part compiled with the right flag for each SIMD and they can choose it at runtime (unless explicitly set to compile time). This feels counterintuitive to the Highway design or am I missing something?

Yes, that is the dynamic dispatch feature of Google Highway. I use the static dispatch in this case because writing dynamic dispatch code is quite painful (you need to construct most things in stateless C-alike instead of C++). And writing SIMD width in constexpr gives me more confidence about compilter optimization (dynamic dispatch may be able to keep the same performance, though). I will consider that in the future

[lorcan]

Thanks, I’m getting great performance with your library on Apple Silicon.
On MSVC using SSE2 however, it’s rather poor unfortunately. It would be great to use dynamic dispatch here, as imposing AVX2 as a requirement is usually not realistic.

Also may I ask why FFT inputs do not take const pointers, as the data is not modified in-place as it seems ?

[zsliu98]

Thanks for trying the library

@lorcan

On MSVC using SSE2 however, it’s rather poor unfortunately.

It cannot out-perform IPP/KFR on x86-64 but should be on par with and likely faster than other libraries. Besdies that, I would guess IPP/KFR are dispatching to AVX2 while zldsp::fft is forced on SSE2, otherwise you shall not see a big gap.

It may get even faster if you use SoA for complex numbers (especially around medium orders). And SoA would also benifit spectrum processing a lot.

@lorcan

It would be great to use dynamic dispatch here, as imposing AVX2 as a requirement is usually not realistic.

Techinically it is possible. And technically it is possible to dispatch the whole processBlock() function. It is just so painful to arrange all things in that way … A much easier solution is to compile two versions and detect the SIMD via app/plugin installer, which should also avoid any runtime dispatch overhead.

@lorcan

Also may I ask why FFT inputs do not take const pointers, as the data is not modified in-place as it seems

Given that input/output might be pointed to the same address, I am not sure whether it is safe to do so, if the compiler does some optimization there?

[lorcan]

It cannot out-perform IPP/KFR on x86-64 but should be on par with and likely faster than other libraries. Besdies that, I would guess IPP/KFR are dispatching to AVX2 while zldsp::fft is forced on SSE2, otherwise you shall not see a big gap.

It may get even faster if you use SoA for complex numbers (especially around medium orders). And SoA would also benifit spectrum processing a lot.

Yes, IPP is dispatching to AVX/AVX2. Here are the results for real FFT I have on my Ryzen 5800 (AoS):

Intel IPP ippSP AVX2 (l9) 2026.0.0 (r0xa7ad6ebc), enabled: avx2, cpu: avx2

|  size |  c/s |  fwd |  inv | GFLOP |  inv | name       | type   | snr dB |
| ----: | ---: | ---: | ---: | ----: | ---: | ---------- | ------ | -----: |
|  8192 |  5.4 |  2.5 |  3.0 |  43.9 | 36.9 | Intel IPP  | float  | -135.0 |
|  8192 |  9.5 |  4.7 |  4.8 |  24.6 | 24.1 | PFFFT      | float  | -135.1 |
|  8192 |  8.2 |  4.0 |  4.2 |  27.6 | 27.7 | zldsp::fft | float  | -134.2 |
|  8192 |  9.2 |  4.3 |  4.9 |  26.5 | 22.8 | Intel IPP  | double | -309.5 |
|  8192 | 16.5 |  8.0 |  8.5 |  12.2 | 13.9 | PFFFT      | double | -309.6 |
|  8192 | 15.4 |  7.5 |  7.9 |  15.2 | 14.1 | zldsp::fft | double | -307.4 |
| 65536 |  6.1 |  2.8 |  3.3 |  48.1 | 44.8 | Intel IPP  | float  | -133.9 |
| 65536 | 11.1 |  5.5 |  5.6 |  25.6 | 24.4 | PFFFT      | float  | -134.0 |
| 65536 | 12.1 |  6.0 |  6.2 |  23.7 | 24.1 | zldsp::fft | float  | -133.2 |
| 65536 | 11.6 |  5.5 |  6.0 |  25.8 | 23.4 | Intel IPP  | double | -307.7 |
| 65536 | 20.8 | 10.2 | 10.6 |  14.1 | 14.0 | PFFFT      | double | -308.6 |
| 65536 | 20.0 |  9.7 | 10.3 |  14.7 | 14.5 | zldsp::fft | double | -306.2 |

SoA is perhaps faster for FFT, but then you generally have to convert to AoS if you do any subsequent processing and want to benefit from accelerated functions which mostly expect AoS.

Techinically it is possible. And technically it is possible to dispatch the whole processBlock() function. It is just so painful to arrange all things in that way … A much easier solution is to compile two versions and detect the SIMD via app/plugin installer, which should also avoid any runtime dispatch overhead.

The runtime dispatch overhead should be negligible for any real-world FFT size, and infinitesmal with regards to the AVX2 speedup. Thats’ what most FFT libs do anyway.

Compiling two versions is technically feasible but a major PITA for developers and users: you have to do double the builds, tests, installers, etc. Plus most users have no idea what AVX even means.

I understand Highway provides a simple macro mechanism to do that, perhaps I could submit a PR.

Given that input/output might be pointed to the same address, I am not sure whether it is safe to do so, if the compiler does some optimization there?

Unless the data is marked `__restrict`, a sane compiler should assume the pointers could alias.
Usually libraries that can work in place (is this the case of zldsp::fft ?), they usually provide two interfaces, one explictely marked ‘in place’ with a single T buffer, and another ‘out of place’ with `const T*, T*`

[zsliu98]

@lorcan

Yes, IPP is dispatching to AVX/AVX2. Here are the results for real FFT I have on my Ryzen 5800 (AoS):

Exactly. If both are AVX2, I would except around 20% max gap.

@lorcan

SoA is perhaps faster for FFT, but then you generally have to convert to AoS if you do any subsequent processing and want to benefit from accelerated functions which mostly expect AoS.

If the functions come from an external library that you cannot control, that makes sense. Though, from RAW SIMD perspective (or if you write the functions in Google Highway), SoA is very easy to load/store while AoS requies expensive shuffling (especially on x86-64). And I think that’s why vDSP chooses split-complex (SoA) in the first place.

@lorcan

Compiling two versions is technically feasible but a major PITA for developers and users: you have to do double the builds, tests, installers, etc. Plus most users have no idea what AVX even means.

My idea is that there is only ONE installer but could detect the CPU SIMD support at the installation time, then users don’t need to bother with it. I haven’t tested it though as I am pretty sure I cannot use pure MSI for this. If you are a big company that uses a installation portal app, that should be super easy.

@lorcan

I understand Highway provides a simple macro mechanism to do that, perhaps I could submit a PR.

Feel free to try, but if you ever use more Highway code on both DSP and UI side, the refactoring won’t be very easy as far as I know. That’s why I stop at that point.

@lorcan

Unless the data is marked `__restrict`, a sane compiler should assume the pointers could alias.
Usually libraries that can work in place (is this the case of zldsp::fft ?), they usually provide two interfaces, one explictely marked ‘in place’ with a single T buffer, and another ‘out of place’ with `const T*, T*`

zldsp::fft supports both cases as the data will be copied to the internal (mostly AoSoA) structure.

[yairadix]

It is just so painful to arrange all things in that way

It shouldn’t be too painful. They made some nice crazy C preprocessor shenanigans to include the .cpp file several times with different compiler flags and different HWY_NAMESPACE values creating symbols for different backends, and a HWY_EXPORT macro to create the dynamic dispatch function for it.
The alternative of several binaries seems like a much bigger pain to the users imho

[zsliu98]

@yairadix

It shouldn’t be too painful. They made some nice crazy C preprocessor shenanigans to include the .cpp file several times with different compiler flags and different HWY_NAMESPACE values creating symbols for different backends, and a HWY_EXPORT macro to create the dynamic dispatch function for it.

Perhaps I misread the Highway document last time I will revisit it if possible, and see which one is easier to implement: runtime dispatch or detection via installer.

[zsliu98]

Done (with quite some help from LLM). I still don’t like the dispatch interface, so I keep both static & dynamic dispatch.

To keep the library header-only and make it possible to combine with other Highway dispatch code, I choose the style of bringing-your-own-dispatch-caller. I have also provided a dynamic dispatch caller as an example See the README/website for details.

[zsliu98]

Optimize high-order performance (a little bit) and add RFFT convolution APIs So for pure RFFT convolution you can avoid one working space & have slightly better performance.

Benchmark code has also been updated. Now benchmark code has much longer execution time & measures the native transformation of each library (i.e., the RFFT results are allowed to be in different formats).

Benchmark - High-performance FFT Library zldsp::fft

[lorcan]

@zsliu98 Thanks for implementing the dynamic dispatch feature !
Unfortunately on my setup this reduces performance dramatically, especially in double precision.

I profiled two RFFT (forward + backward) runs and pasted the results below.
As can be seen most time for AVX2 is spent in reconstruct_macro_twiddle_tile, so perhaps this function doesn’t get optimized properly ?

Static | SSE2

|  size |   c/s |   fwd |   inv | GFLOP |  inv |         name |   type | snr dB |    time |
| ----: | ----: | ----: | ----: | ----: | ---: | ------------ | ------ | -----: | ------: |
| 16384 |   5.6 |   2.6 |   3.0 |  47.1 | 40.7 |    Intel IPP |  float | -134.6 |  26.3us |
| 16384 |   9.9 |   4.9 |   5.0 |  26.3 | 25.0 |        PFFFT |  float | -134.8 |  44.8us |
| 16384 |   8.7 |   4.3 |   4.5 |  29.8 | 28.2 |   zldsp::fft |  float | -133.8 |  39.6us |
| 16384 |  10.0 |   4.7 |   5.3 |  27.0 | 23.9 |    Intel IPP | double | -308.1 |  45.2us |
| 16384 |  17.1 |   8.4 |   8.6 |  15.4 | 14.7 |        PFFFT | double | -309.4 |  76.4us |
| 16384 |  17.9 |   8.7 |   9.2 |  13.3 | 13.2 |   zldsp::fft | double | -307.0 |  86.5us |

Dynamic | AVX2

|  size |   c/s |   fwd |   inv | GFLOP |  inv |         name |   type | snr dB |    time |
| ----: | ----: | ----: | ----: | ----: | ---: | ------------ | ------ | -----: | ------: |
| 16384 |   6.5 |   2.7 |   3.8 |  46.5 | 33.0 |    Intel IPP |  float | -134.6 |  29.7us |
| 16384 |   9.9 |   5.0 |   4.9 |  25.8 | 26.0 |        PFFFT |  float | -134.8 |  44.3us |
| 16384 |  14.1 |   5.5 |   8.6 |  23.0 | 15.2 |   zldsp::fft |  float | -134.3 |  62.6us |
| 16384 |  10.2 |   5.3 |   5.0 |  24.2 | 21.8 |    Intel IPP | double | -308.1 |  50.0us |
| 16384 |  17.1 |   8.4 |   8.6 |  15.3 | 15.0 |        PFFFT | double | -309.4 |  75.7us |
| 16384 |  67.1 |  32.9 |  34.2 |   4.0 |  3.8 |   zldsp::fft | double | -307.2 | 292.5us |

[zsliu98]

If you are still on Ryzen 5800, for 16,384-point RFFT, the reconstruct_macro_twiddle_tile should not present in the float version. And if the AVX2 is significantly slower than SSE2, I would guess the dynamic dispatch is messed up.

May I know whether you can re-benchmark zldsp::fft under AVX2 static dispatch?

• MSVC: /arch:AVX2
• Clang: -march=x86-64-v3 -maes -mpclmul

[lorcan]

Yes reconstruct_macro_twiddle_tileis only indeed used for double.
I don’t have the AVX2 static dispatch exact numbers anymore, but they were on par with Intel IPP.
I can re-run the test if you need the exact numbers.
So as you say it does seem the dynamic dispatch is not working as intended.

[zsliu98]

If you could re-run the AVX2 static dispatch test, it would be really helpful as I can safely conclude the results.

Even given the float32 version results alone, the dynamic dispatch is already very very suspicious. And the No.1 suspect seems to be MSVC as it is not covered by current throughput benchmark code. I will investigate it.

_____

For MSVC, we recommend compiling with /Gv to allow non-inlined functions to pass vector arguments in registers. If intending to use the AVX2 target together with half-width vectors (e.g. for PromoteTo), it is also important to compile with /arch:AVX2. This seems to be the only way to reliably generate VEX-encoded SSE instructions on MSVC. Sometimes MSVC generates VEX-encoded SSE instructions, if they are mixed with AVX, but not always, see DevCom-10618264. Otherwise, mixing VEX-encoded AVX2 instructions and non-VEX SSE may cause severe performance degradation. Unfortunately, with /arch:AVX2 option, the resulting binary will then require AVX2. Note that no such flag is needed for clang and GCC because they support target-specific attributes, which we use to ensure proper VEX code generation for AVX2 targets.

Well, it does seem to match Highway documentation regarding “severe performance degradation” … I will try to find a workaround.

[lorcan]

Attached are the full benchmark results, including the AVX2 static dispatch (had to zip those because this forum doesn’t accept markdown uploads).

Actually in the static case zldsp’s AVX2 performance is a bit worse than with SSE1/SSE2, which looks a little weird too. I wonder if the CPU might be throttling down the clock as can happen with AVX2 on certain CPU’s apparently

It’s unfortunate that MSVC doesn’t support #pragma push(__vectorcall) because that means you have to annotate every function with __vectorcallmanually.

bench_fft_x64.zip (5.5 KB)

[lorcan]

Out of curiosity I tried again using CMake FetchContent to get Google Highway instead of conan,
now I get better results with static dispatch (dynamic stays the same).

bench_fft_x64_static_avx2_cmake.zip (1.8 KB)

[zsliu98]

@lorcan

Out of curiosity I tried again using CMake FetchContent to get Google Highway instead of conan,
now I get better results with static dispatch (dynamic stays the same).

Conan seems to use a pre-built Highway binary, which might not be very suitable for this case? Just an assumption.

Anyway, I do find that under MSVC the dynamic dispatch AVX2 and static dispatch AVX2 path differ a lot while the dynamic path seems to contain quite some SSE2 instructions (which is a BIG red flag). So I would guess the mixture of SSE2 and AVX2 instructions leads to serious performance problem. Clang-cl fortunately does not have such problem.

It seems that for MSVC I need to build the FFT backend twice separately (with different flags) and uses Highway dynamic dispatch only to select the backend. That does seem to be a very very ugly way to solve the issue.