diff --git a/src/sfizz/Wavetables.cpp b/src/sfizz/Wavetables.cpp
index 4562240c..058466d7 100644
--- a/src/sfizz/Wavetables.cpp
+++ b/src/sfizz/Wavetables.cpp
@@ -142,4 +142,36 @@ WavetableRange WavetableRange::getRangeForFrequency(float f)
     return getRangeForOctave(oct);
 }
 
+//------------------------------------------------------------------------------
+WavetableMulti WavetableMulti::createForHarmonicProfile(
+    const HarmonicProfile& hp, unsigned tableSize, double refSampleRate)
+{
+    WavetableMulti wm;
+    constexpr unsigned multiSize = wm.multiSize();
+
+    // amplitude to match ARIA's default generator output
+    constexpr double amplitude = 0.25;
+
+    wm._tableSize = tableSize;
+    wm._multiData.reset(new float[tableSize * multiSize]);
+
+    for (unsigned m = 0; m < multiSize; ++m) {
+        WavetableRange range = WavetableRange::getRangeForOctave(m);
+
+        double freq = range.maxFrequency;
+
+        // A spectrum S of fundamental F has: S[1]=F and S[N/2]=Fs'/2
+        // which lets it generate frequency up to Fs'/2=F*N/2.
+        // Therefore it's desired to cut harmonics at C=0.5*Fs/Fs'=0.5*Fs/(F*N).
+        double cutoff = (0.5 * refSampleRate / tableSize) / freq;
+
+        float* ptr = &wm._multiData[m * tableSize];
+        absl::Span<float> table(ptr, tableSize);
+
+        hp.generate(table, amplitude, cutoff);
+    }
+
+    return wm;
+}
+
 } // namespace sfz
diff --git a/src/sfizz/Wavetables.h b/src/sfizz/Wavetables.h
index 9226ce44..317d88ab 100644
--- a/src/sfizz/Wavetables.h
+++ b/src/sfizz/Wavetables.h
@@ -6,6 +6,7 @@
 
 #pragma once
 #include <absl/types/span.h>
+#include <memory>
 #include <complex>
 
 namespace sfz {
@@ -70,4 +71,44 @@ public:
     //     octave 9: 10240 Hz - 20480 Hz
 };
 
+/**
+   Multisample of a wavetable, which is a collection of FFT-filtered mipmaps
+   adapted for various playback frequencies.
+ */
+class WavetableMulti {
+public:
+    // number of tables in the multisample
+    unsigned tableSize() const { return _tableSize; }
+
+    // number of tables in the multisample
+    static constexpr unsigned multiSize() { return WavetableRange::countOctaves; }
+
+    // get the N-th table in the multisample
+    absl::Span<const float> getTable(unsigned index) const
+    {
+        unsigned size = _tableSize;
+        const float* ptr = &_multiData[index * _tableSize];
+        return { ptr, size };
+    }
+
+    // get the table which is adequate for a given playback frequency
+    absl::Span<const float> getTableForFrequency(float freq) const
+    {
+        return getTable(WavetableRange::getOctaveForFrequency(freq));
+    }
+
+    // create a multisample according to a given harmonic profile
+    // the reference sample rate is the minimum value accepted by the DSP
+    // system (most defavorable wrt. aliasing)
+    static WavetableMulti createForHarmonicProfile(
+        const HarmonicProfile& hp, unsigned tableSize, double refSampleRate = 44100.0);
+
+private:
+    // length of each individual table of the multisample
+    unsigned _tableSize = 0;
+
+    // internal storage, having `multiSize` rows and `tableSize` columns.
+    std::unique_ptr<float[]> _multiData;
+};
+
 } // namespace sfz
diff --git a/tests/PlotWavetables.cpp b/tests/PlotWavetables.cpp
index 6037c549..f37a6897 100644
--- a/tests/PlotWavetables.cpp
+++ b/tests/PlotWavetables.cpp
@@ -11,7 +11,7 @@
 
 static void usage()
 {
-    std::cerr << "Usage: sfizz_plot_wavetables [-w wave] [-a amplitude] [-c cutoff]\n";
+    std::cerr << "Usage: sfizz_plot_wavetables [-w wave] [-a amplitude] [-c cutoff] [-m]\n";
 }
 
 int main(int argc, char* argv[])
@@ -19,6 +19,7 @@ int main(int argc, char* argv[])
     absl::string_view waveName;
     double amplitude = 1.0;
     double cutoff = 0.5;
+    bool generateMulti = false;
 
     for (int i = 1; i < argc; ++i) {
         absl::string_view arg = argv[i];
@@ -33,10 +34,15 @@ int main(int argc, char* argv[])
         } else if (arg == "-c") {
             if (i + 1 >= argc || !absl::SimpleAtod(argv[++i], &cutoff))
                 return usage(), 1;
-        } else
+        } else if (arg == "-m")
+            generateMulti = true;
+        else
             return usage(), 1;
     }
 
+    if (waveName.empty())
+        waveName = "saw";
+
     const sfz::HarmonicProfile* hp = nullptr;
     if (waveName == "sine")
         hp = &sfz::HarmonicProfile::getSine();
@@ -54,10 +60,37 @@ int main(int argc, char* argv[])
     constexpr size_t tableSize = 2048;
     float table[tableSize];
 
-    hp->generate(table, amplitude, cutoff);
+    if (!generateMulti) {
+        hp->generate(table, amplitude, cutoff);
 
-    for (size_t i = 0; i < tableSize; ++i)
-        std::cout << (i * (1.0 / tableSize)) << ' ' << table[i] << '\n';
+        for (size_t i = 0; i < tableSize; ++i)
+            std::cout << (i * (1.0 / tableSize)) << ' ' << table[i] << '\n';
+    } else {
+        sfz::WavetableMulti multi = sfz::WavetableMulti::createForHarmonicProfile(
+            *hp, tableSize);
+        unsigned multiSize = multi.multiSize();
+
+        if (true) {
+            // print all tables one after another
+            for (unsigned m = 0; m < multiSize; ++m) {
+                absl::Span<const float> table = multi.getTable(m);
+                for (size_t i = 0; i < tableSize; ++i) {
+                    std::cout << ((i + m * tableSize) * (1.0 / tableSize))
+                              << ' ' << table[i] << '\n';
+                }
+            }
+        } else {
+            // print all tables separately
+            for (size_t i = 0; i < tableSize; ++i) {
+                std::cout << (i * (1.0 / tableSize));
+                for (unsigned m = 0; m < multiSize; ++m) {
+                    absl::Span<const float> table = multi.getTable(m);
+                    std::cout << ' ' << table[i];
+                }
+                std::cout << '\n';
+            }
+        }
+    }
 
     return 0;
 }