diff --git a/.gitignore b/.gitignore index cc66a45..d20b6bd 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,11 @@ build_ios/ build_macos/ build_windows/ icons -*.png \ No newline at end of file +*.png +vamp-plugin/ +tests/ +ADC-2024/ +*.json +LICENCE +readme.md +vamp-plugin-sdk.cmake \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 76e67a7..c81ab1e 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/libraries/loop-tempo-estimator b/libraries/loop-tempo-estimator new file mode 160000 index 0000000..9369d20 --- /dev/null +++ b/libraries/loop-tempo-estimator @@ -0,0 +1 @@ +Subproject commit 9369d20106e8cbd719b8cb11b3db8b97ac403229 diff --git a/src/AudioEngine.cpp b/src/AudioEngine.cpp index 7e8f28d..e42addb 100644 --- a/src/AudioEngine.cpp +++ b/src/AudioEngine.cpp @@ -11,6 +11,38 @@ #include #include +// 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(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->setHPF(60.0f); // Cut rumble below 60Hz + p->setExpander(1.5f, -50.0f); + 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->setSmoothing(2); // Fast decay + p->setExpander(2.5f, -40.0f); + 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(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(*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 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 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); } \ No newline at end of file diff --git a/src/AudioEngine.h b/src/AudioEngine.h index 7814871..fa3be9c 100644 --- a/src/AudioEngine.h +++ b/src/AudioEngine.h @@ -18,7 +18,8 @@ public: struct FrameData { std::vector freqs; - std::vector db; + std::vector db; // Mixed (Primary + Transient) -> For Bar Height + std::vector 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& data); + void analysisReady(float bpm, float confidence); private slots: void onBufferReady(); diff --git a/src/MainWindow.cpp b/src/MainWindow.cpp index f2724ef..aa5894c 100644 --- a/src/MainWindow.cpp +++ b/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)); } diff --git a/src/MainWindow.h b/src/MainWindow.h index 81a1e39..7f29e2e 100644 --- a/src/MainWindow.h +++ b/src/MainWindow.h @@ -1,3 +1,5 @@ +// src/MainWindow.h + #pragma once #include #include @@ -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(); diff --git a/src/Processor.cpp b/src/Processor.cpp index f61acb4..8cfdee7 100644 --- a/src/Processor.cpp +++ b/src/Processor.cpp @@ -136,9 +136,14 @@ Processor::Spectrum Processor::getSpectrum() { float freq = i * (float)m_sampleRate / m_frameSize; - // HPF: Attenuate frequencies below cutoff - if (freq < m_hpfCutoff) { - mag *= 0.0001f; // -80dB attenuation + // HPF: 2nd Order Butterworth Curve + // Gain = 1 / sqrt(1 + (fc/f)^4) + 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)); diff --git a/src/VisualizerWidget.cpp b/src/VisualizerWidget.cpp index 07ff24f..59a8eff 100644 --- a/src/VisualizerWidget.cpp +++ b/src/VisualizerWidget.cpp @@ -91,6 +91,8 @@ void VisualizerWidget::updateData(const std::vector& 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& 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& 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 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) - // This fixes the "faded highs" issue by ensuring the loudest part of the spectrum - // (even if it's just high hats) is fully opaque. + // 3. Global Normalization for (float& v : vertexEnergy) { v = std::clamp(v / globalMax, 0.0f, 1.0f); } - // 4. Calculate Segment Modifiers (Vote System) - std::vector segmentModifiers(freqs.size() - 1, 0.0f); + // 4. Calculate Segment Modifiers (Procedural Pattern with Competition) + std::vector brightMods(freqs.size() - 1, 0.0f); + std::vector 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); + + // COMPETITION LOGIC (Controlled by Entropy) + // m_entropyStrength (0.0 to 3.0) + // Low Entropy = Smoother, Less Aggressive + // High Entropy = Sharper, More Faceted + + float entropyFactor = std::max(0.1f, m_entropyStrength); + + // Aggressive Contrast: Boost low sharpness + float peakIntensity = std::clamp(std::pow(sharpness * 10.0f * entropyFactor, 0.3f), 0.0f, 1.0f); - float brightBoost = 1.0f * intensity; + // Dynamic Decay: + float decayBase = 0.65f - std::clamp(sharpness * 3.0f * entropyFactor, 0.0f, 0.35f); - if (leftDominant) { - // Left Segment (i-1) gets Brightness - if (i > 0) segmentModifiers[i-1] += brightBoost; + auto applyPattern = [&](int dist, bool isBrightSide, int direction) { + int segIdx = (direction == -1) ? (i - dist) : (i + dist - 1); + if (segIdx < 0 || segIdx >= (int)brightMods.size()) return; + + int cycle = (dist - 1) / 3; + int step = (dist - 1) % 3; - // Right Side (Dark Side) - if (i < segmentModifiers.size()) segmentModifiers[i] -= 0.6f * intensity; // Dark - 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 { - // Right Segment (i) gets Brightness - if (i < segmentModifiers.size()) segmentModifiers[i] += brightBoost; + float decay = std::pow(decayBase, cycle); + float intensity = peakIntensity * decay; - // Left Side (Dark Side) - if (i > 0) segmentModifiers[i-1] -= 0.6f * intensity; // Dark - if (i >= 2) segmentModifiers[i-2] -= 0.9f * intensity; // Darker - if (i >= 3) segmentModifiers[i-3] -= 0.6f * intensity; // Dark - if (i >= 4) segmentModifiers[i-4] -= 0.3f * intensity; // Fade + 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 - float modifier = segmentModifiers[i]; - float multiplier = (modifier >= 0) ? (1.0f + modifier) : (1.0f / (1.0f - modifier * 2.0f)); - - float finalBrightness = std::clamp(baseBrightness * multiplier * m_brightness, 0.0f, 1.0f); + // Apply Brightness Modifier + float bMod = brightMods[i]; + 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 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 --- - // 1. Start with the boosted, normalized energy (avgEnergy) - // 2. Apply the shadow modifier (linear scalar for alpha) - // If modifier is negative (shadow), reduce alpha. - // If positive (peak), boost alpha slightly. - float alphaScalar = (modifier >= 0) ? (1.0f + modifier * 0.2f) : (1.0f + modifier); - + // Apply Alpha Modifier + float aMod = alphaMods[i]; + float aMult = (aMod >= 0) ? (1.0f + aMod * 0.5f) : (1.0f + aMod); + if (aMult < 0.1f) aMult = 0.1f; + float intensity = avgEnergy; float alpha = 0.4f + (intensity - 0.5f) * m_contrast; - - // Apply the shadow scalar to the calculated alpha - alpha = std::clamp(alpha * alphaScalar, 0.0f, 1.0f); + alpha = std::clamp(alpha * aMult, 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); diff --git a/src/VisualizerWidget.h b/src/VisualizerWidget.h index 5720d97..21395eb 100644 --- a/src/VisualizerWidget.h +++ b/src/VisualizerWidget.h @@ -29,7 +29,8 @@ private: struct BinState { std::deque 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 diff --git a/src/complex_block.cpp b/src/complex_block.cpp new file mode 100644 index 0000000..3666cd0 --- /dev/null +++ b/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&, const std::vector&) + * +**/ + +#include "complex_block.h" +#include +#include +#include + +std::pair>, std::vector>> +BlockHilbert::hilbertTransform(const std::vector& left_signal, const std::vector& 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(n), fft_in_L, fft_out_L, FFTW_FORWARD, FFTW_ESTIMATE); + fftw_plan plan_backward_L = fftw_plan_dft_1d(static_cast(n), ifft_in_L, ifft_out_L, FFTW_BACKWARD, FFTW_ESTIMATE); + fftw_plan plan_forward_R = fftw_plan_dft_1d(static_cast(n), fft_in_R, fft_out_R, FFTW_FORWARD, FFTW_ESTIMATE); + fftw_plan plan_backward_R = fftw_plan_dft_1d(static_cast(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> analytic_L(n); + std::vector> analytic_R(n); + for (size_t i = 0; i < n; ++i) { + analytic_L[i].real(ifft_out_L[i][0] / static_cast(n)); + analytic_L[i].imag(ifft_out_L[i][1] / static_cast(n)); + analytic_R[i].real(ifft_out_R[i][0] / static_cast(n)); + analytic_R[i].imag(ifft_out_R[i][1] / static_cast(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}; +} \ No newline at end of file diff --git a/src/complex_block.h b/src/complex_block.h new file mode 100644 index 0000000..8745f56 --- /dev/null +++ b/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 +#include +#include // For std::pair + +class BlockHilbert { +public: + std::pair>, std::vector>> + hilbertTransform(const std::vector& left_signal, const std::vector& right_signal); +}; + +#endif // COMPLEX_BLOCK_H \ No newline at end of file diff --git a/src/complex_frames.cpp b/src/complex_frames.cpp new file mode 100644 index 0000000..4647b65 --- /dev/null +++ b/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 +#include + +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>> +RealtimeHilbert::process(const std::vector& left_input_block, const std::vector& 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> result_L(m_hop_size); + std::vector> 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}; +} \ No newline at end of file diff --git a/src/complex_frames.h b/src/complex_frames.h new file mode 100644 index 0000000..139535c --- /dev/null +++ b/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 +#include +#include +#include + +class RealtimeHilbert { +public: + RealtimeHilbert(); + ~RealtimeHilbert(); + + void reinit(size_t fft_size); + std::pair>, std::vector>> + process(const std::vector& left_input_block, const std::vector& right_input_block); + +private: + void cleanup(); + + size_t m_fft_size; + size_t m_hop_size; + std::vector m_history_buffer_L; + std::vector 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 \ No newline at end of file diff --git a/src/trig_interpolation.cpp b/src/trig_interpolation.cpp new file mode 100644 index 0000000..2df5839 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include // For std::function + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +// Forward declaration for the new helper function +static std::vector create_match_curve( + const std::vector& smooth_source, + const std::vector& smooth_target, + double strength +); + +// --- FIX: Add new static helper function to correctly process magnitude spectrums --- +static std::vector> get_cepstrum_from_magnitude_spectrum( + const std::vector& mag_spectrum, + BlockHilbert& hilbert_processor, + std::function>&, std::vector>&)> 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> 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> cepstrum(n); + ifft_func(log_mag_full, cepstrum); + + // 3. Perform Hilbert transform substitution on the real part to get the analytic cepstrum. + std::vector 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& source_L_in, const std::vector& source_R_in, + const std::vector& target_L_in, const std::vector& target_R_in, + int granularity, int detail_level, double strength, int fir_size, double lr_link, + PhaseMode phase_mode +) -> std::tuple< + std::vector, std::vector, // FIR Taps (L/R) + std::vector, // Frequency Axis + std::vector, std::vector, std::vector, // L(src,tgt,match) + std::vector, std::vector, std::vector // 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 match_curve_L(pre_link_match_L.size()); + std::vector 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 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 create_match_curve( + const std::vector& smooth_source, + const std::vector& 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> 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> 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 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(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(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(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 TrigInterpolation::create_fir_from_curve(const std::vector& match_curve, int fir_size, PhaseMode phase_mode) { + if (match_curve.empty() || fir_size <= 0) return {}; + + size_t n = match_curve.size(); + std::vector> 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> 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> 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> 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> impulse_response(n); + ifft(spectrum, impulse_response); + + // Window and extract FIR taps + std::vector 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(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> TrigInterpolation::get_idealized_cepstrum(const std::vector& 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> analytic_signal = analytic_pair.first; + + // 2. Take the FFT + std::vector> 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 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> cepstrum(n); + ifft(spectrum, cepstrum); + + // 6. Perform the Hilbert transform substitution + std::vector 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 TrigInterpolation::idealize_curve(const std::vector>& analytic_cepstrum, int granularity, int detail_level) { + if (analytic_cepstrum.empty()) return {}; + + size_t n = analytic_cepstrum.size(); + std::vector 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(granularity / 100.0 * 60.0); + int iterations = 1 + static_cast(detail_level / 100.0 * 4.0); + + for (int iter = 0; iter < iterations; ++iter) { + std::vector next_curve(n, 0.0); + std::vector is_point_set(n, false); + + // 1. Binning + std::vector bin_boundaries; + double log_n = log((double)n); + for (int i = 0; i <= num_bins; ++i) { + double ratio = static_cast(i) / num_bins; + size_t boundary = static_cast(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> 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(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 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(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 TrigInterpolation::unwrap_phase(const std::vector>& c_vec) { + std::vector 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>& input, std::vector>& 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>& input, std::vector>& 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 +#include +#include +#include "complex_block.h" + +class TrigInterpolation { +public: + // User select + enum class PhaseMode { + Hilbert, + StandardCepstral + }; + + auto process_and_generate_fir( + const std::vector& source_L_in, const std::vector& source_R_in, + const std::vector& target_L_in, const std::vector& target_R_in, + int granularity, int detail_level, double strength, int fir_size, double lr_link, + PhaseMode phase_mode + ) -> std::tuple< + std::vector, std::vector, // FIR Taps (L/R) + std::vector, // Frequency Axis + std::vector, std::vector, std::vector, // L(src,target,match) + std::vector, std::vector, std::vector // R(src,target,match) + >; + +private: + // Core Smoothing and FIR Generation + std::vector idealize_curve(const std::vector>& analytic_cepstrum, int granularity, int detail_level); + std::vector interpolate_match_curve( + const std::vector& smooth_source, + const std::vector& smooth_target, + double strength + ); + std::vector create_fir_from_curve(const std::vector& match_curve, int fir_size, PhaseMode phase_mode); + + // Signal Processing Helpers + std::vector> get_idealized_cepstrum(const std::vector& signal); + std::vector unwrap_phase(const std::vector>& c_vec); + void fft(const std::vector>& input, std::vector>& output); + void ifft(const std::vector>& input, std::vector>& output); + + BlockHilbert m_block_hilbert; +}; + +#endif // TRIG_INTERPOLATION_H +