Use pffft in spectrogram #5960

Paul-Licameli · 2024-02-15T15:13:38Z

Resolves: #5880

Depends on

Fix misaligned pffft #5950

Serial performance improvement in spectrogram drawing, using pffft

I signed CLA
The title of the pull request describes an issue it addresses
If changes are extensive, then there is a sequence of easily reviewable commits
Each commit's message describes its purpose and effects
There are no behavior changes unnecessary for the stated purpose of the PR

Recommended:

Each commit compiles and runs on my machine without known undesirable changes of behavior

saintmatthieu

I think it'd be a relatively high-impact change for our users, but unfortunately it is a lot to consider, all the more that there are several algorithms to review and test.

saintmatthieu · 2024-02-22T09:49:08Z

src/effects/EqualizationFilter.cpp

-      im=buffer[hFFT->BitReversed[i]+1];
-      scratch[2*i  ] = re*mFilterFuncR[i] - im*mFilterFuncI[i];
-      scratch[2*i+1] = re*mFilterFuncI[i] + im*mFilterFuncR[i];
+   if(0) {


I suppose this if(0) is here because of the missing down-scaling upon IFFT?

I pushed that fix, looks like it works now.

saintmatthieu · 2024-02-22T10:46:38Z

src/tracks/playabletrack/wavetrack/ui/SpectrumCache.cpp

+      size_t i;
+      for (i = 0; i < len; ++i)
+         buffer[i] *= window[i];
+      for( ; i < len; ++i)


That doesn't look right. This second loop won't have any iteration.
On the other hand, looking in SpectrogramSettings.cpp, I find that window already includes the zero-padding, and that hFFT->Points is len/2. So I guess we can safely ignore padding here.

Thank you for noticing that! I don't want to let that stand.

saintmatthieu · 2024-02-22T11:09:13Z

src/tracks/playabletrack/wavetrack/ui/SpectrumCache.cpp

      else
-         out[i] = 10.0*log10f(power);
+         out[index / 2] = 10.0 * log10f(power);


Hey, maybe we could replace this log10f with an optimized log2. We saw in the tempo-estimation optimization PR that it sped up the STFT a great deal. (We probably should not blindly reuse the GetLog2 of MirDsp.cpp, though, because that approximation is only accurate to two decimal places, as you correctly found out, @Paul-Licameli ).

No.

Don’t use approximation for speed. Plot the best spectral data we can calculate.

saintmatthieu · 2024-02-22T11:21:53Z

src/tracks/playabletrack/wavetrack/ui/SpectrumCache.cpp

   const size_t scratchSize =
-      reassignment ? 3u * PffftAlignedCount{ fftLen } : fftLen;
+      (reassignment ? 4u : 2u) * PffftAlignedCount{ fftLen };


We're here in Populate, and within that function it's difficult to find out why there's a need for more than just fftLen. One has to dig CalculateOneSpectrum to find out about that dependency. And in fact, it doesn't seem that the standard spectrogram view requires scratch.

... Now I see that the buffer allocation is done outside a nested loop, which might explain why the scratch size must be determined outside of CalculateOneSpectrum.

I think the comment before this line is sufficient explanation.

... also more type checking that buffers are aligned. Also making sure not to inline any constants or functions in pffft.h, so that the libsoxr version isn't accidentally used by calling code. The full generality of complex time domain data, or the non-reordering calls, are not yet needed or exposed through this wrapper.

saintmatthieu

This is an important improvement and the type safety is good, but I'd like this ambiguity in the PffftFloatVector API to be resolved, at least by making that aligned overload private.

saintmatthieu · 2024-02-27T06:49:54Z

libraries/lib-fft/PffftTransformer.cpp

@@ -3,12 +3,12 @@

  Audacity: A Digital Audio Editor

-  PowerSpectrumGetter.cpp


At this stage you might want to keep the PowerSpectrumGetter in its own file and let PffftTransformer.h/cpop have their new Paul Licameli banners. The functionality you've authored covers 90% of this file. At the very least please add your name to the banner.

I think the power spectrum is a closely enough related utility that it can live in the same file.

saintmatthieu · 2024-02-27T06:58:12Z

src/tracks/playabletrack/wavetrack/ui/SpectrumCache.cpp

   const size_t scratchSize =
-      reassignment ? 3u * PffftAlignedCount{ fftLen } : fftLen;
+      (reassignment ? 4u : 2u) * PffftAlignedCount{ fftLen };


... Now I see that the buffer allocation is done outside a nested loop, which might explain why the scratch size must be determined outside of CalculateOneSpectrum.

saintmatthieu · 2024-02-27T07:31:12Z

libraries/lib-fft/PowerSpectrumGetter.h

   PffftAlignedCount &operator =(const PffftAlignedCount &) = default;

   //! @invariant `result: result % FloatAlignment == 0`
-   operator size_t() const { return value; }
+   operator size_t() const noexcept { return value; }

 private:
   size_t value{};


This should be const, it would seem.
And then, at least on my Windows machine with build settings that are still inconsistent with CI, I can also constexpr operator size_t() const noexcept { return value; } ; or are there platforms where it fails?

Making these const and constexpr allowed me to write this:

static_assert(PffftAlignedCount(0) == 0); static_assert(PffftAlignedCount(1) == 16); static_assert(PffftAlignedCount(16) == 16); static_assert(PffftAlignedCount(17) == 32);

What I understand from this is that this PffftAlignedCount(n) answers the question "How many sizeof(float) bytes do I need for n floats?". This is useful if PffftFloatVector is to be used as a vector of vectors, or matrix, but not as a straightforward vector; something like

PffftFloatVector v(2); v[1] = 2.f; const auto c1 = v.aligned(PffftAlignedCount(1));

doesn't make sense: c1.get() is out of bound and *c1.get() isn't 2.
In fact, I found out that the PffftFloats aligned( overload with one argument is not used anywhere other than PffftFloatVector, so it can - and should - be made private.

I was explicit in comments above the class:

Usual hazards of pointer invalidation or out-of-bounds subscripting still apply. Size and capacity are NOT constrained to be aligned.

I disagree that the one-argument aligned should be made private.

There may be future need.

saintmatthieu · 2024-02-27T07:57:09Z

libraries/lib-fft/PowerSpectrumGetter.cpp

+   pffft_transform_ordered(get(), input.get(), out, work.get(),
+      PFFFT_BACKWARD);
+   if (renormalize)
+      std::transform(out, out + N, out, [N = N](float f){ return f / N; });


Requires #include <algorithm>

saintmatthieu · 2024-02-27T08:32:51Z

src/effects/EqualizationFilter.cpp

-   // Inverse FFT and normalization
-   InverseRealFFTf(scratch, hFFT.get());
-   ReorderToTime(hFFT.get(), scratch, buffer);
+   transformer.InverseTransformOrdered(buf, buf, scr, true);


anonymous boolean

One easily jumps to the function definition to know what the boolean is, if the IDE browser is good enough.

saintmatthieu

Indeed, @Paul-Licameli, using a spectrogram length of 8 in the settings causes an assertion in PFFFT.
Reading more carefully pffft.h:

   - supports only transforms for inputs of length N of the form
   N=(2^a)*(3^b)*(5^c), a >= 5, b >=0, c >= 0 (32, 48, 64, 96, 128,
   144, 160, etc are all acceptable lengths). Performance is best for
   128<=N<=8192.

This would imply that the least size is 32. (Maybe by luck, but I just tried a size of 16 and that worked.)

@LWinterberg do you have an idea of the amount of dissatisfaction raising the least spectrum size to 32 would cause? I never came across a need for such small window sizes in audio applications. Maybe analysis of non-audio waves sampled at extremely low rates ?

PS: I wasn't aware that non-power-of-two sizes were also possible - interesting!

LWinterberg · 2024-02-27T09:53:37Z

@LWinterberg do you have an idea of the amount of dissatisfaction raising the least spectrum size to 32 would cause? I never came across a need for such small window sizes in audio applications. Maybe analysis of non-audio waves sampled at extremely low rates ?

After I didn't find a a spectrogram in promising places like Multi-Beam Spatio-Spectral Beamforming Receiver for Wideband Phased Arrays, Pipeline FFT architectures optimized for FPGAs and Goldstein phase filtering on InSAR data (it's an interferogram), I checked Sonic Visualiser and found:

If even they don't need anything below 32, it's probably fair to say we don't either.

Paul-Licameli · 2024-02-27T11:23:28Z

Indeed, @Paul-Licameli, using a spectrogram length of 8 in the settings causes an assertion in PFFFT. Reading more carefully pffft.h:
   - supports only transforms for inputs of length N of the form
   N=(2^a)*(3^b)*(5^c), a >= 5, b >=0, c >= 0 (32, 48, 64, 96, 128,
   144, 160, etc are all acceptable lengths). Performance is best for
   128<=N<=8192.
This would imply that the least size is 32. (Maybe by luck, but I just tried a size of 16 and that worked.)

@LWinterberg do you have an idea of the amount of dissatisfaction raising the least spectrum size to 32 would cause? I never came across a need for such small window sizes in audio applications. Maybe analysis of non-audio waves sampled at extremely low rates ?

PS: I wasn't aware that non-power-of-two sizes were also possible - interesting!

Interesting, they allow other prime factors than 2. See Lanczos Lemma.

The divide and conquer step of fft.

saintmatthieu

Unsatisfied to see argued suggestions rejected for no apparent good reason, but the remaining time for completion of this PR doesn't allow further argument. Please raise the least allowed frame size to 32 and I'll approve.

saintmatthieu · 2024-02-28T08:29:32Z

libraries/lib-fft/PffftTransformer.cpp

@@ -6,6 +6,7 @@
  PffftTransformer.cpp

  Matthieu Hodgkinson
+  Paul Licameli


Funny character

saintmatthieu · 2024-02-28T08:30:52Z

libraries/lib-fft/PffftTransformer.h

@@ -68,7 +68,7 @@ struct PffftAlignedCount
   PffftAlignedCount &operator =(const PffftAlignedCount &) = default;

   //! @invariant `result: result % FloatAlignment == 0`
-   operator size_t() const noexcept { return value; }
+   constexpr operator size_t() const noexcept { return value; }


value could/should still be const.

saintmatthieu · 2024-02-29T08:35:59Z

libraries/lib-fft/PffftTransformer.cpp

+      auto pOutput = output.get();
+      const auto factor = N / requestedSize;
+      for (size_t ii = 1, nn = requestedSize; ii < nn; ++ii)
+         pOutput[ii] = pOutput[ii * factor];


You overlooked that output here has layout [DC, Nyq, real1, imag1, ...].
But besides that I think it should work.

Paul-Licameli requested review from crsib and saintmatthieu February 15, 2024 15:13

Paul-Licameli self-assigned this Feb 15, 2024

Paul-Licameli force-pushed the Use-pffft-in-spectrogram branch 3 times, most recently from aaad77d to 97d3fe6 Compare February 20, 2024 12:24

saintmatthieu reviewed Feb 22, 2024

View reviewed changes

Paul-Licameli force-pushed the Use-pffft-in-spectrogram branch 2 times, most recently from 9789675 to 9119270 Compare February 24, 2024 11:18

Paul-Licameli added 9 commits February 26, 2024 07:10

constexpr an noexcept in PowerSpectrumGetter.h

9dee5c3

Rename files to PffftTransformer.*

0cb741c

Aligned allocation of input data when computing spectrograms

f3261ef

Aligned allocation of output data when computing spectrograms

4a2a841

Spectrogram drawing uses pffft

4db6edf

Filter wasn't deep-const! Externalize the scratch to fix that

d428a5f

Alignment of buffers used in Equalization

ff14885

Use pffft for the EQ effect processing

305c96d

Paul-Licameli force-pushed the Use-pffft-in-spectrogram branch from 683e60a to 305c96d Compare February 26, 2024 14:02

Use Pffft in SpectrumTransformer (thus in noise reduction)

2c55e9b

saintmatthieu reviewed Feb 27, 2024

View reviewed changes

saintmatthieu self-requested a review February 27, 2024 08:49

saintmatthieu requested changes Feb 27, 2024

View reviewed changes

Paul-Licameli added 5 commits February 27, 2024 16:21

Be sure zero padding is correct

d7607d2

Update banners

07c7694

Fix Windows, Linux builds

4c69c1b

More constexpr

5c4465e

comments

04efbd7

saintmatthieu requested changes Feb 28, 2024

View reviewed changes

Paul-Licameli added 2 commits February 28, 2024 06:08

More constexpr, fixes the Linux build

aabfdbe

Attempt to make the smallest spectrogram window sizes still work

592d80b

saintmatthieu reviewed Feb 29, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Use pffft in spectrogram #5960

Use pffft in spectrogram #5960

Paul-Licameli commented Feb 15, 2024 •

edited by saintmatthieu

saintmatthieu left a comment

saintmatthieu Feb 22, 2024

saintmatthieu Feb 22, 2024

saintmatthieu Feb 22, 2024

Paul-Licameli Feb 27, 2024

saintmatthieu Feb 22, 2024

Paul-Licameli Feb 27, 2024

saintmatthieu Feb 22, 2024

saintmatthieu Feb 27, 2024

Paul-Licameli Feb 27, 2024

saintmatthieu left a comment

saintmatthieu Feb 27, 2024

Paul-Licameli Feb 27, 2024

saintmatthieu Feb 27, 2024

saintmatthieu Feb 27, 2024

saintmatthieu Feb 27, 2024

Paul-Licameli Feb 27, 2024

Paul-Licameli Feb 27, 2024

saintmatthieu Feb 27, 2024

saintmatthieu Feb 27, 2024

Paul-Licameli Feb 27, 2024

saintmatthieu left a comment

LWinterberg commented Feb 27, 2024

Paul-Licameli commented Feb 27, 2024

saintmatthieu left a comment

saintmatthieu Feb 28, 2024

saintmatthieu Feb 28, 2024

saintmatthieu Feb 29, 2024

		@@ -3,12 +3,12 @@

		Audacity: A Digital Audio Editor

		PowerSpectrumGetter.cpp

Use pffft in spectrogram #5960

Are you sure you want to change the base?

Use pffft in spectrogram #5960

Conversation

Paul-Licameli commented Feb 15, 2024 • edited by saintmatthieu

saintmatthieu left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

saintmatthieu left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

saintmatthieu left a comment

Choose a reason for hiding this comment

LWinterberg commented Feb 27, 2024

Paul-Licameli commented Feb 27, 2024

saintmatthieu left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Paul-Licameli commented Feb 15, 2024 •

edited by saintmatthieu