Audacity 3.2.0
SpectrumTransformer.cpp
Go to the documentation of this file.
1/**********************************************************************
2
3Audacity: A Digital Audio Editor
4
5SpectrumTransformer.cpp
6
7Edward Hui
8
9**********************************************************************/
10
11#include "SpectrumTransformer.h"
12
13#include <algorithm>
14#include "FFT.h"
15#include "WaveTrack.h"
16
18 eWindowFunctions inWindowType,
19 eWindowFunctions outWindowType,
20 size_t windowSize, unsigned stepsPerWindow,
21 bool leadingPadding, bool trailingPadding )
22: mWindowSize{ windowSize }
23, mSpectrumSize{ 1 + mWindowSize / 2 }
24, mStepsPerWindow{ stepsPerWindow }
25, mStepSize{ mWindowSize / mStepsPerWindow }
26, mLeadingPadding{ leadingPadding }
27, mTrailingPadding{ trailingPadding }
28, hFFT{ GetFFT(mWindowSize) }
29, mFFTBuffer( mWindowSize )
30, mInWaveBuffer( mWindowSize )
31, mOutOverlapBuffer( mWindowSize )
32, mNeedsOutput{ needsOutput }
33{
34 // Check preconditions
35
36 // Powers of 2 only!
37 wxASSERT(mWindowSize > 0 &&
38 0 == (mWindowSize & (mWindowSize - 1)));
39
40 wxASSERT(mWindowSize % mStepsPerWindow == 0);
41
42 wxASSERT(!(inWindowType == eWinFuncRectangular && outWindowType == eWinFuncRectangular));
43
44 // To do: check that inWindowType, outWindowType, and mStepsPerWindow
45 // are compatible for correct overlap-add reconstruction.
46
47 // Create windows as needed
48 if (inWindowType != eWinFuncRectangular) {
49 mInWindow.resize(mWindowSize);
50 std::fill(mInWindow.begin(), mInWindow.end(), 1.0f);
51 NewWindowFunc(inWindowType, mWindowSize, false, mInWindow.data());
52 }
53 if (outWindowType != eWinFuncRectangular) {
54 mOutWindow.resize(mWindowSize);
55 std::fill(mOutWindow.begin(), mOutWindow.end(), 1.0f);
56 NewWindowFunc(outWindowType, mWindowSize, false, mOutWindow.data());
57 }
58
59 // Must scale one or the other window so overlap-add
60 // comes out right
61 double denom = 0;
62 for (size_t ii = 0; ii < mWindowSize; ii += mStepSize) {
63 denom +=
64 (mInWindow.empty() ? 1.0 : mInWindow[ii])
65 *
66 (mOutWindow.empty() ? 1.0 : mOutWindow[ii]);
67 }
68 // It is ASSUMED that you have chosen window types and
69 // steps per window, so that this sum denom would be the
70 // same, starting the march anywhere from 0 to mStepSize - 1.
71 // Else, your overlap-add won't be right, and the transformer
72 // might not be an identity even when you do nothing to the
73 // spectra.
74
75 float *pWindow = 0;
76 if (!mInWindow.empty())
77 pWindow = mInWindow.data();
78 else if (!mOutWindow.empty())
79 pWindow = mOutWindow.data();
80 else
81 // Can only happen if both window types were rectangular
82 wxASSERT(false);
83 for (size_t ii = 0; ii < mWindowSize; ++ii)
84 *pWindow++ /= denom;
85}
86
87auto SpectrumTransformer::NewWindow(size_t windowSize)
88 -> std::unique_ptr<Window>
89{
90 return std::make_unique<Window>(windowSize);
91}
92
94{
95 return true;
96}
97
99{
100 return true;
101}
102
103void
104TrackSpectrumTransformer::DoOutput(const float *outBuffer, size_t mStepSize)
105{
106 mOutputTrack->Append((constSamplePtr)outBuffer, floatSample, mStepSize);
107}
108
109bool SpectrumTransformer::Start(size_t queueLength)
110{
111 // Prepare clean queue
112 ResizeQueue(queueLength);
113 for (auto &pWindow : mQueue)
114 pWindow->Zero();
115
116 // invoke derived method
117 if (!DoStart())
118 return false;
119
120 // Clean input and output buffers
121 {
122 float *pFill;
123 pFill = mInWaveBuffer.data();
124 std::fill(pFill, pFill + mWindowSize, 0.0f);
125 pFill = mOutOverlapBuffer.data();
126 std::fill(pFill, pFill + mWindowSize, 0.0f);
127 }
128
129 if (!mLeadingPadding)
130 {
131 // We do not want leading zero padded windows
132 mInWavePos = 0;
133 mOutStepCount = -static_cast<int>(queueLength - 1);
134 }
135 else
136 {
137 // So that the queue gets primed with some windows,
138 // zero-padded in front, the first having mStepSize
139 // samples of wave data:
141 // This starts negative, to count up until the queue fills:
142 mOutStepCount = -static_cast<int>(queueLength - 1)
143 // ... and then must pass over the padded windows,
144 // before the first full window:
145 - static_cast<int>(mStepsPerWindow - 1);
146 }
147
148 mInSampleCount = 0;
149
150 return true;
151}
152
154 const float *buffer, size_t len )
155{
156 if (buffer)
157 mInSampleCount += len;
158 bool success = true;
159 while (success && len &&
160 mOutStepCount * static_cast<int>(mStepSize) < mInSampleCount) {
161 auto avail = std::min(len, mWindowSize - mInWavePos);
162 if (buffer)
163 memmove(&mInWaveBuffer[mInWavePos], buffer, avail * sizeof(float));
164 else
165 memset(&mInWaveBuffer[mInWavePos], 0, avail * sizeof(float));
166 if (buffer)
167 buffer += avail;
168 len -= avail;
169 mInWavePos += avail;
170
171 if (mInWavePos == mWindowSize) {
173
174 // invoke derived method
175 if ( (success = processor(*this)), success )
176 OutputStep();
177
180
181 // Shift input.
182 memmove(mInWaveBuffer.data(), &mInWaveBuffer[mStepSize],
183 (mWindowSize - mStepSize) * sizeof(float));
185 }
186 }
187
188 return success;
189}
190
191void SpectrumTransformer::ResizeQueue(size_t queueLength)
192{
193 int oldLen = mQueue.size();
194 mQueue.resize(queueLength);
195 for (size_t ii = oldLen; ii < queueLength; ++ii)
196 // invoke derived method to get a queue element
197 // with appropriate extra fields
199}
200
202{
203 // Transform samples to frequency domain, windowed as needed
204 {
205 auto pFFTBuffer = mFFTBuffer.data(), pInWaveBuffer = mInWaveBuffer.data();
206 if (mInWindow.size() > 0) {
207 auto pInWindow = mInWindow.data();
208 for (size_t ii = 0; ii < mWindowSize; ++ii)
209 *pFFTBuffer++ = *pInWaveBuffer++ * *pInWindow++;
210 }
211 else
212 memmove(pFFTBuffer, pInWaveBuffer, mWindowSize * sizeof(float));
213 }
214 RealFFTf(mFFTBuffer.data(), hFFT.get());
215
216 auto &record = Nth(0);
217
218 // Store real and imaginary parts for later inverse FFT
219 {
220 float *pReal = &record.mRealFFTs[1];
221 float *pImag = &record.mImagFFTs[1];
222 int *pBitReversed = &hFFT->BitReversed[1];
223 const auto last = mSpectrumSize - 1;
224 for (size_t ii = 1; ii < last; ++ii) {
225 const int kk = *pBitReversed++;
226 *pReal++ = mFFTBuffer[kk];
227 *pImag++ = mFFTBuffer[kk + 1];
228 }
229 // DC and Fs/2 bins need to be handled specially
230 const float dc = mFFTBuffer[0];
231 record.mRealFFTs[0] = dc;
232
233 const float nyquist = mFFTBuffer[1];
234 record.mImagFFTs[0] = nyquist; // For Fs/2, not really imaginary
235 }
236}
237
239{
240 std::rotate(mQueue.begin(), mQueue.end() - 1, mQueue.end());
241}
242
244{
245 bool bLoopSuccess = true;
246 if (mTrailingPadding) {
247 // Keep flushing empty input buffers through the history
248 // windows until we've output exactly as many samples as
249 // were input.
250 // Well, not exactly, but not more than one step-size of extra samples
251 // at the end.
252
253 while (bLoopSuccess &&
254 mOutStepCount * static_cast<int>(mStepSize) < mInSampleCount)
255 bLoopSuccess = ProcessSamples(processor, nullptr, mStepSize);
256 }
257
258 if (bLoopSuccess)
259 // invoke derived method
260 bLoopSuccess = DoFinish();
261
262 return bLoopSuccess;
263}
264
266{
267 auto allocSize = mQueue.size();
268 auto size = mOutStepCount + allocSize;
269 if (mLeadingPadding)
270 size += mStepsPerWindow - 1;
271
272 if (size < allocSize)
273 return size.as_size_t();
274 else
275 return allocSize;
276}
277
278// Formerly part of EffectNoiseReduction::Worker::ReduceNoise()
280{
281 if (!mNeedsOutput)
282 return;
283 if (QueueIsFull()) {
284 const auto last = mSpectrumSize - 1;
285 Window &record = **mQueue.rbegin();
286
287 const float *pReal = &record.mRealFFTs[1];
288 const float *pImag = &record.mImagFFTs[1];
289 float *pBuffer = &mFFTBuffer[2];
290 auto nn = mSpectrumSize - 2;
291 for (; nn--;) {
292 *pBuffer++ = *pReal++;
293 *pBuffer++ = *pImag++;
294 }
295 mFFTBuffer[0] = record.mRealFFTs[0];
296 // The Fs/2 component is stored as the imaginary part of the DC component
297 mFFTBuffer[1] = record.mImagFFTs[0];
298
299 // Invert the FFT into the output buffer
300 InverseRealFFTf(mFFTBuffer.data(), hFFT.get());
301
302 // Overlap-add
303 if (mOutWindow.size() > 0) {
304 auto pOut = mOutOverlapBuffer.data();
305 auto pWindow = mOutWindow.data();
306 auto pBitReversed = &hFFT->BitReversed[0];
307 for (size_t jj = 0; jj < last; ++jj) {
308 auto kk = *pBitReversed++;
309 *pOut++ += mFFTBuffer[kk] * (*pWindow++);
310 *pOut++ += mFFTBuffer[kk + 1] * (*pWindow++);
311 }
312 }
313 else {
314 auto pOut = mOutOverlapBuffer.data();
315 auto pBitReversed = &hFFT->BitReversed[0];
316 for (size_t jj = 0; jj < last; ++jj) {
317 auto kk = *pBitReversed++;
318 *pOut++ += mFFTBuffer[kk];
319 *pOut++ += mFFTBuffer[kk + 1];
320 }
321 }
322 auto buffer = mOutOverlapBuffer.data();
323 if (mOutStepCount >= 0) {
324 // Output the first portion of the overlap buffer, they're done
325 DoOutput(buffer, mStepSize);
326 }
327 // Shift the remainder over.
328 memmove(buffer, buffer + mStepSize, sizeof(float)*(mWindowSize - mStepSize));
329 std::fill(buffer + mWindowSize - mStepSize, buffer + mWindowSize, 0.0f);
330 }
331}
332
334{
335 if (mLeadingPadding)
336 return (mOutStepCount >= -static_cast<int>(mStepsPerWindow - 1));
337 else
338 return (mOutStepCount >= 0);
339}
340
342 WaveTrack *track, size_t queueLength, sampleCount start, sampleCount len)
343{
344 if (!track)
345 return false;
346
347 mpTrack = track;
348
349 if (!Start(queueLength))
350 return false;
351
352 mStart = start;
353 mLen = len;
354 auto bufferSize = track->GetMaxBlockSize();
355 FloatVector buffer(bufferSize);
356
357 bool bLoopSuccess = true;
358 auto samplePos = start;
359 while (bLoopSuccess && samplePos < start + len) {
360 //Get a blockSize of samples (smaller than the size of the buffer)
361 const auto blockSize = limitSampleBufferSize(
362 std::min(bufferSize, track->GetBestBlockSize(samplePos)),
363 start + len - samplePos);
364
365 //Get the samples from the track and put them in the buffer
366 track->GetFloats(buffer.data(), samplePos, blockSize);
367 samplePos += blockSize;
368
369 bLoopSuccess = ProcessSamples(processor, buffer.data(), blockSize);
370 }
371
372 if (!Finish(processor))
373 return false;
374
375 return bLoopSuccess;
376}
377
379{
380 if (mOutputTrack) {
381 // Flush the output WaveTrack (since it's buffered)
382 mOutputTrack->Flush();
383
384 // Take the output track and insert it in place of the original
385 // sample data
386 auto t0 = mOutputTrack->LongSamplesToTime(mStart);
387 auto tLen = mOutputTrack->LongSamplesToTime(mLen);
388 // Filtering effects always end up with more data than they started with.
389 // Delete this 'tail'.
390 mOutputTrack->HandleClear(tLen, mOutputTrack->GetEndTime(), false, false);
391 mpTrack->ClearAndPaste(t0, t0 + tLen, &*mOutputTrack, true, false);
392 }
393
394 mOutputTrack.reset();
395 return true;
396}
397
399
401
403
405{
406 mOutputTrack = NeedsOutput() && mpTrack ? mpTrack->EmptyCopy() : nullptr;
408}
int min(int a, int b)
void NewWindowFunc(int whichFunction, size_t NumSamplesIn, bool extraSample, float *in)
Definition: FFT.cpp:368
eWindowFunctions
Definition: FFT.h:110
@ eWinFuncRectangular
Definition: FFT.h:111
void RealFFTf(fft_type *buffer, const FFTParam *h)
Definition: RealFFTf.cpp:161
void InverseRealFFTf(fft_type *buffer, const FFTParam *h)
Definition: RealFFTf.cpp:263
HFFT GetFFT(size_t fftlen)
Definition: RealFFTf.cpp:104
size_t limitSampleBufferSize(size_t bufferSize, sampleCount limit)
Definition: SampleCount.cpp:22
const char * constSamplePtr
Definition: SampleFormat.h:56
bool GetFloats(float *buffer, sampleCount start, size_t len, fillFormat fill=fillZero, bool mayThrow=true, sampleCount *pNumWithinClips=nullptr) const
Retrieve samples from a track in floating-point format, regardless of the storage format.
Definition: SampleTrack.h:71
sampleCount mOutStepCount
sometimes negative
FloatVector mFFTBuffer
These have size mWindowSize:
std::function< bool(SpectrumTransformer &) > WindowProcessor
Type of function that transforms windows in the queue.
bool Finish(const WindowProcessor &processor)
Call once after a sequence of calls to ProcessSamples; flushes the queue and Invokes DoFinish.
void ResizeQueue(size_t queueLength)
virtual bool DoStart()
Called before any calls to ProcessWindow.
virtual void DoOutput(const float *outBuffer, size_t mStepSize)=0
Called within ProcessSamples if output was requested.
std::vector< float > FloatVector
const unsigned mStepsPerWindow
bool ProcessSamples(const WindowProcessor &processor, const float *buffer, size_t len)
Call multiple times.
virtual std::unique_ptr< Window > NewWindow(size_t windowSize)
Allocates a window to place in the queue.
Window & Nth(int n)
Access the queue, so you can inspect and modify any window in it.
std::vector< std::unique_ptr< Window > > mQueue
virtual ~SpectrumTransformer()
virtual bool DoFinish()
Called after the last call to ProcessWindow().
SpectrumTransformer(bool needsOutput, eWindowFunctions inWindowType, eWindowFunctions outWindowType, size_t windowSize, unsigned stepsPerWindow, bool leadingPadding, bool trailingPadding)
FloatVector mInWindow
These have size mWindowSize, or 0 for rectangular window:
size_t CurrentQueueSize() const
How many windows in the queue have been filled?
bool Start(size_t queueLength)
Call once before a sequence of calls to ProcessSamples; Invokes DoStart.
std::shared_ptr< WaveTrack > mOutputTrack
~TrackSpectrumTransformer() override
bool Process(const WindowProcessor &processor, WaveTrack *track, size_t queueLength, sampleCount start, sampleCount len)
Invokes Start(), ProcessSamples(), and Finish()
void DoOutput(const float *outBuffer, size_t mStepSize) override
Called within ProcessSamples if output was requested.
bool DoFinish() override
Called after the last call to ProcessWindow().
bool DoStart() override
Called before any calls to ProcessWindow.
A Track that contains audio waveform data.
Definition: WaveTrack.h:57
size_t GetMaxBlockSize() const override
This returns a nonnegative number of samples meant to size a memory buffer.
Definition: WaveTrack.cpp:1800
size_t GetBestBlockSize(sampleCount t) const override
This returns a nonnegative number of samples meant to size a memory buffer.
Definition: WaveTrack.cpp:1782
void ClearAndPaste(double t0, double t1, const Track *src, bool preserve=true, bool merge=true, const TimeWarper *effectWarper=NULL)
Definition: WaveTrack.cpp:899
Holder EmptyCopy(const SampleBlockFactoryPtr &pFactory={}, bool keepLink=true) const
Definition: WaveTrack.cpp:689
Positions or offsets within audio files need a wide type.
Definition: SampleCount.h:19
Derive this class to add information to the queue.
FloatVector mRealFFTs
index zero holds the dc coefficient, which has no imaginary part
FloatVector mImagFFTs
index zero holds the nyquist frequency coefficient, actually real