Add the wavetable generator

src/CMakeLists.txt

@@ -17,6 +17,7 @@ set (SFIZZ_SOURCES
     sfizz/Logger.cpp
     sfizz/SfzFilter.cpp
     sfizz/Curve.cpp
+    sfizz/Wavetables.cpp
     sfizz/Effects.cpp
     sfizz/effects/Nothing.cpp
     sfizz/effects/Filter.cpp

src/sfizz/Wavetables.cpp

@@ -0,0 +1,115 @@
+// SPDX-License-Identifier: BSD-2-Clause
+
+// This code is part of the sfizz library and is licensed under a BSD 2-clause
+// license. You should have receive a LICENSE.md file along with the code.
+// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz
+
+#include "Wavetables.h"
+#include <kiss_fftr.h>
+#include <memory>
+
+namespace sfz {
+
+void HarmonicProfile::generate(
+    absl::Span<float> table, double amplitude, double cutoff) const
+{
+    size_t size = table.size();
+
+    typedef std::complex<kiss_fft_scalar> cpx;
+
+    // allocate a spectrum of size N/2+1
+    // bins are equispaced in frequency, with index N/2 being nyquist
+    std::unique_ptr<cpx[]> spec(new cpx[size / 2 + 1]());
+
+    kiss_fftr_cfg cfg = kiss_fftr_alloc(size, true, nullptr, nullptr);
+    if (!cfg)
+        throw std::bad_alloc();
+
+    // bins need scaling and phase offset; this IFFT is a sum of cosines
+    const std::complex<double> k = std::polar(amplitude * 0.5, M_PI / 2);
+
+    // start filling at bin index 1; 1 is fundamental, 0 is DC
+    for (size_t index = 1; index < size / 2 + 1; ++index) {
+        if (index * (1.0 / size) > cutoff)
+            break;
+
+        std::complex<double> harmonic = getHarmonic(index);
+        spec[index] = k * harmonic;
+    }
+
+    kiss_fftri(cfg, reinterpret_cast<kiss_fft_cpx*>(spec.get()), table.data());
+    kiss_fftr_free(cfg);
+}
+
+class SineProfile : public HarmonicProfile {
+public:
+    std::complex<double> getHarmonic(size_t index) const
+    {
+        return (index == 1) ? 1.0 : 0.0;
+    }
+};
+
+class TriangleProfile : public HarmonicProfile {
+public:
+    std::complex<double> getHarmonic(size_t index) const
+    {
+        if ((index & 1) == 0)
+            return 0.0;
+
+        bool s = (index >> 1) & 1;
+        return std::polar<double>(
+            (8 / (M_PI * M_PI)) * (1.0 / (index * index)),
+            s ? 0.0 : M_PI);
+    }
+};
+
+class SawProfile : public HarmonicProfile {
+public:
+    std::complex<double> getHarmonic(size_t index) const
+    {
+        if (index < 1)
+            return 0.0;
+
+        return std::polar(
+            (2.0 / M_PI) / index,
+            (index & 1) ? 0.0 : M_PI);
+    }
+};
+
+class SquareProfile : public HarmonicProfile {
+public:
+    std::complex<double> getHarmonic(size_t index) const
+    {
+        if ((index & 1) == 0)
+            return 0.0;
+
+        return std::polar((4.0 / M_PI) / index, M_PI);
+    }
+};
+
+static const SineProfile sineProfile;
+static const TriangleProfile triangleProfile;
+static const SawProfile sawProfile;
+static const SquareProfile squareProfile;
+
+const HarmonicProfile& HarmonicProfile::getSine()
+{
+    return sineProfile;
+}
+
+const HarmonicProfile& HarmonicProfile::getTriangle()
+{
+    return triangleProfile;
+}
+
+const HarmonicProfile& HarmonicProfile::getSaw()
+{
+    return sawProfile;
+}
+
+const HarmonicProfile& HarmonicProfile::getSquare()
+{
+    return squareProfile;
+}
+
+} // namespace sfz

src/sfizz/Wavetables.h

@@ -0,0 +1,41 @@
+// SPDX-License-Identifier: BSD-2-Clause
+
+// This code is part of the sfizz library and is licensed under a BSD 2-clause
+// license. You should have receive a LICENSE.md file along with the code.
+// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz
+
+#pragma once
+#include <absl/types/span.h>
+#include <complex>
+
+namespace sfz {
+
+/**
+   A description of the harmonics of a particular wave form
+ */
+class HarmonicProfile {
+public:
+    virtual ~HarmonicProfile() {}
+
+    static const HarmonicProfile& getSine();
+    static const HarmonicProfile& getTriangle();
+    static const HarmonicProfile& getSaw();
+    static const HarmonicProfile& getSquare();
+
+    /**
+       @brief Get the value at the given index of the frequency spectrum.
+
+       The modulus and the argument of the complex number are equal to the
+       amplitude and the phase of the harmonic component.
+     */
+    virtual std::complex<double> getHarmonic(size_t index) const = 0;
+
+    /**
+       @brief Generate a period of the waveform and store it in the table.
+
+       Do not generate harmonics above cutoff, which is expressed as Fc/Fs.
+     */
+    void generate(absl::Span<float> table, double amplitude, double cutoff) const;
+};
+
+} // namespace sfz

tests/CMakeLists.txt

@@ -60,4 +60,7 @@ target_link_libraries(filter_apply PRIVATE sfizz::sfizz)
 add_executable(sfizz_plot_curve PlotCurve.cpp)
 target_link_libraries(sfizz_plot_curve PRIVATE sfizz::sfizz)
 
+add_executable(sfizz_plot_wavetables PlotWavetables.cpp)
+target_link_libraries(sfizz_plot_wavetables PRIVATE sfizz::sfizz)
+
 file(COPY "." DESTINATION ${CMAKE_BINARY_DIR}/tests)

tests/PlotWavetables.cpp

@@ -0,0 +1,63 @@
+// SPDX-License-Identifier: BSD-2-Clause
+
+// This code is part of the sfizz library and is licensed under a BSD 2-clause
+// license. You should have receive a LICENSE.md file along with the code.
+// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz
+
+#include "sfizz/Wavetables.h"
+#include <absl/strings/numbers.h>
+#include <absl/strings/string_view.h>
+#include <iostream>
+
+static void usage()
+{
+    std::cerr << "Usage: sfizz_plot_wavetables [-w wave] [-a amplitude] [-c cutoff]\n";
+}
+
+int main(int argc, char* argv[])
+{
+    absl::string_view waveName;
+    double amplitude = 1.0;
+    double cutoff = 0.5;
+
+    for (int i = 1; i < argc; ++i) {
+        absl::string_view arg = argv[i];
+
+        if (arg == "-w") {
+            if (i + 1 >= argc)
+                return usage(), 1;
+            waveName = argv[++i];
+        } else if (arg == "-a") {
+            if (i + 1 >= argc || !absl::SimpleAtod(argv[++i], &amplitude))
+                return usage(), 1;
+        } else if (arg == "-c") {
+            if (i + 1 >= argc || !absl::SimpleAtod(argv[++i], &cutoff))
+                return usage(), 1;
+        } else
+            return usage(), 1;
+    }
+
+    const sfz::HarmonicProfile* hp = nullptr;
+    if (waveName == "sine")
+        hp = &sfz::HarmonicProfile::getSine();
+    else if (waveName == "square")
+        hp = &sfz::HarmonicProfile::getSquare();
+    else if (waveName == "triangle")
+        hp = &sfz::HarmonicProfile::getTriangle();
+    else if (waveName == "saw")
+        hp = &sfz::HarmonicProfile::getSaw();
+    else {
+        std::cerr << "Unknown wave: " << waveName << '\n';
+        return 1;
+    }
+
+    constexpr size_t tableSize = 2048;
+    float table[tableSize];
+
+    hp->generate(table, amplitude, cutoff);
+
+    for (size_t i = 0; i < tableSize; ++i)
+        std::cout << (i * (1.0 / tableSize)) << ' ' << table[i] << '\n';
+
+    return 0;
+}