prelude to epicness

.gitignore

@@ -5,3 +5,10 @@ build_macos/
 build_windows/
 icons
 *.png
+vamp-plugin/
+tests/
+ADC-2024/
+*.json
+LICENCE
+readme.md
+vamp-plugin-sdk.cmake

CMakeLists.txt

@@ -67,6 +67,12 @@ else()
     FetchContent_MakeAvailable(fftw3_source)
 endif()
 
+# --- Loop Tempo Estimator ---
+# Disable tests and vamp plugin for the library to speed up build and reduce deps
+set(BUILD_TESTS OFF CACHE BOOL "Build tests" FORCE)
+set(BUILD_VAMP_PLUGIN OFF CACHE BOOL "Build Vamp plugin" FORCE)
+add_subdirectory(libraries/loop-tempo-estimator)
+
 # --- Icon Generation ---
 
 set(ICON_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/assets/icon_source.png")
@@ -151,6 +157,7 @@ endif()
 target_link_libraries(YrCrystals PRIVATE
     Qt6::Core Qt6::Gui Qt6::Widgets Qt6::Multimedia Qt6::OpenGLWidgets
     ${FFTW_TARGET}
+    loop-tempo-estimator
 )
 
 if(BUILD_ANDROID)

libraries/loop-tempo-estimator

@@ -0,0 +1 @@
+Subproject commit 9369d20106e8cbd719b8cb11b3db8b97ac403229

src/AudioEngine.cpp

@@ -11,6 +11,38 @@
 #include <QtGlobal>
 #include <algorithm>
 
+// Include Loop Tempo Estimator
+#include "LoopTempoEstimator/LoopTempoEstimator.h"
+
+// Wrapper for LTE::LteAudioReader to read from our memory buffer
+class MemoryAudioReader : public LTE::LteAudioReader {
+public:
+    MemoryAudioReader(const float* data, long long numFrames, int sampleRate)
+        : m_data(data), m_numFrames(numFrames), m_sampleRate(sampleRate) {}
+
+    double GetSampleRate() const override { return static_cast<double>(m_sampleRate); }
+    long long GetNumSamples() const override { return m_numFrames; }
+
+    void ReadFloats(float* buffer, long long where, size_t numFrames) const override {
+        for (size_t i = 0; i < numFrames; ++i) {
+            long long srcIdx = (where + i) * 2; // Stereo interleaved
+            if (srcIdx + 1 < m_numFrames * 2) {
+                // Mix down to mono for analysis
+                float l = m_data[srcIdx];
+                float r = m_data[srcIdx + 1];
+                buffer[i] = (l + r) * 0.5f;
+            } else {
+                buffer[i] = 0.0f;
+            }
+        }
+    }
+
+private:
+    const float* m_data;
+    long long m_numFrames;
+    int m_sampleRate;
+};
+
 AudioEngine::AudioEngine(QObject* parent) : QObject(parent) {
     // Main Processors (Steady State)
     m_processors.push_back(new Processor(m_frameSize, m_sampleRate));
@@ -18,8 +50,8 @@ AudioEngine::AudioEngine(QObject* parent) : QObject(parent) {
     
     // Configure Main: Expander + HPF + Moderate Smoothing
     for(auto p : m_processors) {
-        p->setExpander(1.5f, -50.0f); // Gentle expansion to clean up noise
+        p->setExpander(1.5f, -50.0f); 
-        p->setHPF(60.0f);             // Cut rumble below 60Hz
+        p->setHPF(80.0f);             // Mid 2nd Order HPF (80Hz)
         p->setSmoothing(3);
     }
 
@@ -30,8 +62,9 @@ AudioEngine::AudioEngine(QObject* parent) : QObject(parent) {
 
     // Configure Transient: Aggressive expansion, light smoothing
     for(auto p : m_transientProcessors) {
-        p->setExpander(2.5f, -40.0f); // Expand dynamics around -40dB
+        p->setExpander(2.5f, -40.0f); 
-        p->setSmoothing(2);           // Fast decay
+        p->setHPF(100.0f);            // Clean up transients
+        p->setSmoothing(2);           
     }
 
     m_processTimer = new QTimer(this);
@@ -165,6 +198,29 @@ void AudioEngine::onFinished() {
         emit trackLoaded(false);
         return;
     }
+
+    // --- Run Tempo Estimation ---
+    const float* rawFloats = reinterpret_cast<const float*>(m_pcmData.constData());
+    long long totalFloats = m_pcmData.size() / sizeof(float);
+    long long totalFrames = totalFloats / 2; // Stereo
+
+    if (totalFrames > 0) {
+        MemoryAudioReader reader(rawFloats, totalFrames, m_sampleRate);
+        
+        // Use Lenient tolerance to get a result more often
+        auto bpmOpt = LTE::GetBpm(reader, LTE::FalsePositiveTolerance::Lenient, nullptr);
+        
+        if (bpmOpt.has_value()) {
+            float bpm = static_cast<float>(*bpmOpt);
+            qDebug() << "AudioEngine: Detected BPM:" << bpm;
+            emit analysisReady(bpm, 1.0f);
+        } else {
+            qDebug() << "AudioEngine: No BPM detected.";
+            emit analysisReady(0.0f, 0.0f);
+        }
+    }
+    // ----------------------------
+
     m_buffer.setData(m_pcmData);
     if (!m_buffer.open(QIODevice::ReadOnly)) {
         emit trackLoaded(false);
@@ -281,18 +337,33 @@ void AudioEngine::onProcessTimer() {
     m_transientProcessors[1]->pushData(tCh1);
 
     std::vector<FrameData> results;
+    
+    // Final Compressor Settings
+    float compThreshold = -15.0f;
+    float compRatio = 4.0f;
+
     for (size_t i = 0; i < m_processors.size(); ++i) {
         auto specMain = m_processors[i]->getSpectrum();
         auto specTrans = m_transientProcessors[i]->getSpectrum();
 
+        // Capture Primary DB (Steady State) for Crystal Pattern
+        std::vector<float> primaryDb = specMain.db;
+
         // Mix: Overlay the expanded transient peaks onto the main spectrum
         if (specMain.db.size() == specTrans.db.size()) {
             for(size_t b = 0; b < specMain.db.size(); ++b) {
-                specMain.db[b] = std::max(specMain.db[b], specTrans.db[b]);
+                float val = std::max(specMain.db[b], specTrans.db[b]);
+                
+                // Final Compressor (Hard Knee)
+                if (val > compThreshold) {
+                    val = compThreshold + (val - compThreshold) / compRatio;
+                }
+                
+                specMain.db[b] = val;
             }
         }
 
-        results.push_back({specMain.freqs, specMain.db});
+        results.push_back({specMain.freqs, specMain.db, primaryDb});
     }
     emit spectrumReady(results);
 }

src/AudioEngine.h

@@ -18,7 +18,8 @@ public:
 
     struct FrameData {
         std::vector<float> freqs;
-        std::vector<float> db;
+        std::vector<float> db;        // Mixed (Primary + Transient) -> For Bar Height
+        std::vector<float> primaryDb; // Primary Only -> For Crystal Pattern
     };
 
 public slots:
@@ -35,6 +36,7 @@ signals:
     void positionChanged(float pos);
     void trackLoaded(bool success);
     void spectrumReady(const std::vector<AudioEngine::FrameData>& data);
+    void analysisReady(float bpm, float confidence);
 
 private slots:
     void onBufferReady();

src/MainWindow.cpp

@@ -46,6 +46,7 @@ MainWindow::MainWindow(QWidget* parent) : QMainWindow(parent) {
     connect(m_engine, &AudioEngine::positionChanged, m_playerPage->playback(), &PlaybackWidget::updateSeek);
 
     connect(m_engine, &AudioEngine::spectrumReady, m_playerPage->visualizer(), &VisualizerWidget::updateData);
+    connect(m_engine, &AudioEngine::analysisReady, this, &MainWindow::onAnalysisReady);
 
     audioThread->start();
 }
@@ -297,6 +298,36 @@ void MainWindow::onTrackLoaded(bool success) {
     }
 }
 
+void MainWindow::onAnalysisReady(float bpm, float confidence) {
+    // Feedback Mechanism:
+    // Adjust Entropy based on BPM.
+    // High BPM (Fast/Punchy) -> Lower Entropy Slider (0.5 - 0.8) to allow transients.
+    // Low BPM (Slow/Ambient) -> Higher Entropy Slider (1.0 - 1.5) to smooth noise.
+    // Default (No BPM) -> 1.0
+
+    float targetEntropy = 1.0f;
+
+    if (bpm > 0.0f) {
+        // Map 60..180 BPM to 1.5..0.5 Entropy
+        // Formula: 1.5 - ((bpm - 60) / 120)
+        // Clamped between 0.5 and 1.5
+        float normalized = (bpm - 60.0f) / 120.0f;
+        targetEntropy = 1.5f - normalized;
+        targetEntropy = std::clamp(targetEntropy, 0.5f, 1.5f);
+        qDebug() << "Feedback: BPM" << bpm << "-> Setting Entropy to" << targetEntropy;
+    } else {
+        qDebug() << "Feedback: No BPM -> Default Entropy 1.0";
+    }
+
+    // Update Settings Widget (which updates Visualizer)
+    SettingsWidget* s = m_playerPage->settings();
+    s->setParams(
+        s->isGlass(), s->isFocus(), s->isTrails(), s->isAlbumColors(), 
+        s->isShadow(), s->isMirrored(), s->getBins(), s->getBrightness(), 
+        targetEntropy
+    );
+}
+
 void MainWindow::onTrackDoubleClicked(QListWidgetItem* item) {
     loadIndex(m_playlist->row(item));
 }

src/MainWindow.h

@@ -1,3 +1,5 @@
+// src/MainWindow.h
+
 #pragma once
 #include <QMainWindow>
 #include <QDockWidget>
@@ -25,6 +27,7 @@ private slots:
     void onTrackFinished();
     void onTrackLoaded(bool success);
     void onTrackDoubleClicked(QListWidgetItem* item);
+    void onAnalysisReady(float bpm, float confidence);
     void play();
     void pause();
     void nextTrack();

src/Processor.cpp

@@ -136,9 +136,14 @@ Processor::Spectrum Processor::getSpectrum() {
         
         float freq = i * (float)m_sampleRate / m_frameSize;
         
-        // HPF: Attenuate frequencies below cutoff
+        // HPF: 2nd Order Butterworth Curve
-        if (freq < m_hpfCutoff) {
+        // Gain = 1 / sqrt(1 + (fc/f)^4)
-            mag *= 0.0001f; // -80dB attenuation
+        if (m_hpfCutoff > 0.0f && freq > 0.0f) {
+            float ratio = m_hpfCutoff / freq;
+            float gain = 1.0f / std::sqrt(1.0f + (ratio * ratio * ratio * ratio));
+            mag *= gain;
+        } else if (freq == 0.0f) {
+            mag = 0.0f;
         }
 
         dbFull[i] = 20.0f * std::log10(std::max(mag, 1e-12f));

src/VisualizerWidget.cpp

@@ -91,6 +91,8 @@ void VisualizerWidget::updateData(const std::vector<AudioEngine::FrameData>& dat
 
     for (size_t ch = 0; ch < data.size(); ++ch) {
         const auto& db = data[ch].db;
+        const auto& primaryDb = data[ch].primaryDb; // Access Primary DB
+
         size_t numBins = db.size();
         auto& bins = m_channels[ch].bins;
         if (bins.size() != numBins) bins.resize(numBins);
@@ -98,6 +100,7 @@ void VisualizerWidget::updateData(const std::vector<AudioEngine::FrameData>& dat
         for (size_t i = 0; i < numBins; ++i) {
             auto& bin = bins[i];
             float rawVal = db[i];
+            float primaryVal = (i < primaryDb.size()) ? primaryDb[i] : rawVal;
 
             // 1. Update History & Calculate Entropy
             bin.history.push_back(rawVal);
@@ -110,24 +113,27 @@ void VisualizerWidget::updateData(const std::vector<AudioEngine::FrameData>& dat
             float p = 3.0f * m_entropyStrength;
             float responsiveness = 0.02f + (0.98f * std::pow(order, p));
 
-            // 3. Update Visual Bar Height (Slew Limiting)
+            // 3. Update Visual Bar Height (Mixed Signal)
             bin.visualDb = (bin.visualDb * (1.0f - responsiveness)) + (rawVal * responsiveness);
 
-            // 4. Trail Physics
+            // 4. Update Primary Visual DB (Steady Signal for Pattern)
+            // Use a fixed, slower responsiveness for the pattern to keep it stable
+            float patternResp = 0.1f; 
+            bin.primaryVisualDb = (bin.primaryVisualDb * (1.0f - patternResp)) + (primaryVal * patternResp);
+
+            // 5. Trail Physics
             bin.trailLife = std::max(bin.trailLife - 0.02f, 0.0f);
             
             float flux = rawVal - bin.lastRawDb;
             bin.lastRawDb = rawVal;
 
             if (flux > 0) {
-                // Impact Multiplier: Boosts the "Kick" of ordered signals
                 float impactMultiplier = 1.0f + (3.0f * order * m_entropyStrength);
                 float jumpTarget = bin.visualDb + (flux * impactMultiplier);
                 
                 if (jumpTarget > bin.trailDb) {
                     bin.trailDb = jumpTarget;
                     bin.trailLife = 1.0f;
-                    // Thickness based on order
                     bin.trailThickness = 1.0f + (order * 4.0f);
                 }
             }
@@ -268,10 +274,11 @@ void VisualizerWidget::drawContent(QPainter& p, int w, int h) {
         const auto& bins = m_channels[ch].bins;
         if (bins.empty()) continue;
 
-        // 1. Calculate Raw Energy
+        // 1. Calculate Raw Energy (Using Primary DB for Pattern)
         std::vector<float> rawEnergy(bins.size());
         for (size_t j = 0; j < bins.size(); ++j) {
-             rawEnergy[j] = std::clamp((bins[j].visualDb + 80.0f) / 80.0f, 0.0f, 1.0f);
+             // Use primaryVisualDb for the pattern calculation to keep it stable
+             rawEnergy[j] = std::clamp((bins[j].primaryVisualDb + 80.0f) / 80.0f, 0.0f, 1.0f);
         }
 
         // 2. Auto-Balance Highs vs Lows (Dynamic Normalization)
@@ -297,21 +304,19 @@ void VisualizerWidget::drawContent(QPainter& p, int w, int h) {
                 val *= boost;
             }
             
-            // Soft-Knee Compression
             float compressed = std::tanh(val);
             vertexEnergy[j] = compressed;
             if (compressed > globalMax) globalMax = compressed;
         }
 
-        // 3. Global Normalization (Ensure peaks hit 1.0)
+        // 3. Global Normalization
-        // This fixes the "faded highs" issue by ensuring the loudest part of the spectrum
-        // (even if it's just high hats) is fully opaque.
         for (float& v : vertexEnergy) {
             v = std::clamp(v / globalMax, 0.0f, 1.0f);
         }
 
-        // 4. Calculate Segment Modifiers (Vote System)
+        // 4. Calculate Segment Modifiers (Procedural Pattern with Competition)
-        std::vector<float> segmentModifiers(freqs.size() - 1, 0.0f);
+        std::vector<float> brightMods(freqs.size() - 1, 0.0f);
+        std::vector<float> alphaMods(freqs.size() - 1, 0.0f);
 
         for (size_t i = 1; i < vertexEnergy.size() - 1; ++i) {
             float curr = vertexEnergy[i];
@@ -322,29 +327,54 @@ void VisualizerWidget::drawContent(QPainter& p, int w, int h) {
             if (curr > prev && curr > next) {
                 bool leftDominant = (prev > next);
                 float sharpness = std::min(curr - prev, curr - next);
-                float intensity = std::clamp(sharpness * 4.0f, 0.0f, 1.0f);
                 
-                float brightBoost = 1.0f * intensity;
+                // COMPETITION LOGIC (Controlled by Entropy)
+                // m_entropyStrength (0.0 to 3.0)
+                // Low Entropy = Smoother, Less Aggressive
+                // High Entropy = Sharper, More Faceted
                 
-                if (leftDominant) {
+                float entropyFactor = std::max(0.1f, m_entropyStrength);
-                    // Left Segment (i-1) gets Brightness
-                    if (i > 0) segmentModifiers[i-1] += brightBoost;
                 
-                    // Right Side (Dark Side)
+                // Aggressive Contrast: Boost low sharpness
-                    if (i < segmentModifiers.size())     segmentModifiers[i]   -= 0.6f * intensity; // Dark
+                float peakIntensity = std::clamp(std::pow(sharpness * 10.0f * entropyFactor, 0.3f), 0.0f, 1.0f);
-                    if (i + 1 < segmentModifiers.size()) segmentModifiers[i+1] -= 0.9f * intensity; // Darker
-                    if (i + 2 < segmentModifiers.size()) segmentModifiers[i+2] -= 0.6f * intensity; // Dark
-                    if (i + 3 < segmentModifiers.size()) segmentModifiers[i+3] -= 0.3f * intensity; // Fade
 
-                } else {
+                // Dynamic Decay: 
-                    // Right Segment (i) gets Brightness
+                float decayBase = 0.65f - std::clamp(sharpness * 3.0f * entropyFactor, 0.0f, 0.35f);
-                    if (i < segmentModifiers.size()) segmentModifiers[i] += brightBoost;
 
-                    // Left Side (Dark Side)
+                auto applyPattern = [&](int dist, bool isBrightSide, int direction) {
-                    if (i > 0) segmentModifiers[i-1] -= 0.6f * intensity; // Dark
+                    int segIdx = (direction == -1) ? (i - dist) : (i + dist - 1);
-                    if (i >= 2) segmentModifiers[i-2] -= 0.9f * intensity; // Darker
+                    if (segIdx < 0 || segIdx >= (int)brightMods.size()) return;
-                    if (i >= 3) segmentModifiers[i-3] -= 0.6f * intensity; // Dark
+
-                    if (i >= 4) segmentModifiers[i-4] -= 0.3f * intensity; // Fade
+                    int cycle = (dist - 1) / 3;
+                    int step = (dist - 1) % 3;
+                    
+                    float decay = std::pow(decayBase, cycle);
+                    float intensity = peakIntensity * decay;
+                    
+                    if (intensity < 0.01f) return;
+
+                    int type = step;
+                    if (isBrightSide) type = (type + 2) % 3;
+
+                    switch (type) {
+                        case 0: // Ghost (Bright + Trans)
+                            brightMods[segIdx] += 0.8f * intensity;
+                            alphaMods[segIdx]  -= 0.8f * intensity;
+                            break;
+                        case 1: // Shadow (Dark + Opaque)
+                            brightMods[segIdx] -= 0.8f * intensity;
+                            alphaMods[segIdx]  += 0.2f * intensity;
+                            break;
+                        case 2: // Highlight (Bright + Opaque)
+                            brightMods[segIdx] += 0.8f * intensity;
+                            alphaMods[segIdx]  += 0.2f * intensity;
+                            break;
+                    }
+                };
+
+                for (int d = 1; d <= 12; ++d) {
+                    applyPattern(d, leftDominant, -1);
+                    applyPattern(d, !leftDominant, 1);
                 }
             }
         }
@@ -369,15 +399,14 @@ void VisualizerWidget::drawContent(QPainter& p, int w, int h) {
                 binColor = QColor::fromHsvF(hue, 1.0f, 1.0f);
             }
 
-            // Base Brightness from Energy
+            // Base Brightness from Energy (Using Primary for stability)
             float avgEnergy = (vertexEnergy[i] + vertexEnergy[i+1]) / 2.0f;
             float baseBrightness = std::pow(avgEnergy, 0.5f); 
 
-            // Apply Shadow/Highlight Modifiers to Brightness
+            // Apply Brightness Modifier
-            float modifier = segmentModifiers[i];
+            float bMod = brightMods[i];
-            float multiplier = (modifier >= 0) ? (1.0f + modifier) : (1.0f / (1.0f - modifier * 2.0f));
+            float bMult = (bMod >= 0) ? (1.0f + bMod) : (1.0f / (1.0f - bMod * 2.0f));
-            
+            float finalBrightness = std::clamp(baseBrightness * bMult * m_brightness, 0.0f, 1.0f);
-            float finalBrightness = std::clamp(baseBrightness * multiplier * m_brightness, 0.0f, 1.0f);
 
             float h_val, s, v, a;
             binColor.getHsvF(&h_val, &s, &v, &a);
@@ -395,18 +424,14 @@ void VisualizerWidget::drawContent(QPainter& p, int w, int h) {
                 lineColor = dynamicBinColor;
             }
 
-            // --- Alpha Calculation with Shadow Logic ---
+            // Apply Alpha Modifier
-            // 1. Start with the boosted, normalized energy (avgEnergy)
+            float aMod = alphaMods[i];
-            // 2. Apply the shadow modifier (linear scalar for alpha)
+            float aMult = (aMod >= 0) ? (1.0f + aMod * 0.5f) : (1.0f + aMod); 
-            //    If modifier is negative (shadow), reduce alpha.
+            if (aMult < 0.1f) aMult = 0.1f; 
-            //    If positive (peak), boost alpha slightly.
-            float alphaScalar = (modifier >= 0) ? (1.0f + modifier * 0.2f) : (1.0f + modifier);
 
             float intensity = avgEnergy; 
             float alpha = 0.4f + (intensity - 0.5f) * m_contrast;
-            
+            alpha = std::clamp(alpha * aMult, 0.0f, 1.0f);
-            // Apply the shadow scalar to the calculated alpha
-            alpha = std::clamp(alpha * alphaScalar, 0.0f, 1.0f);
             
             fillColor.setAlphaF(alpha);
             lineColor.setAlphaF(0.9f);
@@ -422,6 +447,7 @@ void VisualizerWidget::drawContent(QPainter& p, int w, int h) {
             float x1 = getX(freqs[i] * xOffset) * w;
             float x2 = getX(freqs[i+1] * xOffset) * w;
 
+            // Use visualDb (Mixed) for Height
             float barH1 = std::clamp((b.visualDb + 80.0f) / 80.0f * h, 0.0f, (float)h);
             float barH2 = std::clamp((bNext.visualDb + 80.0f) / 80.0f * h, 0.0f, (float)h);
 

src/VisualizerWidget.h

@@ -29,7 +29,8 @@ private:
 
     struct BinState {
         std::deque<float> history; // For entropy calculation
-        float visualDb = -100.0f;  // The dampened/processed value for display
+        float visualDb = -100.0f;  // Mixed (Height)
+        float primaryVisualDb = -100.0f; // Primary (Pattern)
         float lastRawDb = -100.0f; // To calculate flux
         
         // Trail Physics

src/complex_block.cpp

@@ -0,0 +1,102 @@
+/**       Filename: complex_block.cpp
+ *        Location: /src/
+**/
+/**
+ **********************************************************************************
+ *
+ *                                 Description
+ *
+ **********************************************************************************
+ *
+ *  Implements the `BlockHilbert` class defined in `complex_block.h`.
+ *
+ *  This implementation uses the FFTW3 library to perform the necessary
+ *  Fast Fourier Transform, apply the Hilbert filter in the frequency domain,
+ *  and then perform an Inverse FFT to obtain the analytic signal. All FFTW
+ *  resources are allocated and freed within the single function call.
+ *
+**/
+/**
+ **********************************************************************************
+ *
+ *                                 Implements
+ *
+ **********************************************************************************
+ *
+ *  BlockHilbert::hilbertTransform(const std::vector<double>&, const std::vector<double>&)
+ *
+**/
+
+#include "complex_block.h"
+#include <fftw3.h>
+#include <stdexcept>
+#include <cmath>
+
+std::pair<std::vector<std::complex<double>>, std::vector<std::complex<double>>>
+BlockHilbert::hilbertTransform(const std::vector<double>& left_signal, const std::vector<double>& right_signal) {
+    size_t n = left_signal.size();
+    if (n == 0 || n != right_signal.size()) {
+        return {{}, {}};
+    }
+
+    fftw_complex* fft_in_L = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+    fftw_complex* fft_out_L = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+    fftw_complex* ifft_in_L = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+    fftw_complex* ifft_out_L = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+
+    fftw_complex* fft_in_R = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+    fftw_complex* fft_out_R = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+    fftw_complex* ifft_in_R = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+    fftw_complex* ifft_out_R = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+
+    if (!fft_in_L || !fft_out_L || !ifft_in_L || !ifft_out_L || !fft_in_R || !fft_out_R || !ifft_in_R || !ifft_out_R) {
+        throw std::runtime_error("FFTW memory allocation failed in BlockHilbert.");
+    }
+
+    fftw_plan plan_forward_L = fftw_plan_dft_1d(static_cast<int>(n), fft_in_L, fft_out_L, FFTW_FORWARD, FFTW_ESTIMATE);
+    fftw_plan plan_backward_L = fftw_plan_dft_1d(static_cast<int>(n), ifft_in_L, ifft_out_L, FFTW_BACKWARD, FFTW_ESTIMATE);
+    fftw_plan plan_forward_R = fftw_plan_dft_1d(static_cast<int>(n), fft_in_R, fft_out_R, FFTW_FORWARD, FFTW_ESTIMATE);
+    fftw_plan plan_backward_R = fftw_plan_dft_1d(static_cast<int>(n), ifft_in_R, ifft_out_R, FFTW_BACKWARD, FFTW_ESTIMATE);
+
+    if (!plan_forward_L || !plan_backward_L || !plan_forward_R || !plan_backward_R) {
+        fftw_free(fft_in_L); fftw_free(fft_out_L); fftw_free(ifft_in_L); fftw_free(ifft_out_L);
+        fftw_free(fft_in_R); fftw_free(fft_out_R); fftw_free(ifft_in_R); fftw_free(ifft_out_R);
+        throw std::runtime_error("FFTW plan creation failed in BlockHilbert.");
+    }
+
+    for (size_t i = 0; i < n; ++i) {
+        fft_in_L[i][0] = left_signal[i]; fft_in_L[i][1] = 0.0;
+        fft_in_R[i][0] = right_signal[i]; fft_in_R[i][1] = 0.0;
+    }
+
+    fftw_execute(plan_forward_L);
+    fftw_execute(plan_forward_R);
+
+    for (size_t i = 0; i < n; ++i) {
+        double multiplier = 1.0;
+        if (i > 0 && i < n / 2.0) { multiplier = 2.0; } 
+        else if (i > n / 2.0) { multiplier = 0.0; }
+        
+        ifft_in_L[i][0] = fft_out_L[i][0] * multiplier; ifft_in_L[i][1] = fft_out_L[i][1] * multiplier;
+        ifft_in_R[i][0] = fft_out_R[i][0] * multiplier; ifft_in_R[i][1] = fft_out_R[i][1] * multiplier;
+    }
+
+    fftw_execute(plan_backward_L);
+    fftw_execute(plan_backward_R);
+
+    std::vector<std::complex<double>> analytic_L(n);
+    std::vector<std::complex<double>> analytic_R(n);
+    for (size_t i = 0; i < n; ++i) {
+        analytic_L[i].real(ifft_out_L[i][0] / static_cast<double>(n));
+        analytic_L[i].imag(ifft_out_L[i][1] / static_cast<double>(n));
+        analytic_R[i].real(ifft_out_R[i][0] / static_cast<double>(n));
+        analytic_R[i].imag(ifft_out_R[i][1] / static_cast<double>(n));
+    }
+
+    fftw_destroy_plan(plan_forward_L); fftw_destroy_plan(plan_backward_L);
+    fftw_destroy_plan(plan_forward_R); fftw_destroy_plan(plan_backward_R);
+    fftw_free(fft_in_L); fftw_free(fft_out_L); fftw_free(ifft_in_L); fftw_free(ifft_out_L);
+    fftw_free(fft_in_R); fftw_free(fft_out_R); fftw_free(ifft_in_R); fftw_free(ifft_out_R);
+
+    return {analytic_L, analytic_R};
+}

src/complex_block.h

@@ -0,0 +1,24 @@
+/** inc/complex_block.h
+ * 
+
+ *  This class performs a Hilbert transform on an entire data
+ *  block at once. It is stateless and suitable for offline analysis or any
+ *  task where the entire buffer is available and processing latency is not a
+ *  critical concern.
+
+**/
+
+#ifndef COMPLEX_BLOCK_H
+#define COMPLEX_BLOCK_H
+
+#include <vector>
+#include <complex>
+#include <utility> // For std::pair
+
+class BlockHilbert {
+public:
+    std::pair<std::vector<std::complex<double>>, std::vector<std::complex<double>>> 
+    hilbertTransform(const std::vector<double>& left_signal, const std::vector<double>& right_signal);
+};
+
+#endif // COMPLEX_BLOCK_H

src/complex_frames.cpp

@@ -0,0 +1,177 @@
+/** Filename: complex_frames.cpp
+ * Location: /src/
+**/
+/**
+ **********************************************************************************
+ *
+ * Description
+ *
+ **********************************************************************************
+ *
+ * Implements the `RealtimeHilbert` class for streaming audio processing.
+ *
+ * The core logic involves:
+ * 1. `reinit`: Allocating all necessary FFTW plans and memory buffers once for both channels.
+ * 2. `process`:
+ * - Shifting history buffers for L/R to accommodate the new input frames.
+ * - Performing a real-to-complex FFT on the full window for both channels.
+ * - Applying the Hilbert filter in the frequency domain for both channels.
+ * - Performing a complex-to-complex IFFT to get the analytic signal for both channels.
+ * - Extracting only the valid, new frames of data from the end of the
+ *   overlap-add buffers.
+ * 3. Destructor/cleanup: Safely freeing all allocated FFTW resources.
+ *
+**/
+/**
+ **********************************************************************************
+ *
+ * Implements
+ *
+ **********************************************************************************
+ *
+ * All methods of the `RealtimeHilbert` class.
+ *
+**/
+
+#include "complex_frames.h"
+#include <stdexcept>
+#include <cstring>
+
+RealtimeHilbert::RealtimeHilbert() :
+    m_fft_size(0), m_hop_size(0),
+    m_plan_forward_L(nullptr), m_plan_backward_L(nullptr),
+    m_fft_input_real_L(nullptr), m_fft_output_spectrum_L(nullptr),
+    m_ifft_input_spectrum_L(nullptr), m_ifft_output_complex_L(nullptr),
+    m_plan_forward_R(nullptr), m_plan_backward_R(nullptr),
+    m_fft_input_real_R(nullptr), m_fft_output_spectrum_R(nullptr),
+    m_ifft_input_spectrum_R(nullptr), m_ifft_output_complex_R(nullptr)
+{}
+
+RealtimeHilbert::~RealtimeHilbert() {
+    cleanup();
+}
+
+void RealtimeHilbert::cleanup() {
+    if (m_plan_forward_L) fftw_destroy_plan(m_plan_forward_L);
+    if (m_plan_backward_L) fftw_destroy_plan(m_plan_backward_L);
+    if (m_fft_input_real_L) fftw_free(m_fft_input_real_L);
+    if (m_fft_output_spectrum_L) fftw_free(m_fft_output_spectrum_L);
+    if (m_ifft_input_spectrum_L) fftw_free(m_ifft_input_spectrum_L);
+    if (m_ifft_output_complex_L) fftw_free(m_ifft_output_complex_L);
+
+    if (m_plan_forward_R) fftw_destroy_plan(m_plan_forward_R);
+    if (m_plan_backward_R) fftw_destroy_plan(m_plan_backward_R);
+    if (m_fft_input_real_R) fftw_free(m_fft_input_real_R);
+    if (m_fft_output_spectrum_R) fftw_free(m_fft_output_spectrum_R);
+    if (m_ifft_input_spectrum_R) fftw_free(m_ifft_input_spectrum_R);
+    if (m_ifft_output_complex_R) fftw_free(m_ifft_output_complex_R);
+
+    m_plan_forward_L = nullptr; m_plan_backward_L = nullptr;
+    m_fft_input_real_L = nullptr; m_fft_output_spectrum_L = nullptr;
+    m_ifft_input_spectrum_L = nullptr; m_ifft_output_complex_L = nullptr;
+    
+    m_plan_forward_R = nullptr; m_plan_backward_R = nullptr;
+    m_fft_input_real_R = nullptr; m_fft_output_spectrum_R = nullptr;
+    m_ifft_input_spectrum_R = nullptr; m_ifft_output_complex_R = nullptr;
+}
+
+void RealtimeHilbert::reinit(size_t fft_size) {
+    cleanup();
+    m_fft_size = fft_size;
+    m_hop_size = 0;
+
+    if (m_fft_size == 0 || (m_fft_size & (m_fft_size - 1)) != 0) {
+        throw std::invalid_argument("FFT size must be a power of two.");
+    }
+
+    // Allocate for Left Channel
+    m_fft_input_real_L = fftw_alloc_real(m_fft_size);
+    m_fft_output_spectrum_L = fftw_alloc_complex(m_fft_size / 2 + 1);
+    m_ifft_input_spectrum_L = fftw_alloc_complex(m_fft_size);
+    m_ifft_output_complex_L = fftw_alloc_complex(m_fft_size);
+    
+    // Allocate for Right Channel
+    m_fft_input_real_R = fftw_alloc_real(m_fft_size);
+    m_fft_output_spectrum_R = fftw_alloc_complex(m_fft_size / 2 + 1);
+    m_ifft_input_spectrum_R = fftw_alloc_complex(m_fft_size);
+    m_ifft_output_complex_R = fftw_alloc_complex(m_fft_size);
+
+    if (!m_fft_input_real_L || !m_fft_output_spectrum_L || !m_ifft_input_spectrum_L || !m_ifft_output_complex_L ||
+        !m_fft_input_real_R || !m_fft_output_spectrum_R || !m_ifft_input_spectrum_R || !m_ifft_output_complex_R) {
+        cleanup();
+        throw std::runtime_error("Failed to allocate FFTW buffers in RealtimeHilbert.");
+    }
+
+    // Create plans for Left Channel
+    m_plan_forward_L = fftw_plan_dft_r2c_1d(m_fft_size, m_fft_input_real_L, m_fft_output_spectrum_L, FFTW_ESTIMATE);
+    m_plan_backward_L = fftw_plan_dft_1d(m_fft_size, m_ifft_input_spectrum_L, m_ifft_output_complex_L, FFTW_BACKWARD, FFTW_ESTIMATE);
+
+    // Create plans for Right Channel
+    m_plan_forward_R = fftw_plan_dft_r2c_1d(m_fft_size, m_fft_input_real_R, m_fft_output_spectrum_R, FFTW_ESTIMATE);
+    m_plan_backward_R = fftw_plan_dft_1d(m_fft_size, m_ifft_input_spectrum_R, m_ifft_output_complex_R, FFTW_BACKWARD, FFTW_ESTIMATE);
+
+    if (!m_plan_forward_L || !m_plan_backward_L || !m_plan_forward_R || !m_plan_backward_R) {
+        cleanup();
+        throw std::runtime_error("Failed to create FFTW plans in RealtimeHilbert.");
+    }
+    m_history_buffer_L.assign(m_fft_size, 0.0);
+    m_history_buffer_R.assign(m_fft_size, 0.0);
+}
+
+std::pair<std::vector<std::complex<double>>, std::vector<std::complex<double>>>
+RealtimeHilbert::process(const std::vector<double>& left_input_block, const std::vector<double>& right_input_block) {
+    if (m_fft_size == 0) {
+        throw std::runtime_error("Processor has not been initialized. Call reinit() first.");
+    }
+    if (m_hop_size == 0) {
+        m_hop_size = left_input_block.size();
+        if (m_hop_size > m_fft_size) {
+            throw std::invalid_argument("Input block size cannot be larger than FFT size.");
+        }
+    } else if (left_input_block.size() != m_hop_size || right_input_block.size() != m_hop_size) {
+        throw std::runtime_error("Input block size changed. Processor must be re-initialized.");
+    }
+
+    // --- Left Channel Processing ---
+    std::memmove(m_history_buffer_L.data(), m_history_buffer_L.data() + m_hop_size, (m_fft_size - m_hop_size) * sizeof(double));
+    std::memcpy(m_history_buffer_L.data() + (m_fft_size - m_hop_size), left_input_block.data(), m_hop_size * sizeof(double));
+    std::memcpy(m_fft_input_real_L, m_history_buffer_L.data(), m_fft_size * sizeof(double));
+    fftw_execute(m_plan_forward_L);
+
+    // --- Right Channel Processing ---
+    std::memmove(m_history_buffer_R.data(), m_history_buffer_R.data() + m_hop_size, (m_fft_size - m_hop_size) * sizeof(double));
+    std::memcpy(m_history_buffer_R.data() + (m_fft_size - m_hop_size), right_input_block.data(), m_hop_size * sizeof(double));
+    std::memcpy(m_fft_input_real_R, m_history_buffer_R.data(), m_fft_size * sizeof(double));
+    fftw_execute(m_plan_forward_R);
+
+    // --- Hilbert Transform in Frequency Domain (for both channels) ---
+    for (size_t i = 0; i < m_fft_size; ++i) {
+        if (i < m_fft_size / 2 + 1) {
+            double multiplier = 1.0;
+            if (i > 0 && i < m_fft_size / 2.0) multiplier = 2.0;
+            
+            m_ifft_input_spectrum_L[i][0] = multiplier * m_fft_output_spectrum_L[i][0];
+            m_ifft_input_spectrum_L[i][1] = multiplier * m_fft_output_spectrum_L[i][1];
+            m_ifft_input_spectrum_R[i][0] = multiplier * m_fft_output_spectrum_R[i][0];
+            m_ifft_input_spectrum_R[i][1] = multiplier * m_fft_output_spectrum_R[i][1];
+        } else {
+            m_ifft_input_spectrum_L[i][0] = 0.0; m_ifft_input_spectrum_L[i][1] = 0.0;
+            m_ifft_input_spectrum_R[i][0] = 0.0; m_ifft_input_spectrum_R[i][1] = 0.0;
+        }
+    }
+    
+    fftw_execute(m_plan_backward_L);
+    fftw_execute(m_plan_backward_R);
+
+    std::vector<std::complex<double>> result_L(m_hop_size);
+    std::vector<std::complex<double>> result_R(m_hop_size);
+    for (size_t i = 0; i < m_hop_size; ++i) {
+        size_t index = (m_fft_size - m_hop_size) + i;
+        result_L[i].real(m_ifft_output_complex_L[index][0] / m_fft_size);
+        result_L[i].imag(m_ifft_output_complex_L[index][1] / m_fft_size);
+        result_R[i].real(m_ifft_output_complex_R[index][0] / m_fft_size);
+        result_R[i].imag(m_ifft_output_complex_R[index][1] / m_fft_size);
+    }
+
+    return {result_L, result_R};
+}

src/complex_frames.h

@@ -0,0 +1,56 @@
+/** inc/complex_frames.h
+ 
+  Defines the interface for the `RealtimeHilbert` class.
+ 
+    This is a stateful processor designed for continuous, low-latency, real-time
+    audio streaming. It uses an internal history buffer and an overlap-add to prevent
+    artifacts that would occur with block-based processing of the same chunk-sizes.
+
+    The point of this is to provide a complex stream alongside the raw real stream simultaneously.
+    It's kind of the core of this whole engine. This is what allowed the spiral visualizer idea to work in real-time.
+    
+**/
+
+#ifndef COMPLEX_FRAMES_H
+#define COMPLEX_FRAMES_H
+
+#include <vector>
+#include <complex>
+#include <utility>
+#include <fftw3.h>
+
+class RealtimeHilbert {
+public:
+    RealtimeHilbert();
+    ~RealtimeHilbert();
+
+    void reinit(size_t fft_size);
+    std::pair<std::vector<std::complex<double>>, std::vector<std::complex<double>>> 
+    process(const std::vector<double>& left_input_block, const std::vector<double>& right_input_block);
+
+private:
+    void cleanup();
+
+    size_t m_fft_size;
+    size_t m_hop_size;
+    std::vector<double> m_history_buffer_L;
+    std::vector<double> m_history_buffer_R;
+
+    // --- Left Channel Resources ---
+    fftw_plan m_plan_forward_L;
+    fftw_plan m_plan_backward_L;
+    double*       m_fft_input_real_L;
+    fftw_complex* m_fft_output_spectrum_L;
+    fftw_complex* m_ifft_input_spectrum_L;
+    fftw_complex* m_ifft_output_complex_L;
+
+    // --- Right Channel Resources ---
+    fftw_plan m_plan_forward_R;
+    fftw_plan m_plan_backward_R;
+    double*       m_fft_input_real_R;
+    fftw_complex* m_fft_output_spectrum_R;
+    fftw_complex* m_ifft_input_spectrum_R;
+    fftw_complex* m_ifft_output_complex_R;
+};
+
+#endif // COMPLEX_FRAMES_H

src/trig_interpolation.cpp

@@ -0,0 +1,453 @@
+/** Filename: trig_interpolation.cpp
+ * Location: /src/
+ *
+ * Description: Implements the core trigonometric interpolation logic.
+ * This is where our theoretical discussion becomes code.
+**/
+#include "trig_interpolation.h"
+#include "complex_block.h" // For the Hilbert Transform
+#include <fftw3.h>
+#include <cmath>
+#include <numeric>
+#include <algorithm>
+#include <stdexcept>
+#include <functional> // For std::function
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+// Forward declaration for the new helper function
+static std::vector<double> create_match_curve(
+    const std::vector<double>& smooth_source,
+    const std::vector<double>& smooth_target,
+    double strength
+);
+
+// --- FIX: Add new static helper function to correctly process magnitude spectrums ---
+static std::vector<std::complex<double>> get_cepstrum_from_magnitude_spectrum(
+    const std::vector<double>& mag_spectrum,
+    BlockHilbert& hilbert_processor,
+    std::function<void(const std::vector<std::complex<double>>&, std::vector<std::complex<double>>&)> ifft_func
+) {
+    if (mag_spectrum.empty()) return {};
+
+    // The input mag_spectrum is a half-spectrum of size (N/2 + 1).
+    size_t half_n = mag_spectrum.size();
+    size_t n = (half_n > 1) ? (half_n - 1) * 2 : 0;
+    if (n == 0) return {};
+
+    // 1. Take the natural logarithm of the magnitude spectrum.
+    std::vector<std::complex<double>> log_mag_full(n, {0.0, 0.0});
+    for (size_t i = 0; i < half_n; ++i) {
+        log_mag_full[i] = { log(std::max(1e-20, mag_spectrum[i])), 0.0 };
+    }
+    // Create the conjugate symmetric part for the IFFT.
+    for (size_t i = 1; i < half_n - 1; ++i) {
+        log_mag_full[n - i] = log_mag_full[i];
+    }
+
+    // 2. IFFT to get the real cepstrum.
+    std::vector<std::complex<double>> cepstrum(n);
+    ifft_func(log_mag_full, cepstrum);
+
+    // 3. Perform Hilbert transform substitution on the real part to get the analytic cepstrum.
+    std::vector<double> real_part(n);
+    for(size_t i = 0; i < n; ++i) real_part[i] = cepstrum[i].real();
+
+    auto substituted_pair = hilbert_processor.hilbertTransform(real_part, real_part);
+
+    return substituted_pair.first;
+}
+
+
+// --- Main Public Method ---
+
+auto TrigInterpolation::process_and_generate_fir(
+    const std::vector<double>& source_L_in, const std::vector<double>& source_R_in,
+    const std::vector<double>& target_L_in, const std::vector<double>& target_R_in,
+    int granularity, int detail_level, double strength, int fir_size, double lr_link,
+    PhaseMode phase_mode
+) -> std::tuple<
+    std::vector<float>, std::vector<float>, // FIR Taps (L/R)
+    std::vector<double>,                    // Frequency Axis
+    std::vector<double>, std::vector<double>, std::vector<double>, // L(src,tgt,match)
+    std::vector<double>, std::vector<double>, std::vector<double>  // R(src,tgt,match)
+> {
+
+    // --- FIX: Call the new helper function instead of get_idealized_cepstrum ---
+    // The lambda captures `this` to allow calling the member function `ifft`.
+    auto ifft_lambda = [this](const auto& in, auto& out) { this->ifft(in, out); };
+    auto cepstrum_src_L = get_cepstrum_from_magnitude_spectrum(source_L_in, m_block_hilbert, ifft_lambda);
+    auto cepstrum_src_R = get_cepstrum_from_magnitude_spectrum(source_R_in, m_block_hilbert, ifft_lambda);
+    auto cepstrum_tgt_L = get_cepstrum_from_magnitude_spectrum(target_L_in, m_block_hilbert, ifft_lambda);
+    auto cepstrum_tgt_R = get_cepstrum_from_magnitude_spectrum(target_R_in, m_block_hilbert, ifft_lambda);
+
+    // 2. Independently Smooth Each Curve
+    auto smooth_src_L = idealize_curve(cepstrum_src_L, granularity, detail_level);
+    auto smooth_src_R = idealize_curve(cepstrum_src_R, granularity, detail_level);
+    auto smooth_tgt_L = idealize_curve(cepstrum_tgt_L, granularity, detail_level);
+    auto smooth_tgt_R = idealize_curve(cepstrum_tgt_R, granularity, detail_level);
+
+    // 3. Perform the second trigonometric interpolation to get the match curves
+    auto pre_link_match_L = create_match_curve(smooth_src_L, smooth_tgt_L, strength);
+    auto pre_link_match_R = create_match_curve(smooth_src_R, smooth_tgt_R, strength);
+    
+    std::vector<double> match_curve_L(pre_link_match_L.size());
+    std::vector<double> match_curve_R(pre_link_match_R.size());
+
+    // 4. Apply Stereo Link to the final match curves
+    for (size_t i = 0; i < pre_link_match_L.size(); ++i) {
+        double val_L = pre_link_match_L[i];
+        double val_R = pre_link_match_R[i];
+        match_curve_L[i] = val_L * (1.0 - lr_link) + (val_L + val_R) * 0.5 * lr_link;
+        match_curve_R[i] = val_R * (1.0 - lr_link) + (val_L + val_R) * 0.5 * lr_link;
+    }
+    
+    // 5. Generate FIR Filter from Match Curve
+    auto fir_coeffs_L = create_fir_from_curve(match_curve_L, fir_size, phase_mode);
+    auto fir_coeffs_R = create_fir_from_curve(match_curve_R, fir_size, phase_mode);
+
+    // 6. Prepare Plot Data
+    std::vector<double> frequency_axis, source_L_db, target_L_db, matched_L_db,
+                          source_R_db, target_R_db, matched_R_db;
+
+    size_t num_plot_points = match_curve_L.size() > 1 ? match_curve_L.size() / 2 : 0;
+    if (num_plot_points == 0) {
+        // Return empty tuple if no data
+        return { {}, {}, {}, {}, {}, {}, {}, {}, {} };
+    }
+
+    // The frequency axis should correspond to the original half-spectrum size
+    size_t freq_axis_size = source_L_in.size();
+    frequency_axis.resize(freq_axis_size);
+    std::iota(frequency_axis.begin(), frequency_axis.end(), 0.0);
+
+    auto to_db = [](double val) { return 20.0 * log10(std::max(1e-9, val)); };
+
+    source_L_db.resize(freq_axis_size);
+    target_L_db.resize(freq_axis_size);
+    matched_L_db.resize(freq_axis_size);
+    source_R_db.resize(freq_axis_size);
+    target_R_db.resize(freq_axis_size);
+    matched_R_db.resize(freq_axis_size);
+
+    for(size_t i = 0; i < freq_axis_size; ++i) {
+        source_L_db[i] = to_db(smooth_src_L[i]);
+        target_L_db[i] = to_db(smooth_tgt_L[i]);
+        matched_L_db[i] = to_db(match_curve_L[i]);
+        source_R_db[i] = to_db(smooth_src_R[i]);
+        target_R_db[i] = to_db(smooth_tgt_R[i]);
+        matched_R_db[i] = to_db(match_curve_R[i]);
+    }
+
+    return {
+        fir_coeffs_L, fir_coeffs_R,
+        frequency_axis,
+        source_L_db, target_L_db, matched_L_db,
+        source_R_db, target_R_db, matched_R_db
+    };
+}
+
+// --- New Function for Match Curve Interpolation ---
+static std::vector<double> create_match_curve(
+    const std::vector<double>& smooth_source,
+    const std::vector<double>& smooth_target,
+    double strength) {
+    
+    if (smooth_source.empty()) return {};
+    size_t n = smooth_source.size();
+
+    // 1. Find extrema in the source curve to use as anchors
+    std::vector<std::pair<size_t, double>> source_extrema;
+    if (n > 2) {
+        for (size_t i = 1; i < n - 1; ++i) {
+            if ((smooth_source[i] > smooth_source[i-1] && smooth_source[i] > smooth_source[i+1]) ||
+                (smooth_source[i] < smooth_source[i-1] && smooth_source[i] < smooth_source[i+1])) {
+                source_extrema.push_back({i, smooth_source[i]});
+            }
+        }
+    }
+    if (source_extrema.empty()) return smooth_source; // Return original if no features found
+
+    // 2. Create a new set of "match" extrema by interpolating towards the target
+    std::vector<std::pair<size_t, double>> match_extrema;
+    for (const auto& extremum : source_extrema) {
+        size_t index = extremum.first;
+        double source_val = extremum.second;
+        double target_val = smooth_target[index];
+        double match_val = source_val + (target_val - source_val) * strength;
+        match_extrema.push_back({index, match_val});
+    }
+
+    // 3. Reconstruct the final curve from the new match extrema using cosine bell interpolation
+    std::vector<double> final_match_curve(n);
+    for (size_t i = 0; i < n; ++i) {
+        double total_weight = 0.0;
+        double weighted_sum = 0.0;
+        for (size_t j = 0; j < match_extrema.size(); ++j) {
+            size_t extremum_idx = match_extrema[j].first;
+            double dist = static_cast<double>(i) - extremum_idx;
+            
+            // Define the lobe width based on distance to adjacent extrema
+            size_t prev_idx = (j > 0) ? match_extrema[j-1].first : 0;
+            size_t next_idx = (j < match_extrema.size() - 1) ? match_extrema[j+1].first : n - 1;
+            double lobe_width = (static_cast<double>(next_idx - prev_idx)) / 2.0;
+            if (lobe_width == 0) lobe_width = n / (2.0 * match_extrema.size()); // Fallback
+
+            if (std::abs(dist) < lobe_width) {
+                double weight = (cos(dist / lobe_width * M_PI) + 1.0) / 2.0; // Hann window
+                weighted_sum += match_extrema[j].second * weight;
+                total_weight += weight;
+            }
+        }
+        if (total_weight > 0) {
+            final_match_curve[i] = weighted_sum / total_weight;
+        } else {
+             // If a point is not influenced by any lobe, interpolate from nearest extrema
+             auto it = std::lower_bound(match_extrema.begin(), match_extrema.end(), std::make_pair(i, 0.0));
+             if (it == match_extrema.begin()) final_match_curve[i] = it->second;
+             else if (it == match_extrema.end()) final_match_curve[i] = (it-1)->second;
+             else {
+                 auto prev = it - 1;
+                 double progress = static_cast<double>(i - prev->first) / (it->first - prev->first);
+                 final_match_curve[i] = prev->second + (it->second - prev->second) * progress;
+             }
+        }
+    }
+    return final_match_curve;
+}
+
+
+// --- Private Helper Implementations ---
+
+std::vector<float> TrigInterpolation::create_fir_from_curve(const std::vector<double>& match_curve, int fir_size, PhaseMode phase_mode) {
+    if (match_curve.empty() || fir_size <= 0) return {};
+
+    size_t n = match_curve.size();
+    std::vector<std::complex<double>> spectrum(n, {0.0, 0.0});
+
+    // Convert magnitude to complex spectrum (reconstructing phase)
+    if (phase_mode == PhaseMode::Hilbert) {
+        // For minimum phase via Hilbert, we take log, IFFT, causal part, FFT, exp
+        std::vector<std::complex<double>> log_mag(n);
+        for(size_t i = 0; i < n; ++i) log_mag[i] = {log(std::max(1e-20, match_curve[i])), 0.0};
+        
+        std::vector<std::complex<double>> cepstrum(n);
+        ifft(log_mag, cepstrum);
+        
+        // Causal part (zero out negative time)
+        for(size_t i = n / 2; i < n; ++i) cepstrum[i] = {0.0, 0.0};
+        cepstrum[0] *= 1.0; // DC
+        for(size_t i = 1; i < n / 2; ++i) cepstrum[i] *= 2.0; // Positive freqs
+        
+        std::vector<std::complex<double>> min_phase_spectrum(n);
+        fft(cepstrum, min_phase_spectrum);
+        
+        for(size_t i = 0; i < n; ++i) spectrum[i] = std::exp(min_phase_spectrum[i]);
+
+    } else { // Standard Cepstral (results in linear phase)
+        for(size_t i = 0; i < n; ++i) spectrum[i] = {match_curve[i], 0.0};
+    }
+
+    // IFFT to get impulse response
+    std::vector<std::complex<double>> impulse_response(n);
+    ifft(spectrum, impulse_response);
+
+    // Window and extract FIR taps
+    std::vector<float> fir_taps(fir_size);
+    int center = fir_size / 2;
+    for (int i = 0; i < fir_size; ++i) {
+        double hann_window = 0.5 * (1.0 - cos(2.0 * M_PI * i / (fir_size - 1)));
+        int impulse_idx = (i - center + n) % n;
+        fir_taps[i] = static_cast<float>(impulse_response[impulse_idx].real() * hann_window);
+    }
+    
+    // --- FIX: Normalize the filter taps to prevent volume drop ---
+    double tap_sum = 0.0;
+    for (float tap : fir_taps) {
+        tap_sum += tap;
+    }
+    if (std::abs(tap_sum) > 1e-9) {
+        for (float& tap : fir_taps) {
+            tap /= tap_sum;
+        }
+    }
+    
+    return fir_taps;
+}
+
+std::vector<std::complex<double>> TrigInterpolation::get_idealized_cepstrum(const std::vector<double>& signal) {
+    if (signal.empty()) return {};
+    size_t n = signal.size();
+
+    // 1. Hilbert Transform to gain complex input
+    BlockHilbert hilbert;
+    // Use the signal for both L/R as we process mono here
+    auto analytic_pair = hilbert.hilbertTransform(signal, signal);
+    std::vector<std::complex<double>> analytic_signal = analytic_pair.first;
+
+    // 2. Take the FFT
+    std::vector<std::complex<double>> spectrum(n);
+    fft(analytic_signal, spectrum);
+
+    // 3. Take the complex natural logarithm
+    for (auto& val : spectrum) {
+        if (std::abs(val) > 1e-9) {
+            val = std::log(val);
+        } else {
+            val = {0.0, 0.0};
+        }
+    }
+    
+    // 4. Unwrap the phase
+    std::vector<double> unwrapped = unwrap_phase(spectrum);
+    for(size_t i = 0; i < n; ++i) {
+        spectrum[i].imag(unwrapped[i]);
+    }
+    
+    // 5. Perform an IFFT to get complex cepstrum
+    std::vector<std::complex<double>> cepstrum(n);
+    ifft(spectrum, cepstrum);
+
+    // 6. Perform the Hilbert transform substitution
+    std::vector<double> real_part(n);
+    for(size_t i = 0; i < n; ++i) real_part[i] = cepstrum[i].real();
+    
+    auto substituted_pair = hilbert.hilbertTransform(real_part, real_part);
+    
+    return substituted_pair.first; // This is the final analytic cepstrum
+}
+
+// --- Full Implementation of the Idealize Curve Logic ---
+std::vector<double> TrigInterpolation::idealize_curve(const std::vector<std::complex<double>>& analytic_cepstrum, int granularity, int detail_level) {
+    if (analytic_cepstrum.empty()) return {};
+
+    size_t n = analytic_cepstrum.size();
+    std::vector<double> current_curve(n);
+    for(size_t i = 0; i < n; ++i) {
+        current_curve[i] = std::abs(analytic_cepstrum[i]);
+    }
+    
+    // Map slider values to algorithm parameters
+    int num_bins = 4 + static_cast<int>(granularity / 100.0 * 60.0);
+    int iterations = 1 + static_cast<int>(detail_level / 100.0 * 4.0);
+
+    for (int iter = 0; iter < iterations; ++iter) {
+        std::vector<double> next_curve(n, 0.0);
+        std::vector<bool> is_point_set(n, false);
+
+        // 1. Binning
+        std::vector<size_t> bin_boundaries;
+        double log_n = log((double)n);
+        for (int i = 0; i <= num_bins; ++i) {
+            double ratio = static_cast<double>(i) / num_bins;
+            size_t boundary = static_cast<size_t>(exp(ratio * log_n));
+            boundary = std::max((size_t)1, std::min(n - 1, boundary));
+            bin_boundaries.push_back(boundary);
+        }
+        bin_boundaries[0] = 0;
+
+        // 2. Find extrema and process within bins
+        for (int i = 0; i < num_bins; ++i) {
+            size_t bin_start = bin_boundaries[i];
+            size_t bin_end = bin_boundaries[i+1];
+            if (bin_start >= bin_end -1) continue;
+
+            std::vector<std::pair<size_t, double>> extrema;
+            for (size_t j = bin_start + 1; j < bin_end; ++j) {
+                if ((current_curve[j] > current_curve[j-1] && current_curve[j] > current_curve[j+1]) ||
+                    (current_curve[j] < current_curve[j-1] && current_curve[j] < current_curve[j+1])) {
+                    extrema.push_back({j, current_curve[j]});
+                }
+            }
+            if (extrema.empty()) continue;
+
+            // 3. Weaving/Smoothing via cosine bell interpolation
+            for (size_t j = bin_start; j <= bin_end; ++j) {
+                double total_weight = 0.0;
+                double weighted_sum = 0.0;
+                for (const auto& extremum : extrema) {
+                    double dist = static_cast<double>(j) - extremum.first;
+                    double lobe_width = (bin_end - bin_start) / 2.0; // Approximation
+                    if (std::abs(dist) < lobe_width) {
+                        double weight = (cos(dist / lobe_width * M_PI) + 1.0) / 2.0; // Hann window
+                        weighted_sum += extremum.second * weight;
+                        total_weight += weight;
+                    }
+                }
+                if (total_weight > 0) {
+                    next_curve[j] = weighted_sum / total_weight;
+                    is_point_set[j] = true;
+                }
+            }
+        }
+
+        // 4. Gap Filling (linear interpolation)
+        size_t last_set_point = 0;
+        if(is_point_set[0]) last_set_point = 0;
+        else { // Find first set point
+            for(size_t i=0; i<n; ++i) if(is_point_set[i]) { last_set_point = i; break; }
+        }
+
+        for (size_t i = 1; i < n; ++i) {
+            if (is_point_set[i]) {
+                if (i > last_set_point + 1) { // Gap detected
+                    double start_val = next_curve[last_set_point];
+                    double end_val = next_curve[i];
+                    for (size_t j = last_set_point + 1; j < i; ++j) {
+                        double progress = static_cast<double>(j - last_set_point) / (i - last_set_point);
+                        next_curve[j] = start_val + (end_val - start_val) * progress;
+                    }
+                }
+                last_set_point = i;
+            }
+        }
+        current_curve = next_curve;
+    }
+    return current_curve;
+}
+
+std::vector<double> TrigInterpolation::unwrap_phase(const std::vector<std::complex<double>>& c_vec) {
+    std::vector<double> phases(c_vec.size());
+    for(size_t i=0; i < c_vec.size(); ++i) {
+        phases[i] = std::arg(c_vec[i]);
+    }
+
+    double phase_correction = 0.0;
+    for (size_t i = 1; i < phases.size(); ++i) {
+        double diff = phases[i] - phases[i - 1];
+        if (diff > M_PI) {
+            phase_correction -= 2.0 * M_PI;
+        } else if (diff < -M_PI) {
+            phase_correction += 2.0 * M_PI;
+        }
+        phases[i] += phase_correction;
+    }
+    return phases;
+}
+
+// --- FFTW Helpers (similar to complex_block.cpp but for complex-to-complex) ---
+void TrigInterpolation::fft(const std::vector<std::complex<double>>& input, std::vector<std::complex<double>>& output) {
+    size_t n = input.size();
+    fftw_complex* in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+    fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+    for(size_t i=0; i<n; ++i) { in[i][0] = input[i].real(); in[i][1] = input[i].imag(); }
+    fftw_plan p = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
+    fftw_execute(p);
+    for(size_t i=0; i<n; ++i) { output[i] = {out[i][0], out[i][1]}; }
+    fftw_destroy_plan(p);
+    fftw_free(in); fftw_free(out);
+}
+
+void TrigInterpolation::ifft(const std::vector<std::complex<double>>& input, std::vector<std::complex<double>>& output) {
+    size_t n = input.size();
+    fftw_complex* in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+    fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);
+    for(size_t i=0; i<n; ++i) { in[i][0] = input[i].real(); in[i][1] = input[i].imag(); }
+    fftw_plan p = fftw_plan_dft_1d(n, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
+    fftw_execute(p);
+    for(size_t i=0; i<n; ++i) { output[i] = {out[i][0] / (double)n, out[i][1] / (double)n}; }
+    fftw_destroy_plan(p);
+    fftw_free(in); fftw_free(out);
+}

src/trig_interpolation.h

@@ -0,0 +1,49 @@
+#ifndef TRIG_INTERPOLATION_H
+#define TRIG_INTERPOLATION_H
+
+#include <vector>
+#include <complex>
+#include <tuple>
+#include "complex_block.h"
+
+class TrigInterpolation {
+public:
+    // User select
+    enum class PhaseMode {
+        Hilbert,
+        StandardCepstral
+    };
+
+    auto process_and_generate_fir(
+        const std::vector<double>& source_L_in, const std::vector<double>& source_R_in,
+        const std::vector<double>& target_L_in, const std::vector<double>& target_R_in,
+        int granularity, int detail_level, double strength, int fir_size, double lr_link,
+        PhaseMode phase_mode
+    ) -> std::tuple<
+        std::vector<float>, std::vector<float>, // FIR Taps (L/R)
+        std::vector<double>,                    // Frequency Axis
+        std::vector<double>, std::vector<double>, std::vector<double>, // L(src,target,match)
+        std::vector<double>, std::vector<double>, std::vector<double>  // R(src,target,match)
+    >;
+
+private:
+    // Core Smoothing and FIR Generation
+    std::vector<double> idealize_curve(const std::vector<std::complex<double>>& analytic_cepstrum, int granularity, int detail_level);
+    std::vector<double> interpolate_match_curve(
+        const std::vector<double>& smooth_source,
+        const std::vector<double>& smooth_target,
+        double strength
+    );
+    std::vector<float> create_fir_from_curve(const std::vector<double>& match_curve, int fir_size, PhaseMode phase_mode);
+    
+    // Signal Processing Helpers
+    std::vector<std::complex<double>> get_idealized_cepstrum(const std::vector<double>& signal);
+    std::vector<double> unwrap_phase(const std::vector<std::complex<double>>& c_vec);
+    void fft(const std::vector<std::complex<double>>& input, std::vector<std::complex<double>>& output);
+    void ifft(const std::vector<std::complex<double>>& input, std::vector<std::complex<double>>& output);
+
+    BlockHilbert m_block_hilbert;
+};
+
+#endif // TRIG_INTERPOLATION_H
+