New library for math...

... note the swap of target_link_libraries lines in src/CMakeLists.txt,
needed to build at least on macOS, becuase FFT.h must be looked up first in
lib-math, not in lib-src/twolame

Also making a dependency cycle of SampleFormat and Dither!  But we will tolerate
that within one small library.

libraries/CMakeLists.txt

@@ -11,6 +11,7 @@ set( LIBRARIES
    lib-basic-ui
    lib-exceptions
    lib-preferences
+   lib-math
 )
 
 if ( ${_OPT}has_networking )

libraries/lib-math/CMakeLists.txt

@@ -0,0 +1,31 @@
+#[[
+A library of mathematical utilities and manipulation of samples
+]]#
+
+set( SOURCES
+   Dither.cpp
+   Dither.h
+   FFT.cpp
+   FFT.h
+   InterpolateAudio.cpp
+   InterpolateAudio.h
+   Matrix.cpp
+   Matrix.h
+   RealFFTf.cpp
+   RealFFTf.h
+   SampleCount.cpp
+   SampleCount.h
+   SampleFormat.cpp
+   SampleFormat.h
+   Spectrum.cpp
+   Spectrum.h
+   float_cast.h
+)
+set( LIBRARIES
+   lib-preferences-interface
+   PRIVATE
+   wxBase
+)
+audacity_library( lib-math "${SOURCES}" "${LIBRARIES}"
+   "" ""
+)

libraries/lib-math/Dither.cpp

@@ -0,0 +1,441 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  Dither.cpp
+
+  Steve Harris
+  Markus Meyer
+
+*******************************************************************//*!
+
+\class Dither
+\brief
+
+  This class implements various functions for dithering and is derived
+  from the dither code in the Ardour project, written by Steve Harris.
+
+  Dithering is only done if it really is necessary. Otherwise (e.g.
+  when the source and destination format of the samples is the same),
+  the samples are only copied or converted. However, copied samples
+  are always checked for out-of-bounds values and possibly clipped
+  accordingly.
+
+  These dither algorithms are currently implemented:
+  - No dithering at all
+  - Rectangle dithering
+  - Triangle dithering
+  - Noise-shaped dithering
+
+Dither class. You must construct an instance because it keeps
+state. Call Dither::Apply() to apply the dither. You can call
+Reset() between subsequent dithers to reset the dither state
+and get deterministic behaviour.
+
+*//*******************************************************************/
+
+
+#include "Dither.h"
+
+#include "Internat.h"
+#include "Prefs.h"
+
+// Erik de Castro Lopo's header file that
+// makes sure that we have lrint and lrintf
+// (Note: this file should be included first)
+#include "float_cast.h"
+
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+//#include <sys/types.h>
+//#include <memory.h>
+//#include <assert.h>
+
+#include <wx/defs.h>
+
+//////////////////////////////////////////////////////////////////////////
+
+// Constants for the noise shaping buffer
+const int Dither::BUF_MASK = 7;
+const int Dither::BUF_SIZE = 8;
+
+// Lipshitz's minimally audible FIR
+const float Dither::SHAPED_BS[] = { 2.033f, -2.165f, 1.959f, -1.590f, 0.6149f };
+
+// This is supposed to produce white noise and no dc
+#define DITHER_NOISE (rand() / (float)RAND_MAX - 0.5f)
+
+// The following is a rather ugly, but fast implementation
+// of a dither loop. The macro "DITHER" is expanded to an implementation
+// of a dithering algorithm, which contains no branches in the inner loop
+// except the branches for clipping the sample, and therefore should
+// be quite fast.
+
+#if 0
+// To assist in understanding what the macros are doing, here's an example of what
+// the result would be for Shaped dither:
+//
+// DITHER(ShapedDither, dest, destFormat, destStride, source, sourceFormat, sourceStride, len);
+
+    do {
+        if (sourceFormat == int24Sample && destFormat == int16Sample)
+            do {
+                char *d, *s;
+                unsigned int i;
+                int x;
+                for (d = (char*)dest, s = (char*)source, i = 0; i < len; i++, d += (int16Sample >> 16) * destStride, s += (int24Sample >> 16) * sourceStride)
+                    do {
+                        x = lrintf((ShapedDither((((*(( int*)(s)) / float(1<<23))) * float(1<<15)))));
+                        if (x>(32767))
+                            *((short*)((((d)))))=(32767);
+                        else if (x<(-32768))
+                            *((short*)((((d)))))=(-32768);
+                        else *((short*)((((d)))))=(short)x;
+                    } while (0);
+            } while (0);
+        else if (sourceFormat == floatSample && destFormat == int16Sample)
+            do {
+                char *d, *s;
+                unsigned int i;
+                int x;
+                for (d = (char*)dest, s = (char*)source, i = 0; i < len; i++, d += (int16Sample >> 16) * destStride, s += (floatSample >> 16) * sourceStride)
+                    do {
+                        x = lrintf((ShapedDither((((*((float*)(s)) > 1.0 ? 1.0 : *((float*)(s)) < -1.0 ? -1.0 : *((float*)(s)))) * float(1<<15)))));
+                        if (x>(32767))
+                            *((short*)((((d)))))=(32767);
+                        else if (x<(-32768))
+                            *((short*)((((d)))))=(-32768);
+                        else *((short*)((((d)))))=(short)x;
+                    } while (0);
+            } while (0);
+        else if (sourceFormat == floatSample && destFormat == int24Sample)
+            do {
+                char *d, *s;
+                unsigned int i;
+                int x;
+                for (d = (char*)dest, s = (char*)source, i = 0; i < len; i++, d += (int24Sample >> 16) * destStride, s += (floatSample >> 16) * sourceStride)
+                    do {
+                        x = lrintf((ShapedDither((((*((float*)(s)) > 1.0 ? 1.0 : *((float*)(s)) < -1.0 ? -1.0 : *((float*)(s)))) * float(1<<23)))));
+                        if (x>(8388607))
+                            *((int*)((((d)))))=(8388607);
+                        else if (x<(-8388608))
+                            *((int*)((((d)))))=(-8388608);
+                        else *((int*)((((d)))))=(int)x;
+                    } while (0);
+            } while (0);
+        else {
+            if ( false )
+                ;
+            else wxOnAssert(L"c:\\users\\yam\\documents\\audacity\\mixer\\n\\audacity\\src\\dither.cpp", 348,  __FUNCTION__  , L"false", 0);
+        }
+    } while (0);
+#endif
+
+// Defines for sample conversion
+#define CONVERT_DIV16 float(1<<15)
+#define CONVERT_DIV24 float(1<<23)
+
+// Dereference sample pointer and convert to float sample
+#define FROM_INT16(ptr) (*((short*)(ptr)) / CONVERT_DIV16)
+#define FROM_INT24(ptr) (*((  int*)(ptr)) / CONVERT_DIV24)
+
+// For float, we internally allow values greater than 1.0, which
+// would blow up the dithering to int values.  FROM_FLOAT is
+// only used to dither to int, so clip here.
+#define FROM_FLOAT(ptr) (*((float*)(ptr)) >  1.0 ?  1.0 : \
+                         *((float*)(ptr)) < -1.0 ? -1.0 : \
+                         *((float*)(ptr)))
+
+// Promote sample to range of specified type, keep it float, though
+#define PROMOTE_TO_INT16(sample) ((sample) * CONVERT_DIV16)
+#define PROMOTE_TO_INT24(sample) ((sample) * CONVERT_DIV24)
+
+// Store float sample 'sample' into pointer 'ptr', clip it, if necessary
+// Note: This assumes, a variable 'x' of type int is valid which is
+//       used by this macro.
+#define IMPLEMENT_STORE(ptr, sample, ptr_type, min_bound, max_bound) \
+    do { \
+    x = lrintf(sample); \
+    if (x>(max_bound)) *((ptr_type*)(ptr))=(max_bound); \
+    else if (x<(min_bound)) *((ptr_type*)(ptr))=(min_bound); \
+    else *((ptr_type*)(ptr))=(ptr_type)x; } while (0)
+
+#define STORE_INT16(ptr, sample) IMPLEMENT_STORE((ptr), (sample), short, -32768, 32767)
+#define STORE_INT24(ptr, sample) IMPLEMENT_STORE((ptr), (sample), int, -8388608, 8388607)
+
+// Dither single float 'sample' and store it in pointer 'dst', using 'dither' as algorithm
+#define DITHER_TO_INT16(dither, dst, sample) STORE_INT16((dst), dither(PROMOTE_TO_INT16(sample)))
+#define DITHER_TO_INT24(dither, dst, sample) STORE_INT24((dst), dither(PROMOTE_TO_INT24(sample)))
+
+// Implement one single dither step
+#define DITHER_STEP(dither, store, load, dst, src) \
+    store(dither, (dst), load(src))
+
+// Implement a dithering loop
+// Note: The variable 'x' is needed for the STORE_... macros
+#define DITHER_LOOP(dither, store, load, dst, dstFormat, dstStride, src, srcFormat, srcStride, len) \
+    do { \
+       char *d, *s; \
+       unsigned int ii; \
+       int x; \
+       for (d = (char*)dst, s = (char*)src, ii = 0; \
+            ii < len; \
+            ii++, d += SAMPLE_SIZE(dstFormat) * dstStride, \
+                 s += SAMPLE_SIZE(srcFormat) * srcStride) \
+          DITHER_STEP(dither, store, load, d, s); \
+   } while (0)
+
+// Shortcuts to dithering loops
+#define DITHER_INT24_TO_INT16(dither, dst, dstStride, src, srcStride, len) \
+    DITHER_LOOP(dither, DITHER_TO_INT16, FROM_INT24, dst, int16Sample, dstStride, src, int24Sample, srcStride, len)
+#define DITHER_FLOAT_TO_INT16(dither, dst, dstStride, src, srcStride, len) \
+    DITHER_LOOP(dither, DITHER_TO_INT16, FROM_FLOAT, dst, int16Sample, dstStride, src, floatSample, srcStride, len)
+#define DITHER_FLOAT_TO_INT24(dither, dst, dstStride, src, srcStride, len) \
+    DITHER_LOOP(dither, DITHER_TO_INT24, FROM_FLOAT, dst, int24Sample, dstStride, src, floatSample, srcStride, len)
+
+// Implement a dither. There are only 3 cases where we must dither,
+// in all other cases, no dithering is necessary.
+#define DITHER(dither, dst, dstFormat, dstStride, src, srcFormat, srcStride, len) \
+    do { if (srcFormat == int24Sample && dstFormat == int16Sample) \
+        DITHER_INT24_TO_INT16(dither, dst, dstStride, src, srcStride, len); \
+    else if (srcFormat == floatSample && dstFormat == int16Sample) \
+        DITHER_FLOAT_TO_INT16(dither, dst, dstStride, src, srcStride, len); \
+    else if (srcFormat == floatSample && dstFormat == int24Sample) \
+        DITHER_FLOAT_TO_INT24(dither, dst, dstStride, src, srcStride, len); \
+    else { wxASSERT(false); } \
+    } while (0)
+
+
+Dither::Dither()
+{
+    // On startup, initialize dither by resetting values
+    Reset();
+}
+
+void Dither::Reset()
+{
+    mTriangleState = 0;
+    mPhase = 0;
+    memset(mBuffer, 0, sizeof(float) * BUF_SIZE);
+}
+
+// This only decides if we must dither at all, the dithers
+// are all implemented using macros.
+//
+// "source" and "dest" can contain either interleaved or non-interleaved
+// samples.  They do not have to be the same...one can be interleaved while
+// the otner non-interleaved.
+//
+// The "len" argument specifies the number of samples to process.
+//
+// If either stride value equals 1 then the corresponding buffer contains
+// non-interleaved samples.
+//
+// If either stride value is greater than 1 then the corresponding buffer
+// contains interleaved samples and they will be processed by skipping every
+// stride number of samples.
+
+void Dither::Apply(enum DitherType ditherType,
+                   constSamplePtr source, sampleFormat sourceFormat,
+                   samplePtr dest, sampleFormat destFormat,
+                   unsigned int len,
+                   unsigned int sourceStride /* = 1 */,
+                   unsigned int destStride /* = 1 */)
+{
+    unsigned int i;
+
+    // This code is not designed for 16-bit or 64-bit machine
+    wxASSERT(sizeof(int) == 4);
+    wxASSERT(sizeof(short) == 2);
+
+    // Check parameters
+    wxASSERT(source);
+    wxASSERT(dest);
+    wxASSERT(len >= 0);
+    wxASSERT(sourceStride > 0);
+    wxASSERT(destStride > 0);
+
+    if (len == 0)
+        return; // nothing to do
+
+    if (destFormat == sourceFormat)
+    {
+        // No need to dither, because source and destination
+        // format are the same. Just copy samples.
+        if (destStride == 1 && sourceStride == 1)
+            memcpy(dest, source, len * SAMPLE_SIZE(destFormat));
+        else
+        {
+            if (sourceFormat == floatSample)
+            {
+                auto d = (float*)dest;
+                auto s = (const float*)source;
+
+                for (i = 0; i < len; i++, d += destStride, s += sourceStride)
+                    *d = *s;
+            } else
+            if (sourceFormat == int24Sample)
+            {
+                auto d = (int*)dest;
+                auto s = (const int*)source;
+
+                for (i = 0; i < len; i++, d += destStride, s += sourceStride)
+                    *d = *s;
+            } else
+            if (sourceFormat == int16Sample)
+            {
+                auto d = (short*)dest;
+                auto s = (const short*)source;
+
+                for (i = 0; i < len; i++, d += destStride, s += sourceStride)
+                    *d = *s;
+            } else {
+                wxASSERT(false); // source format unknown
+            }
+        }
+    } else
+    if (destFormat == floatSample)
+    {
+        // No need to dither, just convert samples to float.
+        // No clipping should be necessary.
+        auto d = (float*)dest;
+
+        if (sourceFormat == int16Sample)
+        {
+            auto s = (const short*)source;
+            for (i = 0; i < len; i++, d += destStride, s += sourceStride)
+                *d = FROM_INT16(s);
+        } else
+        if (sourceFormat == int24Sample)
+        {
+            auto s = (const int*)source;
+            for (i = 0; i < len; i++, d += destStride, s += sourceStride)
+                *d = FROM_INT24(s);
+        } else {
+            wxASSERT(false); // source format unknown
+        }
+    } else
+    if (destFormat == int24Sample && sourceFormat == int16Sample)
+    {
+        // Special case when promoting 16 bit to 24 bit
+        auto d = (int*)dest;
+        auto s = (const short*)source;
+        for (i = 0; i < len; i++, d += destStride, s += sourceStride)
+            *d = ((int)*s) << 8;
+    } else
+    {
+        // We must do dithering
+        switch (ditherType)
+        {
+        case DitherType::none:
+            DITHER(NoDither, dest, destFormat, destStride, source, sourceFormat, sourceStride, len);
+            break;
+        case DitherType::rectangle:
+            DITHER(RectangleDither, dest, destFormat, destStride, source, sourceFormat, sourceStride, len);
+            break;
+        case DitherType::triangle:
+            Reset(); // reset dither filter for this NEW conversion
+            DITHER(TriangleDither, dest, destFormat, destStride, source, sourceFormat, sourceStride, len);
+            break;
+        case DitherType::shaped:
+            Reset(); // reset dither filter for this NEW conversion
+            DITHER(ShapedDither, dest, destFormat, destStride, source, sourceFormat, sourceStride, len);
+            break;
+        default:
+            wxASSERT(false); // unknown dither algorithm
+        }
+    }
+}
+
+// Dither implementations
+
+// No dither, just return sample
+inline float Dither::NoDither(float sample)
+{
+    return sample;
+}
+
+// Rectangle dithering, apply one-step noise
+inline float Dither::RectangleDither(float sample)
+{
+    return sample - DITHER_NOISE;
+}
+
+// Triangle dither - high pass filtered
+inline float Dither::TriangleDither(float sample)
+{
+    float r = DITHER_NOISE;
+    float result = sample + r - mTriangleState;
+    mTriangleState = r;
+
+    return result;
+}
+
+// Shaped dither
+inline float Dither::ShapedDither(float sample)
+{
+    // Generate triangular dither, +-1 LSB, flat psd
+    float r = DITHER_NOISE + DITHER_NOISE;
+    if(sample != sample)  // test for NaN
+       sample = 0; // and do the best we can with it
+
+    // Run FIR
+    float xe = sample + mBuffer[mPhase] * SHAPED_BS[0]
+        + mBuffer[(mPhase - 1) & BUF_MASK] * SHAPED_BS[1]
+        + mBuffer[(mPhase - 2) & BUF_MASK] * SHAPED_BS[2]
+        + mBuffer[(mPhase - 3) & BUF_MASK] * SHAPED_BS[3]
+        + mBuffer[(mPhase - 4) & BUF_MASK] * SHAPED_BS[4];
+
+    // Accumulate FIR and triangular noise
+    float result = xe + r;
+
+    // Roll buffer and store last error
+    mPhase = (mPhase + 1) & BUF_MASK;
+    mBuffer[mPhase] = xe - lrintf(result);
+
+    return result;
+}
+
+static const std::initializer_list<EnumValueSymbol> choicesDither{
+   { XO("None") },
+   { XO("Rectangle") },
+   { XO("Triangle") },
+   { XO("Shaped") },
+};
+static auto intChoicesDither = {
+   DitherType::none,
+   DitherType::rectangle,
+   DitherType::triangle,
+   DitherType::shaped,
+};
+
+EnumSetting< DitherType > Dither::FastSetting{
+   wxT("Quality/DitherAlgorithmChoice"),
+   choicesDither,
+   0, // none
+
+   // for migrating old preferences:
+   intChoicesDither,
+   wxT("Quality/DitherAlgorithm")
+};
+
+EnumSetting< DitherType > Dither::BestSetting{
+   wxT("Quality/HQDitherAlgorithmChoice"),
+   choicesDither,
+   3, // shaped
+
+   // for migrating old preferences:
+   intChoicesDither,
+   wxT("Quality/HQDitherAlgorithm")
+};
+
+DitherType Dither::FastDitherChoice()
+{
+   return (DitherType) FastSetting.ReadEnum();
+}
+
+DitherType Dither::BestDitherChoice()
+{
+   return (DitherType) BestSetting.ReadEnum();
+}

libraries/lib-math/Dither.h

@@ -0,0 +1,66 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  Steve Harris
+  Markus Meyer
+
+**********************************************************************/
+
+#ifndef __AUDACITY_DITHER_H__
+#define __AUDACITY_DITHER_H__
+
+#include "SampleFormat.h"
+
+template< typename Enum > class EnumSetting;
+
+
+/// These ditherers are currently available:
+enum DitherType : unsigned {
+   none = 0, rectangle = 1, triangle = 2, shaped = 3 };
+
+class MATH_API Dither
+{
+public:
+    static DitherType FastDitherChoice();
+    static DitherType BestDitherChoice();
+
+    static EnumSetting< DitherType > FastSetting;
+    static EnumSetting< DitherType > BestSetting;
+
+    /// Default constructor
+    Dither();
+
+    /// Reset state of the dither.
+    void Reset();
+
+    /// Apply the actual dithering. Expects the source sample in the
+    /// 'source' variable, the destination sample in the 'dest' variable,
+    /// and hints to the formats of the samples. Even if the sample formats
+    /// are the same, samples are clipped, if necessary.
+    void Apply(DitherType ditherType,
+               constSamplePtr source, sampleFormat sourceFormat,
+               samplePtr dest, sampleFormat destFormat,
+               unsigned int len,
+               unsigned int sourceStride = 1,
+               unsigned int destStride = 1);
+
+private:
+    // Dither methods
+    float NoDither(float sample);
+    float RectangleDither(float sample);
+    float TriangleDither(float sample);
+    float ShapedDither(float sample);
+
+    // Dither constants
+    static const int BUF_SIZE; /* = 8 */
+    static const int BUF_MASK; /* = 7 */
+    static const float SHAPED_BS[];
+
+    // Dither state
+    int mPhase;
+    float mTriangleState;
+    float mBuffer[8 /* = BUF_SIZE */];
+};
+
+#endif /* __AUDACITY_DITHER_H__ */

libraries/lib-math/FFT.cpp

@@ -0,0 +1,700 @@
+/**********************************************************************
+
+  FFT.cpp
+
+  Dominic Mazzoni
+
+  September 2000
+
+*******************************************************************//*!
+
+\file FFT.cpp
+\brief Fast Fourier Transform routines.
+
+  This file contains a few FFT routines, including a real-FFT
+  routine that is almost twice as fast as a normal complex FFT,
+  and a power spectrum routine when you know you don't care
+  about phase information.
+
+  Some of this code was based on a free implementation of an FFT
+  by Don Cross, available on the web at:
+
+    http://www.intersrv.com/~dcross/fft.html
+
+  The basic algorithm for his code was based on Numerican Recipes
+  in Fortran.  I optimized his code further by reducing array
+  accesses, caching the bit reversal table, and eliminating
+  float-to-double conversions, and I added the routines to
+  calculate a real FFT and a real power spectrum.
+
+*//*******************************************************************/
+/*
+  Salvo Ventura - November 2006
+  Added more window functions:
+    * 4: Blackman
+    * 5: Blackman-Harris
+    * 6: Welch
+    * 7: Gaussian(a=2.5)
+    * 8: Gaussian(a=3.5)
+    * 9: Gaussian(a=4.5)
+*/
+
+#include "FFT.h"
+
+#include "Internat.h"
+
+#include "SampleFormat.h"
+
+#include <wx/wxcrtvararg.h>
+#include <wx/intl.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "RealFFTf.h"
+
+static ArraysOf<int> gFFTBitTable;
+static const size_t MaxFastBits = 16;
+
+/* Declare Static functions */
+static void InitFFT();
+
+static bool IsPowerOfTwo(size_t x)
+{
+   if (x < 2)
+      return false;
+
+   if (x & (x - 1))             /* Thanks to 'byang' for this cute trick! */
+      return false;
+
+   return true;
+}
+
+static size_t NumberOfBitsNeeded(size_t PowerOfTwo)
+{
+   if (PowerOfTwo < 2) {
+      wxFprintf(stderr, "Error: FFT called with size %ld\n", PowerOfTwo);
+      exit(1);
+   }
+
+   size_t i = 0;
+   while (PowerOfTwo > 1)
+      PowerOfTwo >>= 1, ++i;
+
+   return i;
+}
+
+int ReverseBits(size_t index, size_t NumBits)
+{
+   size_t i, rev;
+
+   for (i = rev = 0; i < NumBits; i++) {
+      rev = (rev << 1) | (index & 1);
+      index >>= 1;
+   }
+
+   return rev;
+}
+
+void InitFFT()
+{
+   gFFTBitTable.reinit(MaxFastBits);
+
+   size_t len = 2;
+   for (size_t b = 1; b <= MaxFastBits; b++) {
+      auto &array = gFFTBitTable[b - 1];
+      array.reinit(len);
+      for (size_t i = 0; i < len; i++)
+         array[i] = ReverseBits(i, b);
+
+      len <<= 1;
+   }
+}
+
+void DeinitFFT()
+{
+   gFFTBitTable.reset();
+}
+
+static inline size_t FastReverseBits(size_t i, size_t NumBits)
+{
+   if (NumBits <= MaxFastBits)
+      return gFFTBitTable[NumBits - 1][i];
+   else
+      return ReverseBits(i, NumBits);
+}
+
+/*
+ * Complex Fast Fourier Transform
+ */
+
+void FFT(size_t NumSamples,
+         bool InverseTransform,
+         const float *RealIn, const float *ImagIn,
+	 float *RealOut, float *ImagOut)
+{
+   double angle_numerator = 2.0 * M_PI;
+   double tr, ti;                /* temp real, temp imaginary */
+
+   if (!IsPowerOfTwo(NumSamples)) {
+      wxFprintf(stderr, "%ld is not a power of two\n", NumSamples);
+      exit(1);
+   }
+
+   if (!gFFTBitTable)
+      InitFFT();
+
+   if (!InverseTransform)
+      angle_numerator = -angle_numerator;
+
+   /* Number of bits needed to store indices */
+   auto NumBits = NumberOfBitsNeeded(NumSamples);
+
+   /*
+    **   Do simultaneous data copy and bit-reversal ordering into outputs...
+    */
+
+   for (size_t i = 0; i < NumSamples; i++) {
+      auto j = FastReverseBits(i, NumBits);
+      RealOut[j] = RealIn[i];
+      ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
+   }
+
+   /*
+    **   Do the FFT itself...
+    */
+
+   size_t BlockEnd = 1;
+   for (size_t BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {
+
+      double delta_angle = angle_numerator / (double) BlockSize;
+
+      double sm2 = sin(-2 * delta_angle);
+      double sm1 = sin(-delta_angle);
+      double cm2 = cos(-2 * delta_angle);
+      double cm1 = cos(-delta_angle);
+      double w = 2 * cm1;
+      double ar0, ar1, ar2, ai0, ai1, ai2;
+
+      for (size_t i = 0; i < NumSamples; i += BlockSize) {
+         ar2 = cm2;
+         ar1 = cm1;
+
+         ai2 = sm2;
+         ai1 = sm1;
+
+         for (size_t j = i, n = 0; n < BlockEnd; j++, n++) {
+            ar0 = w * ar1 - ar2;
+            ar2 = ar1;
+            ar1 = ar0;
+
+            ai0 = w * ai1 - ai2;
+            ai2 = ai1;
+            ai1 = ai0;
+
+            size_t k = j + BlockEnd;
+            tr = ar0 * RealOut[k] - ai0 * ImagOut[k];
+            ti = ar0 * ImagOut[k] + ai0 * RealOut[k];
+
+            RealOut[k] = RealOut[j] - tr;
+            ImagOut[k] = ImagOut[j] - ti;
+
+            RealOut[j] += tr;
+            ImagOut[j] += ti;
+         }
+      }
+
+      BlockEnd = BlockSize;
+   }
+
+   /*
+      **   Need to normalize if inverse transform...
+    */
+
+   if (InverseTransform) {
+      float denom = (float) NumSamples;
+
+      for (size_t i = 0; i < NumSamples; i++) {
+         RealOut[i] /= denom;
+         ImagOut[i] /= denom;
+      }
+   }
+}
+
+/*
+ * Real Fast Fourier Transform
+ *
+ * This is merely a wrapper of RealFFTf() from RealFFTf.h.
+ */
+
+void RealFFT(size_t NumSamples, const float *RealIn, float *RealOut, float *ImagOut)
+{
+   auto hFFT = GetFFT(NumSamples);
+   Floats pFFT{ NumSamples };
+   // Copy the data into the processing buffer
+   for(size_t i = 0; i < NumSamples; i++)
+      pFFT[i] = RealIn[i];
+
+   // Perform the FFT
+   RealFFTf(pFFT.get(), hFFT.get());
+
+   // Copy the data into the real and imaginary outputs
+   for (size_t i = 1; i<(NumSamples / 2); i++) {
+      RealOut[i]=pFFT[hFFT->BitReversed[i]  ];
+      ImagOut[i]=pFFT[hFFT->BitReversed[i]+1];
+   }
+   // Handle the (real-only) DC and Fs/2 bins
+   RealOut[0] = pFFT[0];
+   RealOut[NumSamples / 2] = pFFT[1];
+   ImagOut[0] = ImagOut[NumSamples / 2] = 0;
+   // Fill in the upper half using symmetry properties
+   for(size_t i = NumSamples / 2 + 1; i < NumSamples; i++) {
+      RealOut[i] =  RealOut[NumSamples-i];
+      ImagOut[i] = -ImagOut[NumSamples-i];
+   }
+}
+
+/*
+ * InverseRealFFT
+ *
+ * This function computes the inverse of RealFFT, above.
+ * The RealIn and ImagIn is assumed to be conjugate-symmetric
+ * and as a result the output is purely real.
+ * Only the first half of RealIn and ImagIn are used due to this
+ * symmetry assumption.
+ *
+ * This is merely a wrapper of InverseRealFFTf() from RealFFTf.h.
+ */
+void InverseRealFFT(size_t NumSamples, const float *RealIn, const float *ImagIn,
+		    float *RealOut)
+{
+   auto hFFT = GetFFT(NumSamples);
+   Floats pFFT{ NumSamples };
+   // Copy the data into the processing buffer
+   for (size_t i = 0; i < (NumSamples / 2); i++)
+      pFFT[2*i  ] = RealIn[i];
+   if(ImagIn == NULL) {
+      for (size_t i = 0; i < (NumSamples / 2); i++)
+         pFFT[2*i+1] = 0;
+   } else {
+      for (size_t i = 0; i < (NumSamples / 2); i++)
+         pFFT[2*i+1] = ImagIn[i];
+   }
+   // Put the fs/2 component in the imaginary part of the DC bin
+   pFFT[1] = RealIn[NumSamples / 2];
+
+   // Perform the FFT
+   InverseRealFFTf(pFFT.get(), hFFT.get());
+
+   // Copy the data to the (purely real) output buffer
+   ReorderToTime(hFFT.get(), pFFT.get(), RealOut);
+}
+
+/*
+ * PowerSpectrum
+ *
+ * This function uses RealFFTf() from RealFFTf.h to perform the real
+ * FFT computation, and then squares the real and imaginary part of
+ * each coefficient, extracting the power and throwing away the phase.
+ *
+ * For speed, it does not call RealFFT, but duplicates some
+ * of its code.
+ */
+
+void PowerSpectrum(size_t NumSamples, const float *In, float *Out)
+{
+   auto hFFT = GetFFT(NumSamples);
+   Floats pFFT{ NumSamples };
+   // Copy the data into the processing buffer
+   for (size_t i = 0; i<NumSamples; i++)
+      pFFT[i] = In[i];
+
+   // Perform the FFT
+   RealFFTf(pFFT.get(), hFFT.get());
+
+   // Copy the data into the real and imaginary outputs
+   for (size_t i = 1; i<NumSamples / 2; i++) {
+      Out[i]= (pFFT[hFFT->BitReversed[i]  ]*pFFT[hFFT->BitReversed[i]  ])
+         + (pFFT[hFFT->BitReversed[i]+1]*pFFT[hFFT->BitReversed[i]+1]);
+   }
+   // Handle the (real-only) DC and Fs/2 bins
+   Out[0] = pFFT[0]*pFFT[0];
+   Out[NumSamples / 2] = pFFT[1]*pFFT[1];
+}
+
+/*
+ * Windowing Functions
+ */
+
+int NumWindowFuncs()
+{
+   return eWinFuncCount;
+}
+
+const TranslatableString WindowFuncName(int whichFunction)
+{
+   switch (whichFunction) {
+   default:
+   case eWinFuncRectangular:
+      return XO("Rectangular");
+   case eWinFuncBartlett:
+      /* i18n-hint a proper name */
+      return XO("Bartlett");
+   case eWinFuncHamming:
+      /* i18n-hint a proper name */
+      return XO("Hamming");
+   case eWinFuncHann:
+      /* i18n-hint a proper name */
+      return XO("Hann");
+   case eWinFuncBlackman:
+      /* i18n-hint a proper name */
+      return XO("Blackman");
+   case eWinFuncBlackmanHarris:
+      /* i18n-hint two proper names */
+      return XO("Blackman-Harris");
+   case eWinFuncWelch:
+      /* i18n-hint a proper name */
+      return XO("Welch");
+   case eWinFuncGaussian25:
+      /* i18n-hint a mathematical function named for C. F. Gauss */
+      return XO("Gaussian(a=2.5)");
+   case eWinFuncGaussian35:
+      /* i18n-hint a mathematical function named for C. F. Gauss */
+      return XO("Gaussian(a=3.5)");
+   case eWinFuncGaussian45:
+      /* i18n-hint a mathematical function named for C. F. Gauss */
+      return XO("Gaussian(a=4.5)");
+   }
+}
+
+void NewWindowFunc(int whichFunction, size_t NumSamplesIn, bool extraSample, float *in)
+{
+   int NumSamples = (int)NumSamplesIn;
+   if (extraSample) {
+      wxASSERT(NumSamples > 0);
+      --NumSamples;
+   }
+   wxASSERT(NumSamples > 0);
+
+   switch (whichFunction) {
+   default:
+      wxFprintf(stderr, "FFT::WindowFunc - Invalid window function: %d\n", whichFunction);
+      break;
+   case eWinFuncRectangular:
+      // Multiply all by 1.0f -- do nothing
+      break;
+
+   case eWinFuncBartlett:
+   {
+      // Bartlett (triangular) window
+      const int nPairs = (NumSamples - 1) / 2; // whether even or odd NumSamples, this is correct
+      const float denom = NumSamples / 2.0f;
+      in[0] = 0.0f;
+      for (int ii = 1;
+           ii <= nPairs; // Yes, <=
+           ++ii) {
+         const float value = ii / denom;
+         in[ii] *= value;
+         in[NumSamples - ii] *= value;
+      }
+      // When NumSamples is even, in[half] should be multiplied by 1.0, so unchanged
+      // When odd, the value of 1.0 is not reached
+   }
+      break;
+   case eWinFuncHamming:
+   {
+      // Hamming
+      const double multiplier = 2 * M_PI / NumSamples;
+      static const double coeff0 = 0.54, coeff1 = -0.46;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= coeff0 + coeff1 * cos(ii * multiplier);
+   }
+      break;
+   case eWinFuncHann:
+   {
+      // Hann
+      const double multiplier = 2 * M_PI / NumSamples;
+      static const double coeff0 = 0.5, coeff1 = -0.5;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= coeff0 + coeff1 * cos(ii * multiplier);
+   }
+      break;
+   case eWinFuncBlackman:
+   {
+      // Blackman
+      const double multiplier = 2 * M_PI / NumSamples;
+      const double multiplier2 = 2 * multiplier;
+      static const double coeff0 = 0.42, coeff1 = -0.5, coeff2 = 0.08;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= coeff0 + coeff1 * cos(ii * multiplier) + coeff2 * cos(ii * multiplier2);
+   }
+      break;
+   case eWinFuncBlackmanHarris:
+   {
+      // Blackman-Harris
+      const double multiplier = 2 * M_PI / NumSamples;
+      const double multiplier2 = 2 * multiplier;
+      const double multiplier3 = 3 * multiplier;
+      static const double coeff0 = 0.35875, coeff1 = -0.48829, coeff2 = 0.14128, coeff3 = -0.01168;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= coeff0 + coeff1 * cos(ii * multiplier) + coeff2 * cos(ii * multiplier2) + coeff3 * cos(ii * multiplier3);
+   }
+      break;
+   case eWinFuncWelch:
+   {
+      // Welch
+      const float N = NumSamples;
+      for (int ii = 0; ii < NumSamples; ++ii) {
+         const float iOverN = ii / N;
+         in[ii] *= 4 * iOverN * (1 - iOverN);
+      }
+   }
+      break;
+   case eWinFuncGaussian25:
+   {
+      // Gaussian (a=2.5)
+      // Precalculate some values, and simplify the fmla to try and reduce overhead
+      static const double A = -2 * 2.5*2.5;
+      const float N = NumSamples;
+      for (int ii = 0; ii < NumSamples; ++ii) {
+         const float iOverN = ii / N;
+         // full
+         // in[ii] *= exp(-0.5*(A*((ii-NumSamples/2)/NumSamples/2))*(A*((ii-NumSamples/2)/NumSamples/2)));
+         // reduced
+         in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN));
+      }
+   }
+      break;
+   case eWinFuncGaussian35:
+   {
+      // Gaussian (a=3.5)
+      static const double A = -2 * 3.5*3.5;
+      const float N = NumSamples;
+      for (int ii = 0; ii < NumSamples; ++ii) {
+         const float iOverN = ii / N;
+         in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN));
+      }
+   }
+      break;
+   case eWinFuncGaussian45:
+   {
+      // Gaussian (a=4.5)
+      static const double A = -2 * 4.5*4.5;
+      const float N = NumSamples;
+      for (int ii = 0; ii < NumSamples; ++ii) {
+         const float iOverN = ii / N;
+         in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN));
+      }
+   }
+      break;
+   }
+
+   if (extraSample && whichFunction != eWinFuncRectangular) {
+      double value = 0.0;
+      switch (whichFunction) {
+      case eWinFuncHamming:
+         value = 0.08;
+         break;
+      case eWinFuncGaussian25:
+         value = exp(-2 * 2.5 * 2.5 * 0.25);
+         break;
+      case eWinFuncGaussian35:
+         value = exp(-2 * 3.5 * 3.5 * 0.25);
+         break;
+      case eWinFuncGaussian45:
+         value = exp(-2 * 4.5 * 4.5 * 0.25);
+         break;
+      default:
+         break;
+      }
+      in[NumSamples] *= value;
+   }
+}
+
+// See cautions in FFT.h !
+void WindowFunc(int whichFunction, size_t NumSamples, float *in)
+{
+   bool extraSample = false;
+   switch (whichFunction)
+   {
+   case eWinFuncHamming:
+   case eWinFuncHann:
+   case eWinFuncBlackman:
+   case eWinFuncBlackmanHarris:
+      extraSample = true;
+      break;
+   default:
+      break;
+   case eWinFuncBartlett:
+      // PRL:  Do nothing here either
+      // But I want to comment that the old function did this case
+      // wrong in the second half of the array, in case NumSamples was odd
+      // but I think that never happened, so I am not bothering to preserve that
+      break;
+   }
+   NewWindowFunc(whichFunction, NumSamples, extraSample, in);
+}
+
+void DerivativeOfWindowFunc(int whichFunction, size_t NumSamples, bool extraSample, float *in)
+{
+   if (eWinFuncRectangular == whichFunction)
+   {
+      // Rectangular
+      // There are deltas at the ends
+      wxASSERT(NumSamples > 0);
+      --NumSamples;
+      // in[0] *= 1.0f;
+      for (int ii = 1; ii < (int)NumSamples; ++ii)
+         in[ii] = 0.0f;
+      in[NumSamples] *= -1.0f;
+      return;
+   }
+
+   if (extraSample) {
+      wxASSERT(NumSamples > 0);
+      --NumSamples;
+   }
+
+   wxASSERT(NumSamples > 0);
+
+   double A;
+   switch (whichFunction) {
+   case eWinFuncBartlett:
+   {
+      // Bartlett (triangular) window
+      // There are discontinuities in the derivative at the ends, and maybe at the midpoint
+      const int nPairs = (NumSamples - 1) / 2; // whether even or odd NumSamples, this is correct
+      const float value = 2.0f / NumSamples;
+      in[0] *=
+         // Average the two limiting values of discontinuous derivative
+         value / 2.0f;
+      for (int ii = 1;
+         ii <= nPairs; // Yes, <=
+         ++ii) {
+         in[ii] *= value;
+         in[NumSamples - ii] *= -value;
+      }
+      if (NumSamples % 2 == 0)
+         // Average the two limiting values of discontinuous derivative
+         in[NumSamples / 2] = 0.0f;
+      if (extraSample)
+         in[NumSamples] *=
+         // Average the two limiting values of discontinuous derivative
+            -value / 2.0f;
+      else
+         // Halve the multiplier previously applied
+         // Average the two limiting values of discontinuous derivative
+         in[NumSamples - 1] *= 0.5f;
+   }
+      break;
+   case eWinFuncHamming:
+   {
+      // Hamming
+      // There are deltas at the ends
+      const double multiplier = 2 * M_PI / NumSamples;
+      static const double coeff0 = 0.54, coeff1 = -0.46 * multiplier;
+      // TODO This code should be more explicit about the precision it intends.
+      // For now we get C4305 warnings, truncation from 'const double' to 'float' 
+      in[0] *= coeff0;
+      if (!extraSample)
+         --NumSamples;
+      for (int ii = 0; ii < (int)NumSamples; ++ii)
+         in[ii] *= - coeff1 * sin(ii * multiplier);
+      if (extraSample)
+         in[NumSamples] *= - coeff0;
+      else
+         // slightly different
+         in[NumSamples] *= - coeff0 - coeff1 * sin(NumSamples * multiplier);
+   }
+      break;
+   case eWinFuncHann:
+   {
+      // Hann
+      const double multiplier = 2 * M_PI / NumSamples;
+      const double coeff1 = -0.5 * multiplier;
+      for (int ii = 0; ii < (int)NumSamples; ++ii)
+         in[ii] *= - coeff1 * sin(ii * multiplier);
+      if (extraSample)
+         in[NumSamples] = 0.0f;
+   }
+      break;
+   case eWinFuncBlackman:
+   {
+      // Blackman
+      const double multiplier = 2 * M_PI / NumSamples;
+      const double multiplier2 = 2 * multiplier;
+      const double coeff1 = -0.5 * multiplier, coeff2 = 0.08 * multiplier2;
+      for (int ii = 0; ii < (int)NumSamples; ++ii)
+         in[ii] *= - coeff1 * sin(ii * multiplier) - coeff2 * sin(ii * multiplier2);
+      if (extraSample)
+         in[NumSamples] = 0.0f;
+   }
+      break;
+   case eWinFuncBlackmanHarris:
+   {
+      // Blackman-Harris
+      const double multiplier = 2 * M_PI / NumSamples;
+      const double multiplier2 = 2 * multiplier;
+      const double multiplier3 = 3 * multiplier;
+      const double coeff1 = -0.48829 * multiplier,
+         coeff2 = 0.14128 * multiplier2, coeff3 = -0.01168 * multiplier3;
+      for (int ii = 0; ii < (int)NumSamples; ++ii)
+         in[ii] *= - coeff1 * sin(ii * multiplier) - coeff2 * sin(ii * multiplier2) - coeff3 * sin(ii * multiplier3);
+      if (extraSample)
+         in[NumSamples] = 0.0f;
+   }
+      break;
+   case eWinFuncWelch:
+   {
+      // Welch
+      const float N = NumSamples;
+      const float NN = NumSamples * NumSamples;
+      for (int ii = 0; ii < (int)NumSamples; ++ii) {
+         in[ii] *= 4 * (N - ii - ii) / NN;
+      }
+      if (extraSample)
+         in[NumSamples] = 0.0f;
+      // Average the two limiting values of discontinuous derivative
+      in[0] /= 2.0f;
+      in[NumSamples - 1] /= 2.0f;
+   }
+      break;
+   case eWinFuncGaussian25:
+      // Gaussian (a=2.5)
+      A = -2 * 2.5*2.5;
+      goto Gaussian;
+   case eWinFuncGaussian35:
+      // Gaussian (a=3.5)
+      A = -2 * 3.5*3.5;
+      goto Gaussian;
+   case eWinFuncGaussian45:
+      // Gaussian (a=4.5)
+      A = -2 * 4.5*4.5;
+      goto Gaussian;
+   Gaussian:
+   {
+      // Gaussian (a=2.5)
+      // There are deltas at the ends
+      const float invN = 1.0f / NumSamples;
+      const float invNN = invN * invN;
+      // Simplify formula from the loop for ii == 0, add term for the delta
+      in[0] *= exp(A * 0.25) * (1 - invN);
+      if (!extraSample)
+         --NumSamples;
+      for (int ii = 1; ii < (int)NumSamples; ++ii) {
+         const float iOverN = ii * invN;
+         in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN)) * (2 * ii * invNN - invN);
+      }
+      if (extraSample)
+         in[NumSamples] *= exp(A * 0.25) * (invN - 1);
+      else {
+         // Slightly different
+         const float iOverN = NumSamples * invN;
+         in[NumSamples] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN)) * (2 * NumSamples * invNN - invN - 1);
+      }
+   }
+      break;
+   default:
+      wxFprintf(stderr, "FFT::DerivativeOfWindowFunc - Invalid window function: %d\n", whichFunction);
+   }
+}

libraries/lib-math/FFT.h

@@ -0,0 +1,160 @@
+/**********************************************************************
+
+  FFT.h
+
+  Dominic Mazzoni
+
+  September 2000
+
+  This file contains a few FFT routines, including a real-FFT
+  routine that is almost twice as fast as a normal complex FFT,
+  and a power spectrum routine which is more convenient when
+  you know you don't care about phase information.  It now also
+  contains a few basic windowing functions.
+
+  Some of this code was based on a free implementation of an FFT
+  by Don Cross, available on the web at:
+
+    http://www.intersrv.com/~dcross/fft.html
+
+  The basic algorithm for his code was based on Numerical Recipes
+  in Fortran.  I optimized his code further by reducing array
+  accesses, caching the bit reversal table, and eliminating
+  float-to-double conversions, and I added the routines to
+  calculate a real FFT and a real power spectrum.
+
+  Note: all of these routines use single-precision floats.
+  I have found that in practice, floats work well until you
+  get above 8192 samples.  If you need to do a larger FFT,
+  you need to use doubles.
+
+**********************************************************************/
+#ifndef __AUDACITY_FFT_H__
+#define __AUDACITY_FFT_H__
+
+#include <wx/defs.h>
+
+class TranslatableString;
+
+/*
+  Salvo Ventura - November 2006
+  Added more window functions:
+    * 4: Blackman
+    * 5: Blackman-Harris
+    * 6: Welch
+    * 7: Gaussian(a=2.5)
+    * 8: Gaussian(a=3.5)
+    * 9: Gaussian(a=4.5)
+*/
+
+#include <wx/defs.h>
+#include <wx/wxchar.h>
+
+#ifndef M_PI
+#define	M_PI		3.14159265358979323846  /* pi */
+#endif
+
+/*
+ * This is the function you will use the most often.
+ * Given an array of floats, this will compute the power
+ * spectrum by doing a Real FFT and then computing the
+ * sum of the squares of the real and imaginary parts.
+ * Note that the output array is half the length of the
+ * input array, and that NumSamples must be a power of two.
+ */
+
+MATH_API
+void PowerSpectrum(size_t NumSamples, const float *In, float *Out);
+
+/*
+ * Computes an FFT when the input data is real but you still
+ * want complex data as output.  The output arrays are the
+ * same length as the input, but will be conjugate-symmetric
+ * NumSamples must be a power of two.
+ */
+
+MATH_API
+void RealFFT(size_t NumSamples,
+             const float *RealIn, float *RealOut, float *ImagOut);
+
+/*
+ * Computes an Inverse FFT when the input data is conjugate symmetric
+ * so the output is purely real.  NumSamples must be a power of
+ * two.
+ */
+MATH_API
+void InverseRealFFT(size_t NumSamples,
+		    const float *RealIn, const float *ImagIn, float *RealOut);
+
+/*
+ * Computes a FFT of complex input and returns complex output.
+ * Currently this is the only function here that supports the
+ * inverse transform as well.
+ */
+
+MATH_API
+void FFT(size_t NumSamples,
+         bool InverseTransform,
+         const float *RealIn, const float *ImagIn, float *RealOut, float *ImagOut);
+
+/*
+ * Multiply values in data by values of the chosen function
+ * DO NOT REUSE!  Prefer NewWindowFunc instead
+ * This version was inconsistent whether the window functions were
+ * symmetrical about NumSamples / 2, or about (NumSamples - 1) / 2
+ * It remains for compatibility until we decide to upgrade all the old uses
+ * All functions have 0 in data[0] except Rectangular, Hamming and Gaussians
+ */
+
+enum eWindowFunctions
+{
+   eWinFuncRectangular,
+   eWinFuncBartlett,
+   eWinFuncHamming,
+   eWinFuncHann,
+   eWinFuncBlackman,
+   eWinFuncBlackmanHarris,
+   eWinFuncWelch,
+   eWinFuncGaussian25,
+   eWinFuncGaussian35,
+   eWinFuncGaussian45,
+   eWinFuncCount
+};
+
+MATH_API
+void WindowFunc(int whichFunction, size_t NumSamples, float *data);
+
+/*
+ * Multiply values in data by values of the chosen function
+ * All functions are symmetrical about NumSamples / 2 if extraSample is false,
+ * otherwise about (NumSamples - 1) / 2
+ * All functions have 0 in data[0] except Rectangular, Hamming and Gaussians
+ */
+MATH_API
+void NewWindowFunc(int whichFunction, size_t NumSamples, bool extraSample, float *data);
+
+/*
+ * Multiply values in data by derivative of the chosen function, assuming
+ * sampling interval is unit
+ * All functions are symmetrical about NumSamples / 2 if extraSample is false,
+ * otherwise about (NumSamples - 1) / 2
+ * All functions have 0 in data[0] except Rectangular, Hamming and Gaussians
+ */
+MATH_API
+void DerivativeOfWindowFunc(int whichFunction, size_t NumSamples, bool extraSample, float *data);
+
+/*
+ * Returns the name of the windowing function (for UI display)
+ */
+
+MATH_API const TranslatableString WindowFuncName(int whichFunction);
+
+/*
+ * Returns the number of windowing functions supported
+ */
+
+MATH_API int NumWindowFuncs();
+
+MATH_API void DeinitFFT();
+
+#endif

libraries/lib-math/InterpolateAudio.cpp

@@ -0,0 +1,204 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  InterpolateAudio.cpp
+
+  Dominic Mazzoni
+
+**********************************************************************/
+
+#include "InterpolateAudio.h"
+
+#include <math.h>
+#include <stdlib.h>
+
+#include <wx/defs.h>
+
+#include "Matrix.h"
+
+static inline int imin(int x, int y)
+{
+   return x<y? x: y;
+}
+
+static inline int imax(int x, int y)
+{
+   return x>y? x: y;
+}
+
+// This function is a really dumb, simple way to interpolate audio,
+// if the more general InterpolateAudio function below doesn't have
+// enough data to work with.  If the bad samples are in the middle,
+// it's literally linear.  If it's on either edge, we add some decay
+// back to zero.
+static void LinearInterpolateAudio(float *buffer, int len,
+                                   int firstBad, int numBad)
+{
+   int i;
+
+   float decay = 0.9f;
+
+   if (firstBad==0) {
+      float delta = buffer[numBad] - buffer[numBad+1];
+      float value = buffer[numBad];
+      i = numBad - 1;
+      while (i >= 0) {
+         value += delta;
+         buffer[i] = value;
+         value *= decay;
+         delta *= decay;
+         i--;
+      }
+   }
+   else if (firstBad + numBad == len) {
+      float delta = buffer[firstBad-1] - buffer[firstBad-2];
+      float value = buffer[firstBad-1];
+      i = firstBad;
+      while (i < firstBad + numBad) {
+         value += delta;
+         buffer[i] = value;
+         value *= decay;
+         delta *= decay;
+         i++;
+      }
+   }
+   else {
+      float v1 = buffer[firstBad-1];
+      float v2 = buffer[firstBad+numBad];
+      float value = v1;
+      float delta = (v2 - v1) / (numBad+1);
+      i = firstBad;
+      while (i < firstBad + numBad) {
+         value += delta;
+         buffer[i] = value;
+         i++;
+      }
+   }
+}
+
+// Here's the main interpolate function, using
+// Least Squares AutoRegression (LSAR):
+void InterpolateAudio(float *buffer, const size_t len,
+                      size_t firstBad, size_t numBad)
+{
+   const auto N = len;
+
+   wxASSERT(len > 0 &&
+            firstBad >= 0 &&
+            numBad < len &&
+            firstBad+numBad <= len);
+
+   if(numBad >= len)
+      return;  //should never have been called!
+
+   if (firstBad == 0) {
+      // The algorithm below has a weird asymmetry in that it
+      // performs poorly when interpolating to the left.  If
+      // we're asked to interpolate the left side of a buffer,
+      // we just reverse the problem and try it that way.
+      Floats buffer2{ len };
+      for(size_t i=0; i<len; i++)
+         buffer2[len-1-i] = buffer[i];
+      InterpolateAudio(buffer2.get(), len, len-numBad, numBad);
+      for(size_t i=0; i<len; i++)
+         buffer[len-1-i] = buffer2[i];
+      return;
+   }
+
+   Vector s(len, buffer);
+
+   // Choose P, the order of the autoregression equation
+   const int IP =
+      imin(imin(numBad * 3, 50), imax(firstBad - 1, len - (firstBad + numBad) - 1));
+
+   if (IP < 3 || IP >= (int)N) {
+      LinearInterpolateAudio(buffer, len, firstBad, numBad);
+      return;
+   }
+
+   size_t P(IP);
+
+   // Add a tiny amount of random noise to the input signal -
+   // this sounds like a bad idea, but the amount we're adding
+   // is only about 1 bit in 16-bit audio, and it's an extremely
+   // effective way to avoid nearly-singular matrices.  If users
+   // run it more than once they get slightly different results;
+   // this is sometimes even advantageous.
+   for(size_t i=0; i<N; i++)
+      s[i] += (rand()-(RAND_MAX/2))/(RAND_MAX*10000.0);
+
+   // Solve for the best autoregression coefficients
+   // using a least-squares fit to all of the non-bad
+   // data we have in the buffer
+   Matrix X(P, P);
+   Vector b(P);
+
+   for(size_t i = 0; i + P < len; i++)
+      if (i+P < firstBad || i >= (firstBad + numBad))
+         for(size_t row=0; row<P; row++) {
+            for(size_t col=0; col<P; col++)
+               X[row][col] += (s[i+row] * s[i+col]);
+            b[row] += s[i+P] * s[i+row];
+         }
+
+   Matrix Xinv(P, P);
+   if (!InvertMatrix(X, Xinv)) {
+      // The matrix is singular!  Fall back on linear...
+      // In practice I have never seen this happen if
+      // we add the tiny bit of random noise.
+      LinearInterpolateAudio(buffer, len, firstBad, numBad);
+      return;
+   }
+
+   // This vector now contains the autoregression coefficients
+   const Vector &a = Xinv * b;
+
+   // Create a matrix (a "Toeplitz" matrix, as it turns out)
+   // which encodes the autoregressive relationship between
+   // elements of the sequence.
+   Matrix A(N-P, N);
+   for(size_t row=0; row<N-P; row++) {
+      for(size_t col=0; col<P; col++)
+         A[row][row+col] = -a[col];
+      A[row][row+P] = 1;
+   }
+
+   // Split both the Toeplitz matrix and the signal into
+   // two pieces.  Note that this code could be made to
+   // work even in the case where the "bad" samples are
+   // not contiguous, but currently it assumes they are.
+   //   "u" is for unknown (bad)
+   //   "k" is for known (good)
+   Matrix Au = MatrixSubset(A, 0, N-P, firstBad, numBad);
+   Matrix A_left = MatrixSubset(A, 0, N-P, 0, firstBad);
+   Matrix A_right = MatrixSubset(A, 0, N-P,
+                                 firstBad+numBad, N-(firstBad+numBad));
+   Matrix Ak = MatrixConcatenateCols(A_left, A_right);
+
+   const Vector &s_left = VectorSubset(s, 0, firstBad);
+   const Vector &s_right = VectorSubset(s, firstBad+numBad,
+                                 N-(firstBad+numBad));
+   const Vector &sk = VectorConcatenate(s_left, s_right);
+
+   // Do some linear algebra to find the best possible
+   // values that fill in the "bad" area
+   Matrix AuT = TransposeMatrix(Au);
+   Matrix X1 = MatrixMultiply(AuT, Au);
+   Matrix X2(X1.Rows(), X1.Cols());
+   if (!InvertMatrix(X1, X2)) {
+      // The matrix is singular!  Fall back on linear...
+      LinearInterpolateAudio(buffer, len, firstBad, numBad);
+      return;
+   }
+   Matrix X2b = X2 * -1.0;
+   Matrix X3 = MatrixMultiply(X2b, AuT);
+   Matrix X4 = MatrixMultiply(X3, Ak);
+   // This vector contains our best guess as to the
+   // unknown values
+   const Vector &su = X4 * sk;
+
+   // Put the results into the return buffer
+   for(size_t i=0; i<numBad; i++)
+      buffer[firstBad+i] = (float)su[i];
+}

libraries/lib-math/InterpolateAudio.h

@@ -0,0 +1,41 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  InterpolateAudio.h
+
+  Dominic Mazzoni
+
+*******************************************************************//*!
+
+\file Matrix.h
+\brief General routine to interpolate (or even extrapolate small amounts)
+ audio when a few of the samples are bad.  Works great for a few
+ dozen bad samples, but not so well with hundreds.  Uses the
+ least-squares autoregression (LSAR) algorithm, as described in:
+
+ Simon Godsill, Peter Rayner, and Olivier Cappe.  Digital Audio Restoration.
+ Berlin: Springer, 1998.
+
+ This is the same work used by Gnome Wave Cleaner (GWC), however this
+ implementation is original.
+
+*//*******************************************************************/
+
+#ifndef __AUDACITY_INTERPOLATE_AUDIO__
+#define __AUDACITY_INTERPOLATE_AUDIO__
+
+#include <cstddef>
+
+// See top of file for a description of the algorithm.  Interpolates
+// the samples from buffer[firstBad] through buffer[firstBad+numBad-1],
+// ignoring whatever value was there previously, and replacing them with
+// values determined from the surrounding audio.  Works best if the bad
+// samples are in the middle, with several times as much data on either
+// side (6x the number of bad samples on either side is great).  However,
+// it will work with less data, and with the bad samples on one end or
+// the other.
+void MATH_API InterpolateAudio(float *buffer, size_t len,
+                                       size_t firstBad, size_t numBad);
+
+#endif // __AUDACITY_INTERPOLATE_AUDIO__

libraries/lib-math/Matrix.cpp

@@ -0,0 +1,345 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  Matrix.cpp
+
+  Dominic Mazzoni
+
+**********************************************************************/
+
+#include "Matrix.h"
+
+#include <stdlib.h>
+#include <math.h>
+
+#include <wx/defs.h>
+
+Vector::Vector()
+{
+}
+
+Vector::Vector(unsigned len, double *data)
+   : mN{ len }
+   , mData(len)
+{
+   if (data)
+      std::copy(data, data + len, mData.get());
+   else
+      std::fill(mData.get(), mData.get() + len, 0.0);
+}
+
+Vector::Vector(unsigned len, float *data)
+   : mN{ len }
+   , mData{ len }
+{
+   if (data)
+      std::copy(data, data + len, mData.get());
+   else
+      std::fill(mData.get(), mData.get() + len, 0.0);
+}
+
+Vector& Vector::operator=(const Vector &other)
+{
+   wxASSERT(Len() == other.Len());
+   std::copy(other.mData.get(), other.mData.get() + mN, mData.get());
+   return *this;
+}
+
+Vector::Vector(const Vector &other)
+   : mN{ other.Len() }
+   , mData{ mN }
+{
+   std::copy(other.mData.get(), other.mData.get() + mN, mData.get());
+}
+
+Vector::~Vector()
+{
+}
+
+void Vector::Reinit(unsigned len)
+{
+   Vector temp(len);
+   Swap(temp);
+}
+
+void Vector::Swap(Vector &that)
+{
+   std::swap(mN, that.mN);
+   mData.swap(that.mData);
+}
+
+double Vector::Sum() const
+{
+   double sum = 0.0;
+   for(unsigned i = 0; i < Len(); i++)
+      sum += mData[i];
+   return sum;
+}
+
+Matrix::Matrix(unsigned rows, unsigned cols, double **data)
+   : mRows{ rows }
+   , mCols{ cols }
+   , mRowVec{ mRows }
+{
+   for(unsigned i = 0; i < mRows; i++) {
+      mRowVec[i].Reinit( mCols );
+      for(unsigned j = 0; j < mCols; j++) {
+         if (data)
+            (*this)[i][j] = data[i][j];
+         else
+            (*this)[i][j] = 0.0;
+      }
+   }
+}
+
+Matrix& Matrix::operator=(const Matrix &other)
+{
+   CopyFrom(other);
+   return *this;
+}
+
+Matrix::Matrix(const Matrix &other)
+{
+   CopyFrom(other);
+}
+
+void Matrix::CopyFrom(const Matrix &other)
+{
+   mRows = other.mRows;
+   mCols = other.mCols;
+   mRowVec.reinit(mRows);
+   for (unsigned i = 0; i < mRows; i++) {
+      mRowVec[i].Reinit( mCols );
+      mRowVec[i] = other.mRowVec[i];
+   }
+}
+
+Matrix::~Matrix()
+{
+}
+
+void Matrix::SwapRows(unsigned i, unsigned j)
+{
+   mRowVec[i].Swap(mRowVec[j]);
+}
+
+Matrix IdentityMatrix(unsigned N)
+{
+   Matrix M(N, N);
+   for(unsigned i = 0; i < N; i++)
+      M[i][i] = 1.0;
+   return M;
+}
+
+Vector operator+(const Vector &left, const Vector &right)
+{
+   wxASSERT(left.Len() == right.Len());
+   Vector v(left.Len());
+   for(unsigned i = 0; i < left.Len(); i++)
+      v[i] = left[i] + right[i];
+   return v;
+}
+
+Vector operator-(const Vector &left, const Vector &right)
+{
+   wxASSERT(left.Len() == right.Len());
+   Vector v(left.Len());
+   for(unsigned i = 0; i < left.Len(); i++)
+      v[i] = left[i] - right[i];
+   return v;
+}
+
+Vector operator*(const Vector &left, const Vector &right)
+{
+   wxASSERT(left.Len() == right.Len());
+   Vector v(left.Len());
+   for(unsigned i = 0; i < left.Len(); i++)
+      v[i] = left[i] * right[i];
+   return v;
+}
+
+Vector operator*(const Vector &left, double right)
+{
+   Vector v(left.Len());
+   for(unsigned i = 0; i < left.Len(); i++)
+      v[i] = left[i] * right;
+   return v;
+}
+
+Vector VectorSubset(const Vector &other, unsigned start, unsigned len)
+{
+   Vector v(len);
+   for(unsigned i = 0; i < len; i++)
+      v[i] = other[start+i];
+   return v;
+}
+
+Vector VectorConcatenate(const Vector& left, const Vector& right)
+{
+   Vector v(left.Len() + right.Len());
+   for(unsigned i = 0; i < left.Len(); i++)
+      v[i] = left[i];
+   for(unsigned i = 0; i < right.Len(); i++)
+      v[i + left.Len()] = right[i];
+   return v;
+}
+
+Vector operator*(const Vector &left, const Matrix &right)
+{
+   wxASSERT(left.Len() == right.Rows());
+   Vector v(right.Cols());
+   for(unsigned i = 0; i < right.Cols(); i++) {
+      v[i] = 0.0;
+      for(unsigned j = 0; j < right.Rows(); j++)
+         v[i] += left[j] * right[j][i];
+   }
+   return v;
+}
+
+Vector operator*(const Matrix &left, const Vector &right)
+{
+   wxASSERT(left.Cols() == right.Len());
+   Vector v(left.Rows());
+   for(unsigned i = 0; i < left.Rows(); i++) {
+      v[i] = 0.0;
+      for(unsigned j = 0; j < left.Cols(); j++)
+         v[i] += left[i][j] * right[j];
+   }
+   return v;
+}
+
+Matrix operator+(const Matrix &left, const Matrix &right)
+{
+   wxASSERT(left.Rows() == right.Rows());
+   wxASSERT(left.Cols() == right.Cols());
+   Matrix M(left.Rows(), left.Cols());
+   for(unsigned i = 0; i < left.Rows(); i++)
+      for(unsigned j = 0; j < left.Cols(); j++)
+         M[i][j] = left[i][j] + right[i][j];
+   return M;
+}
+
+Matrix operator*(const Matrix &left, const double right)
+{
+   Matrix M(left.Rows(), left.Cols());
+   for(unsigned i = 0; i < left.Rows(); i++)
+      for(unsigned j = 0; j < left.Cols(); j++)
+         M[i][j] = left[i][j] * right;
+   return M;
+}
+
+Matrix ScalarMultiply(const Matrix &left, const Matrix &right)
+{
+   wxASSERT(left.Rows() == right.Rows());
+   wxASSERT(left.Cols() == right.Cols());
+   Matrix M(left.Rows(), left.Cols());
+   for(unsigned i = 0; i < left.Rows(); i++)
+      for(unsigned j = 0; j < left.Cols(); j++)
+         M[i][j] = left[i][j] * right[i][j];
+   return M;
+}
+
+Matrix MatrixMultiply(const Matrix &left, const Matrix &right)
+{
+   wxASSERT(left.Cols() == right.Rows());
+   Matrix M(left.Rows(), right.Cols());
+   for(unsigned i = 0; i < left.Rows(); i++)
+      for(unsigned j = 0; j < right.Cols(); j++) {
+         M[i][j] = 0.0;
+         for(unsigned k = 0; k < left.Cols(); k++)
+            M[i][j] += left[i][k] * right[k][j];
+      }
+   return M;
+}
+
+Matrix MatrixSubset(const Matrix &input,
+                    unsigned startRow, unsigned numRows,
+                    unsigned startCol, unsigned numCols)
+{
+   Matrix M(numRows, numCols);
+   for(unsigned i = 0; i < numRows; i++)
+      for(unsigned j = 0; j < numCols; j++)
+         M[i][j] = input[startRow+i][startCol+j];
+   return M;
+}
+
+Matrix MatrixConcatenateCols(const Matrix& left, const Matrix& right)
+{
+   wxASSERT(left.Rows() == right.Rows());
+   Matrix M(left.Rows(), left.Cols() + right.Cols());
+   for(unsigned i = 0; i < left.Rows(); i++) {
+      for(unsigned j = 0; j < left.Cols(); j++)
+         M[i][j] = left[i][j];
+      for(unsigned j = 0; j < right.Cols(); j++)
+         M[i][j+left.Cols()] = right[i][j];
+   }
+   return M;
+}
+
+Matrix TransposeMatrix(const Matrix& other)
+{
+   Matrix M(other.Cols(), other.Rows());
+   for(unsigned i = 0; i < other.Rows(); i++)
+      for(unsigned j = 0; j < other.Cols(); j++)
+         M[j][i] = other[i][j];
+   return M;
+}
+
+bool InvertMatrix(const Matrix& input, Matrix& Minv)
+{
+   // Very straightforward implementation of
+   // Gauss-Jordan elimination to invert a matrix.
+   // Returns true if successful
+
+   wxASSERT(input.Rows() == input.Cols());
+   auto N = input.Rows();
+
+   Matrix M = input;
+   Minv = IdentityMatrix(N);
+
+   // Do the elimination one column at a time
+   for(unsigned i = 0; i < N; i++) {
+      // Pivot the row with the largest absolute value in
+      // column i, into row i
+      double absmax = 0.0;
+      unsigned int argmax = 0;
+
+      for(unsigned j = i; j < N; j++)
+         if (fabs(M[j][i]) > absmax) {
+            absmax = fabs(M[j][i]);
+            argmax = j;
+         }
+
+      // If no row has a nonzero value in that column,
+      // the matrix is singular and we have to give up.
+      if (absmax == 0)
+         return false;
+
+      if (i != argmax) {
+         M.SwapRows(i, argmax);
+         Minv.SwapRows(i, argmax);
+      }
+
+      // Divide this row by the value of M[i][i]
+      double factor = 1.0 / M[i][i];
+      M[i] = M[i] * factor;
+      Minv[i] = Minv[i] * factor;
+
+      // Eliminate the rest of the column
+      for(unsigned j = 0; j < N; j++) {
+         if (j == i)
+            continue;
+         if (fabs(M[j][i]) > 0) {
+            // Subtract a multiple of row i from row j
+            factor = M[j][i];
+            for(unsigned k = 0; k < N; k++) {
+               M[j][k] -= (M[i][k] * factor);
+               Minv[j][k] -= (Minv[i][k] * factor);
+            }
+         }
+      }
+   }
+
+   return true;
+}

libraries/lib-math/Matrix.h

@@ -0,0 +1,113 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  Matrix.h
+
+  Dominic Mazzoni
+
+*******************************************************************//*!
+
+\file Matrix.h
+\brief Holds both the Matrix and Vector classes, supporting
+  linear algebra operations, including matrix inversion.
+  Used by InterpolateAudio.
+
+\class Matrix
+\brief Holds a matrix of doubles and supports arithmetic, subsetting,
+  and matrix inversion.  Used by InterpolateAudio.
+
+\class Vector
+\brief Holds a matrix of doubles and supports arithmetic operations,
+  including Vector-Matrix operations.  Used by InterpolateAudio.
+
+*//*******************************************************************/
+
+#ifndef __AUDACITY_MATRIX__
+#define __AUDACITY_MATRIX__
+
+#include "SampleFormat.h"
+
+class Matrix;
+
+class Vector
+{
+ public:
+   Vector();
+   Vector(const Vector& copyFrom);
+   Vector(unsigned len, double *data=NULL);
+   Vector(unsigned len, float *data);
+   Vector& operator=(const Vector &other);
+   ~Vector();
+
+   void Reinit(unsigned len);
+   void Swap(Vector &that);
+
+   inline double& operator[](unsigned i) { return mData[i]; }
+   inline double operator[](unsigned i) const { return mData[i]; }
+   inline unsigned Len() const { return mN; }
+
+   double Sum() const;
+
+ private:
+   unsigned mN{ 0 };
+   Doubles mData;
+};
+
+class Matrix
+{
+ public:
+   Matrix(const Matrix& copyFrom);
+   Matrix(unsigned rows, unsigned cols, double **data=NULL);
+   ~Matrix();
+
+   Matrix& operator=(const Matrix& other);
+
+   inline Vector& operator[](unsigned i) { return mRowVec[i]; }
+   inline Vector& operator[](unsigned i) const { return mRowVec[i]; }
+   inline unsigned Rows() const { return mRows; }
+   inline unsigned Cols() const { return mCols; }
+
+   void SwapRows(unsigned i, unsigned j);
+
+ private:
+   void CopyFrom(const Matrix& other);
+
+   unsigned mRows;
+   unsigned mCols;
+   ArrayOf<Vector> mRowVec;
+};
+
+bool InvertMatrix(const Matrix& input, Matrix& Minv);
+
+Matrix TransposeMatrix(const Matrix& M);
+
+Matrix IdentityMatrix(unsigned N);
+
+Vector operator+(const Vector &left, const Vector &right);
+Vector operator-(const Vector &left, const Vector &right);
+Vector operator*(const Vector &left, const Vector &right);
+Vector operator*(const Vector &left, double right);
+
+Vector VectorSubset(const Vector &other, unsigned start, unsigned len);
+Vector VectorConcatenate(const Vector& left, const Vector& right);
+
+Vector operator*(const Vector &left, const Matrix &right);
+Vector operator*(const Matrix &left, const Vector &right);
+
+Matrix operator+(const Matrix &left, const Matrix &right);
+Matrix operator*(const Matrix &left, const double right);
+
+// No operator* on matrices due to ambiguity
+Matrix ScalarMultiply(const Matrix &left, const Matrix &right);
+Matrix MatrixMultiply(const Matrix &left, const Matrix &right);
+
+Matrix MatrixSubset(const Matrix &M,
+                    unsigned startRow, unsigned numRows,
+                    unsigned startCol, unsigned numCols);
+
+Matrix MatrixConcatenateCols(const Matrix& left, const Matrix& right);
+
+bool InvertMatrix(const Matrix& M, Matrix& Minv);
+
+#endif // __AUDACITY_MATRIX__

libraries/lib-math/RealFFTf.cpp

@@ -0,0 +1,368 @@
+/*
+*     Program: REALFFTF.C
+*      Author: Philip Van Baren
+*        Date: 2 September 1993
+*
+* Description: These routines perform an FFT on real data to get a conjugate-symmetric
+*              output, and an inverse FFT on conjugate-symmetric input to get a real
+*              output sequence.
+*
+*              This code is for floating point data.
+*
+*              Modified 8/19/1998 by Philip Van Baren
+*                 - made the InitializeFFT and EndFFT routines take a structure
+*                   holding the length and pointers to the BitReversed and SinTable
+*                   tables.
+*              Modified 5/23/2009 by Philip Van Baren
+*                 - Added GetFFT and ReleaseFFT routines to retain common SinTable
+*                   and BitReversed tables so they don't need to be reallocated
+*                   and recomputed on every call.
+*                 - Added Reorder* functions to undo the bit-reversal
+*
+*  Copyright (C) 2009  Philip VanBaren
+*
+*  This program is free software; you can redistribute it and/or modify
+*  it under the terms of the GNU General Public License as published by
+*  the Free Software Foundation; either version 2 of the License, or
+*  (at your option) any later version.
+*
+*  This program is distributed in the hope that it will be useful,
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*  GNU General Public License for more details.
+*
+*  You should have received a copy of the GNU General Public License
+*  along with this program; if not, write to the Free Software
+*  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#include "RealFFTf.h"
+
+#include <vector>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include <wx/thread.h>
+
+#ifndef M_PI
+#define	M_PI		3.14159265358979323846  /* pi */
+#endif
+
+/*
+*  Initialize the Sine table and Twiddle pointers (bit-reversed pointers)
+*  for the FFT routine.
+*/
+HFFT InitializeFFT(size_t fftlen)
+{
+   int temp;
+   HFFT h{ safenew FFTParam };
+
+   /*
+   *  FFT size is only half the number of data points
+   *  The full FFT output can be reconstructed from this FFT's output.
+   *  (This optimization can be made since the data is real.)
+   */
+   h->Points = fftlen / 2;
+
+   h->SinTable.reinit(2*h->Points);
+
+   h->BitReversed.reinit(h->Points);
+
+   for(size_t i = 0; i < h->Points; i++)
+   {
+      temp = 0;
+      for(size_t mask = h->Points / 2; mask > 0; mask >>= 1)
+         temp = (temp >> 1) + (i & mask ? h->Points : 0);
+
+      h->BitReversed[i] = temp;
+   }
+
+   for(size_t i = 0; i < h->Points; i++)
+   {
+      h->SinTable[h->BitReversed[i]  ]=(fft_type)-sin(2*M_PI*i/(2*h->Points));
+      h->SinTable[h->BitReversed[i]+1]=(fft_type)-cos(2*M_PI*i/(2*h->Points));
+   }
+
+#ifdef EXPERIMENTAL_EQ_SSE_THREADED
+   // NEW SSE FFT routines work on live data
+   for(size_t i = 0; i < 32; i++)
+      if((1 << i) & fftlen)
+         h->pow2Bits = i;
+#endif
+
+   return h;
+}
+
+enum : size_t { MAX_HFFT = 10 };
+
+// Maintain a pool:
+static std::vector< std::unique_ptr<FFTParam> > hFFTArray(MAX_HFFT);
+wxCriticalSection getFFTMutex;
+
+/* Get a handle to the FFT tables of the desired length */
+/* This version keeps common tables rather than allocating a NEW table every time */
+HFFT GetFFT(size_t fftlen)
+{
+   // To do:  smarter policy about when to retain in the pool and when to
+   // allocate a unique instance.
+
+   wxCriticalSectionLocker locker{ getFFTMutex };
+   
+   size_t h = 0;
+   auto n = fftlen/2;
+   auto size = hFFTArray.size();
+   for(;
+       (h < size) && hFFTArray[h] && (n != hFFTArray[h]->Points);
+       h++)
+      ;
+   if(h < size) {
+      if(hFFTArray[h] == NULL) {
+         hFFTArray[h].reset( InitializeFFT(fftlen).release() );
+      }
+      return HFFT{ hFFTArray[h].get() };
+   } else {
+      // All buffers used, so fall back to allocating a NEW set of tables
+      return InitializeFFT(fftlen);
+   }
+}
+
+/* Release a previously requested handle to the FFT tables */
+void FFTDeleter::operator() (FFTParam *hFFT) const
+{
+   wxCriticalSectionLocker locker{ getFFTMutex };
+
+   auto it = hFFTArray.begin(), end = hFFTArray.end();
+   while (it != end && it->get() != hFFT)
+      ++it;
+   if ( it != end )
+      ;
+   else
+      delete hFFT;
+}
+
+/*
+*  Forward FFT routine.  Must call GetFFT(fftlen) first!
+*
+*  Note: Output is BIT-REVERSED! so you must use the BitReversed to
+*        get legible output, (i.e. Real_i = buffer[ h->BitReversed[i] ]
+*                                  Imag_i = buffer[ h->BitReversed[i]+1 ] )
+*        Input is in normal order.
+*
+* Output buffer[0] is the DC bin, and output buffer[1] is the Fs/2 bin
+* - this can be done because both values will always be real only
+* - this allows us to not have to allocate an extra complex value for the Fs/2 bin
+*
+*  Note: The scaling on this is done according to the standard FFT definition,
+*        so a unit amplitude DC signal will output an amplitude of (N)
+*        (Older revisions would progressively scale the input, so the output
+*        values would be similar in amplitude to the input values, which is
+*        good when using fixed point arithmetic)
+*/
+void RealFFTf(fft_type *buffer, const FFTParam *h)
+{
+   fft_type *A,*B;
+   const fft_type *sptr;
+   const fft_type *endptr1,*endptr2;
+   const int *br1,*br2;
+   fft_type HRplus,HRminus,HIplus,HIminus;
+   fft_type v1,v2,sin,cos;
+
+   auto ButterfliesPerGroup = h->Points/2;
+
+   /*
+   *  Butterfly:
+   *     Ain-----Aout
+   *         \ /
+   *         / \
+   *     Bin-----Bout
+   */
+
+   endptr1 = buffer + h->Points * 2;
+
+   while(ButterfliesPerGroup > 0)
+   {
+      A = buffer;
+      B = buffer + ButterfliesPerGroup * 2;
+      sptr = h->SinTable.get();
+
+      while(A < endptr1)
+      {
+         sin = *sptr;
+         cos = *(sptr+1);
+         endptr2 = B;
+         while(A < endptr2)
+         {
+            v1 = *B * cos + *(B + 1) * sin;
+            v2 = *B * sin - *(B + 1) * cos;
+            *B = (*A + v1);
+            *(A++) = *(B++) - 2 * v1;
+            *B = (*A - v2);
+            *(A++) = *(B++) + 2 * v2;
+         }
+         A = B;
+         B += ButterfliesPerGroup * 2;
+         sptr += 2;
+      }
+      ButterfliesPerGroup >>= 1;
+   }
+   /* Massage output to get the output for a real input sequence. */
+   br1 = h->BitReversed.get() + 1;
+   br2 = h->BitReversed.get() + h->Points - 1;
+
+   while(br1<br2)
+   {
+      sin=h->SinTable[*br1];
+      cos=h->SinTable[*br1+1];
+      A=buffer+*br1;
+      B=buffer+*br2;
+      HRplus = (HRminus = *A     - *B    ) + (*B     * 2);
+      HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) * 2);
+      v1 = (sin*HRminus - cos*HIplus);
+      v2 = (cos*HRminus + sin*HIplus);
+      *A = (HRplus  + v1) * (fft_type)0.5;
+      *B = *A - v1;
+      *(A+1) = (HIminus + v2) * (fft_type)0.5;
+      *(B+1) = *(A+1) - HIminus;
+
+      br1++;
+      br2--;
+   }
+   /* Handle the center bin (just need a conjugate) */
+   A=buffer+*br1+1;
+   *A=-*A;
+   /* Handle DC bin separately - and ignore the Fs/2 bin
+   buffer[0]+=buffer[1];
+   buffer[1]=(fft_type)0;*/
+   /* Handle DC and Fs/2 bins separately */
+   /* Put the Fs/2 value into the imaginary part of the DC bin */
+   v1=buffer[0]-buffer[1];
+   buffer[0]+=buffer[1];
+   buffer[1]=v1;
+}
+
+
+/* Description: This routine performs an inverse FFT to real data.
+*              This code is for floating point data.
+*
+*  Note: Output is BIT-REVERSED! so you must use the BitReversed to
+*        get legible output, (i.e. wave[2*i]   = buffer[ BitReversed[i] ]
+*                                  wave[2*i+1] = buffer[ BitReversed[i]+1 ] )
+*        Input is in normal order, interleaved (real,imaginary) complex data
+*        You must call GetFFT(fftlen) first to initialize some buffers!
+*
+* Input buffer[0] is the DC bin, and input buffer[1] is the Fs/2 bin
+* - this can be done because both values will always be real only
+* - this allows us to not have to allocate an extra complex value for the Fs/2 bin
+*
+*  Note: The scaling on this is done according to the standard FFT definition,
+*        so a unit amplitude DC signal will output an amplitude of (N)
+*        (Older revisions would progressively scale the input, so the output
+*        values would be similar in amplitude to the input values, which is
+*        good when using fixed point arithmetic)
+*/
+void InverseRealFFTf(fft_type *buffer, const FFTParam *h)
+{
+   fft_type *A,*B;
+   const fft_type *sptr;
+   const fft_type *endptr1,*endptr2;
+   const int *br1;
+   fft_type HRplus,HRminus,HIplus,HIminus;
+   fft_type v1,v2,sin,cos;
+
+   auto ButterfliesPerGroup = h->Points / 2;
+
+   /* Massage input to get the input for a real output sequence. */
+   A = buffer + 2;
+   B = buffer + h->Points * 2 - 2;
+   br1 = h->BitReversed.get() + 1;
+   while(A<B)
+   {
+      sin=h->SinTable[*br1];
+      cos=h->SinTable[*br1+1];
+      HRplus = (HRminus = *A     - *B    ) + (*B     * 2);
+      HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) * 2);
+      v1 = (sin*HRminus + cos*HIplus);
+      v2 = (cos*HRminus - sin*HIplus);
+      *A = (HRplus  + v1) * (fft_type)0.5;
+      *B = *A - v1;
+      *(A+1) = (HIminus - v2) * (fft_type)0.5;
+      *(B+1) = *(A+1) - HIminus;
+
+      A+=2;
+      B-=2;
+      br1++;
+   }
+   /* Handle center bin (just need conjugate) */
+   *(A+1)=-*(A+1);
+   /* Handle DC bin separately - this ignores any Fs/2 component
+   buffer[1]=buffer[0]=buffer[0]/2;*/
+   /* Handle DC and Fs/2 bins specially */
+   /* The DC bin is passed in as the real part of the DC complex value */
+   /* The Fs/2 bin is passed in as the imaginary part of the DC complex value */
+   /* (v1+v2) = buffer[0] == the DC component */
+   /* (v1-v2) = buffer[1] == the Fs/2 component */
+   v1=0.5f*(buffer[0]+buffer[1]);
+   v2=0.5f*(buffer[0]-buffer[1]);
+   buffer[0]=v1;
+   buffer[1]=v2;
+
+   /*
+   *  Butterfly:
+   *     Ain-----Aout
+   *         \ /
+   *         / \
+   *     Bin-----Bout
+   */
+
+   endptr1 = buffer + h->Points * 2;
+
+   while(ButterfliesPerGroup > 0)
+   {
+      A = buffer;
+      B = buffer + ButterfliesPerGroup * 2;
+      sptr = h->SinTable.get();
+
+      while(A < endptr1)
+      {
+         sin = *(sptr++);
+         cos = *(sptr++);
+         endptr2 = B;
+         while(A < endptr2)
+         {
+            v1 = *B * cos - *(B + 1) * sin;
+            v2 = *B * sin + *(B + 1) * cos;
+            *B = (*A + v1) * (fft_type)0.5;
+            *(A++) = *(B++) - v1;
+            *B = (*A + v2) * (fft_type)0.5;
+            *(A++) = *(B++) - v2;
+         }
+         A = B;
+         B += ButterfliesPerGroup * 2;
+      }
+      ButterfliesPerGroup >>= 1;
+   }
+}
+
+void ReorderToFreq(const FFTParam *hFFT, const fft_type *buffer,
+		   fft_type *RealOut, fft_type *ImagOut)
+{
+   // Copy the data into the real and imaginary outputs
+   for(size_t i = 1; i < hFFT->Points; i++) {
+      RealOut[i] = buffer[hFFT->BitReversed[i]  ];
+      ImagOut[i] = buffer[hFFT->BitReversed[i]+1];
+   }
+   RealOut[0] = buffer[0]; // DC component
+   ImagOut[0] = 0;
+   RealOut[hFFT->Points] = buffer[1]; // Fs/2 component
+   ImagOut[hFFT->Points] = 0;
+}
+
+void ReorderToTime(const FFTParam *hFFT, const fft_type *buffer, fft_type *TimeOut)
+{
+   // Copy the data into the real outputs
+   for(size_t i = 0; i < hFFT->Points; i++) {
+      TimeOut[i*2  ]=buffer[hFFT->BitReversed[i]  ];
+      TimeOut[i*2+1]=buffer[hFFT->BitReversed[i]+1];
+   }
+}

libraries/lib-math/RealFFTf.h

@@ -0,0 +1,32 @@
+#ifndef __realfftf_h
+#define __realfftf_h
+
+#include "MemoryX.h"
+
+using fft_type = float;
+struct FFTParam {
+   ArrayOf<int> BitReversed;
+   ArrayOf<fft_type> SinTable;
+   size_t Points;
+#ifdef EXPERIMENTAL_EQ_SSE_THREADED
+   int pow2Bits;
+#endif
+};
+
+struct MATH_API FFTDeleter{
+   void operator () (FFTParam *p) const;
+};
+
+using HFFT = std::unique_ptr<
+   FFTParam, FFTDeleter
+>;
+
+MATH_API HFFT GetFFT(size_t);
+MATH_API void RealFFTf(fft_type *, const FFTParam *);
+MATH_API void InverseRealFFTf(fft_type *, const FFTParam *);
+MATH_API void ReorderToTime(const FFTParam *hFFT, const fft_type *buffer, fft_type *TimeOut);
+MATH_API void ReorderToFreq(const FFTParam *hFFT, const fft_type *buffer,
+		   fft_type *RealOut, fft_type *ImagOut);
+
+#endif
+

libraries/lib-math/SampleCount.cpp

@@ -0,0 +1,26 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  @file SampleCount.cpp
+
+  Paul Licameli split from audacity/Types.h
+
+**********************************************************************/
+#include "SampleCount.h"
+
+#include <algorithm>
+#include <wx/debug.h>
+
+size_t sampleCount::as_size_t() const {
+   wxASSERT(value >= 0);
+   wxASSERT(static_cast<std::make_unsigned<type>::type>(value) <= std::numeric_limits<size_t>::max());
+   return value;
+}
+
+size_t limitSampleBufferSize( size_t bufferSize, sampleCount limit )
+{
+   return
+      std::min( sampleCount( bufferSize ), std::max( sampleCount(0), limit ) )
+         .as_size_t();
+}

libraries/lib-math/SampleCount.h

@@ -0,0 +1,134 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  @file SampleCount.h
+
+  Paul Licameli split from audacity/Types.h
+
+**********************************************************************/
+#ifndef __AUDACITY_SAMPLE_COUNT__
+#define __AUDACITY_SAMPLE_COUNT__
+
+#include <cstddef>
+
+//! Positions or offsets within audio files need a wide type
+/*! This type disallows implicit interconversions with narrower types */
+class MATH_API sampleCount
+{
+public:
+   using type = long long;
+   static_assert(sizeof(type) == 8, "Wrong width of sampleCount");
+
+   sampleCount () : value { 0 } {}
+
+   // Allow implicit conversion from integral types
+   sampleCount ( type v ) : value { v } {}
+   sampleCount ( unsigned long long v ) : value ( v ) {}
+   sampleCount ( int v ) : value { v } {}
+   sampleCount ( unsigned v ) : value { v } {}
+   sampleCount ( long v ) : value { v } {}
+
+   // unsigned long is 64 bit on some platforms.  Let it narrow.
+   sampleCount ( unsigned long v ) : value ( v ) {}
+
+   // Beware implicit conversions from floating point values!
+   // Otherwise the meaning of binary operators with sampleCount change
+   // their meaning when sampleCount is not an alias!
+   explicit sampleCount ( float f ) : value ( f ) {}
+   explicit sampleCount ( double d ) : value ( d ) {}
+
+   sampleCount ( const sampleCount& ) = default;
+   sampleCount &operator= ( const sampleCount& ) = default;
+
+   float as_float() const { return value; }
+   double as_double() const { return value; }
+
+   long long as_long_long() const { return value; }
+
+   size_t as_size_t() const;
+
+   sampleCount &operator += (sampleCount b) { value += b.value; return *this; }
+   sampleCount &operator -= (sampleCount b) { value -= b.value; return *this; }
+   sampleCount &operator *= (sampleCount b) { value *= b.value; return *this; }
+   sampleCount &operator /= (sampleCount b) { value /= b.value; return *this; }
+   sampleCount &operator %= (sampleCount b) { value %= b.value; return *this; }
+
+   sampleCount operator - () const { return -value; }
+
+   sampleCount &operator ++ () { ++value; return *this; }
+   sampleCount operator ++ (int)
+      { sampleCount result{ *this }; ++value; return result; }
+
+   sampleCount &operator -- () { --value; return *this; }
+   sampleCount operator -- (int)
+      { sampleCount result{ *this }; --value; return result; }
+
+private:
+   type value;
+};
+
+inline bool operator == (sampleCount a, sampleCount b)
+{
+   return a.as_long_long() == b.as_long_long();
+}
+
+inline bool operator != (sampleCount a, sampleCount b)
+{
+   return !(a == b);
+}
+
+inline bool operator < (sampleCount a, sampleCount b)
+{
+   return a.as_long_long() < b.as_long_long();
+}
+
+inline bool operator >= (sampleCount a, sampleCount b)
+{
+   return !(a < b);
+}
+
+inline bool operator > (sampleCount a, sampleCount b)
+{
+   return b < a;
+}
+
+inline bool operator <= (sampleCount a, sampleCount b)
+{
+   return !(b < a);
+}
+
+inline sampleCount operator + (sampleCount a, sampleCount b)
+{
+   return sampleCount{ a } += b;
+}
+
+inline sampleCount operator - (sampleCount a, sampleCount b)
+{
+   return sampleCount{ a } -= b;
+}
+
+inline sampleCount operator * (sampleCount a, sampleCount b)
+{
+   return sampleCount{ a } *= b;
+}
+
+inline sampleCount operator / (sampleCount a, sampleCount b)
+{
+   return sampleCount{ a } /= b;
+}
+
+inline sampleCount operator % (sampleCount a, sampleCount b)
+{
+   return sampleCount{ a } %= b;
+}
+
+// ----------------------------------------------------------------------------
+// Function returning the minimum of a sampleCount and a size_t,
+// hiding the casts
+// ----------------------------------------------------------------------------
+
+MATH_API
+size_t limitSampleBufferSize( size_t bufferSize, sampleCount limit );
+
+#endif

libraries/lib-math/SampleFormat.cpp

@@ -0,0 +1,120 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  SampleFormat.h
+
+  Dominic Mazzoni
+
+*******************************************************************//*!
+
+\file SampleFormat.cpp
+\brief Functions that work with Dither and initialise it.
+
+
+  This file handles converting between all of the different
+  sample formats that Audacity supports, such as 16-bit,
+  24-bit (packed into a 32-bit int), and 32-bit float.
+
+  Floating-point samples use the range -1.0...1.0, inclusive.
+  Integer formats use the full signed range of their data type,
+  for example 16-bit samples use the range -32768...32767.
+  This means that reading in a wav file and writing it out again
+  ('round tripping'), via floats, is lossless; -32768 equates to -1.0f
+  and 32767 equates to +1.0f - (a little bit).
+  It also means (unfortunately) that writing out +1.0f leads to
+  clipping by 1 LSB.  This creates some distortion, but I (MJS) have
+  not been able to measure it, it's so small.  Zero is preserved.
+
+  http://limpet.net/audacity/bugzilla/show_bug.cgi?id=200
+  leads to some of the discussions that were held about this.
+
+   Note: These things are now handled by the Dither class, which
+         also replaces the CopySamples() method (msmeyer)
+
+*//*******************************************************************/
+
+#include "SampleFormat.h"
+#include "Dither.h" // CYCLE
+
+#include <wx/intl.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "Prefs.h"
+#include "Internat.h"
+
+DitherType gLowQualityDither = DitherType::none;
+DitherType gHighQualityDither = DitherType::shaped;
+static Dither gDitherAlgorithm;
+
+void InitDitherers()
+{
+   // Read dither preferences
+   gLowQualityDither = Dither::FastDitherChoice();
+   gHighQualityDither = Dither::BestDitherChoice();
+}
+
+TranslatableString GetSampleFormatStr(sampleFormat format)
+{
+   switch(format) {
+   case int16Sample:
+      /* i18n-hint: Audio data bit depth (precision): 16-bit integers */
+      return XO("16-bit PCM");
+   case int24Sample:
+      /* i18n-hint: Audio data bit depth (precision): 24-bit integers */
+      return XO("24-bit PCM");
+   case floatSample:
+      /* i18n-hint: Audio data bit depth (precision): 32-bit floating point */
+      return XO("32-bit float");
+   }
+   return XO("Unknown format"); // compiler food
+}
+
+// TODO: Risky?  Assumes 0.0f is represented by 0x00000000;
+void ClearSamples(samplePtr dst, sampleFormat format,
+                  size_t start, size_t len)
+{
+   auto size = SAMPLE_SIZE(format);
+   memset(dst + start*size, 0, len*size);
+}
+
+void ReverseSamples(samplePtr dst, sampleFormat format,
+                  int start, int len)
+{
+   auto size = SAMPLE_SIZE(format);
+   samplePtr first = dst + start * size;
+   samplePtr last = dst + (start + len - 1) * size;
+   enum : size_t { fixedSize = SAMPLE_SIZE(floatSample) };
+   wxASSERT(static_cast<size_t>(size) <= fixedSize);
+   char temp[fixedSize];
+   while (first < last) {
+      memcpy(temp, first, size);
+      memcpy(first, last, size);
+      memcpy(last, temp, size);
+      first += size;
+      last -= size;
+   }
+}
+
+void  SamplesToFloats(constSamplePtr src, sampleFormat srcFormat,
+   float *dst, size_t len, size_t srcStride, size_t dstStride)
+{
+   CopySamples( src, srcFormat,
+      reinterpret_cast<samplePtr>(dst), floatSample, len,
+      DitherType::none,
+      srcStride, dstStride);
+}
+
+void CopySamples(constSamplePtr src, sampleFormat srcFormat,
+   samplePtr dst, sampleFormat dstFormat, size_t len,
+   DitherType ditherType, /* = gHighQualityDither */
+   unsigned int srcStride /* = 1 */,
+   unsigned int dstStride /* = 1 */)
+{
+   gDitherAlgorithm.Apply(
+      ditherType,
+      src, srcFormat, dst, dstFormat, len, srcStride, dstStride);
+}

libraries/lib-math/SampleFormat.h

@@ -0,0 +1,190 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  @file SampleFormat.h
+
+  Dominic Mazzoni
+
+**********************************************************************/
+
+#ifndef __AUDACITY_SAMPLE_FORMAT__
+#define __AUDACITY_SAMPLE_FORMAT__
+
+#include "MemoryX.h"
+#include <cstdlib>
+
+//
+// Definitions / Meta-Information
+//
+
+enum DitherType : unsigned;
+//! These global variables are assigned at application startup or after change of preferences.
+extern MATH_API DitherType gLowQualityDither, gHighQualityDither;
+
+// ----------------------------------------------------------------------------
+// Supported sample formats
+// ----------------------------------------------------------------------------
+enum sampleFormat : unsigned
+{
+   //! The increasing sequence of these enum values must correspond to the increasing data type width
+   //! These values persist in saved project files, so must not be changed in later program versions
+   int16Sample = 0x00020001,
+   int24Sample = 0x00040001,
+   floatSample = 0x0004000F,
+
+   //! Two synonyms for previous values that might change if more values were added
+   narrowestSampleFormat = int16Sample,
+   widestSampleFormat = floatSample,
+};
+
+// ----------------------------------------------------------------------------
+// Provide the number of bytes a specific sample will take
+// ----------------------------------------------------------------------------
+#define SAMPLE_SIZE(SampleFormat) (SampleFormat >> 16)
+
+// ----------------------------------------------------------------------------
+// Generic pointer to sample data
+// ----------------------------------------------------------------------------
+using samplePtr = char *;
+using constSamplePtr = const char *;
+
+// Used to determine how to fill in empty areas of audio.
+typedef enum {
+   fillZero = 0,
+   fillTwo = 2
+}fillFormat;
+
+/** \brief Return the size on disk of one uncompressed sample (bytes) */
+#define SAMPLE_SIZE_DISK(SampleFormat) (((SampleFormat) == int24Sample) ? \
+   size_t{ 3 } : SAMPLE_SIZE(SampleFormat) )
+
+class TranslatableString;
+MATH_API TranslatableString GetSampleFormatStr(sampleFormat format);
+
+//
+// Allocating/Freeing Samples
+//
+
+class SampleBuffer {
+
+public:
+   SampleBuffer()
+      : mPtr(0)
+   {}
+   SampleBuffer(size_t count, sampleFormat format)
+      : mPtr((samplePtr)malloc(count * SAMPLE_SIZE(format)))
+   {}
+   ~SampleBuffer()
+   {
+      Free();
+   }
+
+   // WARNING!  May not preserve contents.
+   SampleBuffer &Allocate(size_t count, sampleFormat format)
+   {
+      Free();
+      mPtr = (samplePtr)malloc(count * SAMPLE_SIZE(format));
+      return *this;
+   }
+
+
+   void Free()
+   {
+      free(mPtr);
+      mPtr = 0;
+   }
+
+   samplePtr ptr() const { return mPtr; }
+
+
+private:
+   samplePtr mPtr;
+};
+
+class GrowableSampleBuffer : private SampleBuffer
+{
+public:
+   GrowableSampleBuffer()
+      : SampleBuffer()
+      , mCount(0)
+   {}
+
+   GrowableSampleBuffer(size_t count, sampleFormat format)
+      : SampleBuffer(count, format)
+      , mCount(count)
+   {}
+
+   GrowableSampleBuffer &Resize(size_t count, sampleFormat format)
+   {
+      if (!ptr() || mCount < count) {
+         Allocate(count, format);
+         mCount = count;
+      }
+      return *this;
+   }
+
+   void Free()
+   {
+      SampleBuffer::Free();
+      mCount = 0;
+   }
+
+   using SampleBuffer::ptr;
+
+private:
+   size_t mCount;
+};
+
+//
+// Copying, Converting and Clearing Samples
+//
+
+MATH_API
+//! Copy samples from any format into the widest format, which is 32 bit float, with no dithering
+/*!
+ @param src address of source samples
+ @param srcFormat format of source samples, determines sizeof each one
+ @param dst address of floating-point numbers
+ @param len count of samples to copy
+ @param srcStride how many samples to advance src after copying each one
+ @param dstString how many samples to advance dst after copying each one
+ */
+void SamplesToFloats(constSamplePtr src, sampleFormat srcFormat,
+    float *dst, size_t len, size_t srcStride = 1, size_t dstStride = 1);
+
+MATH_API
+//! Copy samples from any format to any other format; apply dithering only if narrowing the format
+/*!
+ @copydetails SamplesToFloats()
+ @param dstFormat format of destination samples, determines sizeof each one
+ @param ditherType choice of dithering algorithm to use if narrowing the format
+ */
+void CopySamples(constSamplePtr src, sampleFormat srcFormat,
+   samplePtr dst, sampleFormat dstFormat, size_t len,
+   DitherType ditherType = gHighQualityDither, //!< default is loaded from a global variable
+   unsigned int srcStride=1, unsigned int dstStride=1);
+
+MATH_API
+void      ClearSamples(samplePtr buffer, sampleFormat format,
+                       size_t start, size_t len);
+
+MATH_API
+void      ReverseSamples(samplePtr buffer, sampleFormat format,
+                         int start, int len);
+
+//
+// This must be called on startup and everytime NEW ditherers
+// are set in preferences.
+//
+
+MATH_API
+void      InitDitherers();
+
+// These are so commonly done for processing samples in floating point form in memory,
+// let's have abbreviations.
+using Floats = ArrayOf<float>;
+using FloatBuffers = ArraysOf<float>;
+using Doubles = ArrayOf<double>;
+
+#endif

libraries/lib-math/Spectrum.cpp

@@ -0,0 +1,125 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  Spectrum.cpp
+
+  Dominic Mazzoni
+
+*******************************************************************//*!
+
+\file Spectrum.cpp
+\brief Functions for computing Spectra.
+
+*//*******************************************************************/
+
+#include "Spectrum.h"
+
+#include <math.h>
+
+#include "SampleFormat.h"
+
+bool ComputeSpectrum(const float * data, size_t width,
+                     size_t windowSize,
+                     double WXUNUSED(rate), float *output,
+                     bool autocorrelation, int windowFunc)
+{
+   if (width < windowSize)
+      return false;
+
+   if (!data || !output)
+      return true;
+
+   Floats processed{ windowSize };
+
+   for (size_t i = 0; i < windowSize; i++)
+      processed[i] = float(0.0);
+   auto half = windowSize / 2;
+
+   Floats in{ windowSize };
+   Floats out{ windowSize };
+   Floats out2{ windowSize };
+
+   size_t start = 0;
+   unsigned windows = 0;
+   while (start + windowSize <= width) {
+      for (size_t i = 0; i < windowSize; i++)
+         in[i] = data[start + i];
+
+      WindowFunc(windowFunc, windowSize, in.get());
+
+      if (autocorrelation) {
+         // Take FFT
+         RealFFT(windowSize, in.get(), out.get(), out2.get());
+         // Compute power
+         for (size_t i = 0; i < windowSize; i++)
+            in[i] = (out[i] * out[i]) + (out2[i] * out2[i]);
+
+         // Tolonen and Karjalainen recommend taking the cube root
+         // of the power, instead of the square root
+
+         for (size_t i = 0; i < windowSize; i++)
+            in[i] = powf(in[i], 1.0f / 3.0f);
+
+         // Take FFT
+         RealFFT(windowSize, in.get(), out.get(), out2.get());
+      }
+      else
+         PowerSpectrum(windowSize, in.get(), out.get());
+
+      // Take real part of result
+      for (size_t i = 0; i < half; i++)
+        processed[i] += out[i];
+
+      start += half;
+      windows++;
+   }
+
+   if (autocorrelation) {
+
+      // Peak Pruning as described by Tolonen and Karjalainen, 2000
+      /*
+       Combine most of the calculations in a single for loop.
+       It should be safe, as indexes refer only to current and previous elements,
+       that have already been clipped, etc...
+      */
+      for (size_t i = 0; i < half; i++) {
+        // Clip at zero, copy to temp array
+        if (processed[i] < 0.0)
+            processed[i] = float(0.0);
+        out[i] = processed[i];
+        // Subtract a time-doubled signal (linearly interp.) from the original
+        // (clipped) signal
+        if ((i % 2) == 0)
+           processed[i] -= out[i / 2];
+        else
+           processed[i] -= ((out[i / 2] + out[i / 2 + 1]) / 2);
+
+        // Clip at zero again
+        if (processed[i] < 0.0)
+            processed[i] = float(0.0);
+      }
+
+      // Reverse and scale
+      for (size_t i = 0; i < half; i++)
+         in[i] = processed[i] / (windowSize / 4);
+      for (size_t i = 0; i < half; i++)
+         processed[half - 1 - i] = in[i];
+   } else {
+      // Convert to decibels
+      // But do it safely; -Inf is nobody's friend
+      for (size_t i = 0; i < half; i++){
+         float temp=(processed[i] / windowSize / windows);
+         if (temp > 0.0)
+            processed[i] = 10 * log10(temp);
+         else
+            processed[i] = 0;
+      }
+   }
+
+   for(size_t i = 0; i < half; i++)
+      output[i] = processed[i];
+
+   return true;
+}
+

libraries/lib-math/Spectrum.h

@@ -0,0 +1,29 @@
+/**********************************************************************
+
+  Audacity: A Digital Audio Editor
+
+  Spectrum.h
+
+  Dominic Mazzoni
+
+**********************************************************************/
+
+#ifndef __AUDACITY_SPECTRUM__
+#define __AUDACITY_SPECTRUM__
+
+#include "FFT.h"
+
+/*
+  This function computes the power (mean square amplitude) as
+  a function of frequency, for some block of audio data.
+
+  width: the number of samples
+  calculates windowSize/2 frequency samples
+*/
+
+MATH_API
+bool ComputeSpectrum(const float * data, size_t width, size_t windowSize,
+                     double rate, float *out, bool autocorrelation,
+                     int windowFunc = eWinFuncHann);
+
+#endif

libraries/lib-math/float_cast.h

@@ -0,0 +1,172 @@
+/*
+** Copyright (C) 2001 Erik de Castro Lopo <erikd AT mega-nerd DOT com>
+**
+** Permission to use, copy, modify, distribute, and sell this file for any
+** purpose is hereby granted without fee, provided that the above copyright
+** and this permission notice appear in all copies.  No representations are
+** made about the suitability of this software for any purpose.  It is
+** provided "as is" without express or implied warranty.
+*/
+
+/* Version 1.1 */
+
+
+/*============================================================================
+**	On Intel Pentium processors (especially PIII and probably P4), converting
+**	from float to int is very slow. To meet the C specs, the code produced by
+**	most C compilers targeting Pentium needs to change the FPU rounding mode
+**	before the float to int conversion is performed.
+**
+**	Changing the FPU rounding mode causes the FPU pipeline to be flushed. It
+**	is this flushing of the pipeline which is so slow.
+**
+**	Fortunately the ISO C99 specifications define the functions lrint, lrintf,
+**	llrint and llrintf which fix this problem as a side effect.
+**
+**	On Unix-like systems, the configure process should have detected the
+**	presence of these functions. If they weren't found we have to replace them
+**	here with a standard C cast.
+*/
+
+/*
+**	The C99 prototypes for lrint and lrintf are as follows:
+**
+**		long int lrintf (float x) ;
+**		long int lrint  (double x) ;
+*/
+
+/*	The presence of the required functions are detected during the configure
+**	process and the values HAVE_LRINT and HAVE_LRINTF are set accordingly in
+**	the config.h file.
+*/
+#if (defined (WIN32) || defined (_WIN32)) && defined(_MSC_VER) && defined(_M_IX86)
+	// As of Visual Studio 2019 16.9, these functions have been made intrinsic and the build
+	// will fail. Unfortunately, the intrinsic versions run a LOT slower than the ones
+   // below, so force the compiler to use ours instead.
+	#pragma function( lrint, lrintf )
+
+   // Including math.h allows us to use the inline assembler versions without
+   // producing errors in newer Visual Studio versions.
+   // Without the include, we get different linkage error messages.
+   // Without the inline assembler versions, these functions are VERY slow.
+   // I also see that the include was part of the original source for this file:
+   //    http://www.mega-nerd.com/FPcast/
+   
+   #include <math.h>
+
+   /*	Win32 doesn't seem to have these functions.
+	**	Therefore implement inline versions of these functions here.
+	*/
+	__inline long int
+	lrint (double flt)
+	{	int intgr;
+
+		_asm
+		{	fld flt
+			fistp intgr
+			} ;
+
+		return intgr ;
+	}
+
+	__inline long int
+	lrintf (float flt)
+	{	int intgr;
+
+		_asm
+		{	fld flt
+			fistp intgr
+			} ;
+
+		return intgr ;
+	}
+
+	__inline long long int
+	llrint (double flt)
+	{	long long int intgr;
+
+		_asm
+		{	fld flt
+			fistp intgr
+			} ;
+
+		return intgr ;
+	}
+
+	__inline long long int
+	llrintf (float flt)
+	{	long long int intgr;
+
+		_asm
+		{	fld flt
+			fistp intgr
+			} ;
+
+		return intgr ;
+	}
+#elif (defined (WIN32) || defined (_WIN32)) && defined(_M_X64)
+
+	#include <math.h>
+	#include <immintrin.h>
+	#include <emmintrin.h>
+
+	#ifdef _MSC_VER
+		#pragma function(lrint, lrintf)
+	#endif
+
+	__inline 
+	long int lrint(double flt)
+	{
+       return _mm_cvtsd_si32(_mm_set_sd(flt));
+    }
+
+	__inline 
+	long int lrintf (float flt)
+	{
+		return _mm_cvtss_si32(_mm_set_ss(flt));
+	}
+
+	__inline 
+	long long int llrint(double flt)
+    {
+        return _mm_cvtsd_si64(_mm_set_sd(flt));
+    }
+
+    __inline 
+    long long int llrintf(float flt)
+    {
+       return _mm_cvtss_si64(_mm_set_ss(flt));
+    }
+
+#elif (HAVE_LRINT && HAVE_LRINTF)
+
+	/*	These defines enable functionality introduced with the 1999 ISO C
+	**	standard. They must be defined before the inclusion of math.h to
+	**	engage them. If optimisation is enabled, these functions will be
+	**	inlined. With optimisation switched off, you have to link in the
+	**	maths library using -lm.
+	*/
+
+	#define	_ISOC9X_SOURCE 1
+	#define  _ISOC99_SOURCE 1
+
+	#define	__USE_ISOC9X   1
+	#define	__USE_ISOC99   1
+
+	#include	<math.h>
+#else
+
+   /* dmazzoni: modified these to do a proper rounding, even though
+    * it's slower.  Correctness and consistency is more important
+    * than speed, especially since lrint/lrintf are certainly not
+    * available everywhere.
+    *
+    * MM: Now uses internal math.h rint() function
+    */
+
+   #include	<math.h>
+
+   #define	lrint(dbl)		((int)rint(dbl))
+   #define	lrintf(flt)		((int)rint(flt))
+
+#endif