diff --git a/source/dsp/PitchSmoother.h b/source/dsp/PitchSmoother.h index 1c0513f..3cbf1dc 100644 --- a/source/dsp/PitchSmoother.h +++ b/source/dsp/PitchSmoother.h @@ -41,14 +41,17 @@ class PitchSmoother return std::exp2(smoothed); } - smoothed += alpha * (logFreq - smoothed); + float delta = std::abs(logFreq - smoothed); + float effectiveAlpha = delta > 0.08f ? 1.0f : alpha; + + smoothed += effectiveAlpha * (logFreq - smoothed); return std::exp2(smoothed); } private: inline void recomputeAlpha() { - float tau = smoothingAmount * 0.2f; + float tau = smoothingAmount * 0.05f; if (tau < 1e-6f) alpha = 1.0f; else diff --git a/source/dsp/YinPitchDetector.cpp b/source/dsp/YinPitchDetector.cpp index e1d9c9d..8230cda 100644 --- a/source/dsp/YinPitchDetector.cpp +++ b/source/dsp/YinPitchDetector.cpp @@ -4,40 +4,46 @@ void YinPitchDetector::prepare(double sampleRate) { - decimation = (sampleRate > 50000.0) ? 4 : 2; - decimationCounter = 0; - decimationAccum = 0.0f; - analysisSR = sampleRate / decimation; + analysisSR = sampleRate; - windowSize = 2048; - halfWindow = windowSize / 2; + halfWindow = static_cast(std::ceil(sampleRate / 70.0)); + windowSize = 2 * halfWindow; + hopSize = static_cast(std::ceil(sampleRate * 0.003)); + + fftOrder = static_cast(std::ceil(std::log2(2.0 * windowSize))); + fftSize = 1 << fftOrder; + fft = std::make_unique(fftOrder); + fftInput.resize(static_cast(fftSize * 2), 0.0f); + fftOutput.resize(static_cast(fftSize * 2), 0.0f); buffer.assign(static_cast(windowSize), 0.0f); + linearBuffer.resize(static_cast(windowSize)); diff.resize(static_cast(halfWindow)); cmndf.resize(static_cast(halfWindow)); writePos = 0; - bufferFull = false; + hopCounter = 0; + windowFilled = false; lastResult = {}; } void YinPitchDetector::feedSample(float sample) { - decimationAccum += sample; - if (++decimationCounter < decimation) - return; - - float decimatedSample = decimationAccum / static_cast(decimation); - decimationAccum = 0.0f; - decimationCounter = 0; + buffer[static_cast(writePos)] = sample; + writePos = (writePos + 1) % windowSize; + ++hopCounter; - buffer[static_cast(writePos)] = decimatedSample; - ++writePos; + if (!windowFilled) + { + if (writePos == 0) + windowFilled = true; + else + return; + } - if (writePos >= windowSize) + if (hopCounter >= hopSize) { - writePos = 0; - bufferFull = true; + hopCounter = 0; analyse(); } } @@ -49,24 +55,47 @@ PitchResult YinPitchDetector::getResult() const void YinPitchDetector::analyse() { - if (!bufferFull) - return; - auto n = static_cast(halfWindow); - // Step 1: Difference function - for (size_t tau = 0; tau < n; ++tau) + for (int i = 0; i < windowSize; ++i) + linearBuffer[static_cast(i)] = buffer[static_cast((writePos + i) % windowSize)]; + + std::fill(fftInput.begin(), fftInput.end(), 0.0f); + for (size_t i = 0; i < n; ++i) + fftInput[i] = linearBuffer[i]; + fft->performRealOnlyForwardTransform(fftInput.data()); + + std::fill(fftOutput.begin(), fftOutput.end(), 0.0f); + for (int i = 0; i < windowSize; ++i) + fftOutput[static_cast(i)] = linearBuffer[static_cast(i)]; + fft->performRealOnlyForwardTransform(fftOutput.data()); + + for (int k = 0; k < fftSize; ++k) { - float sum = 0.0f; - for (size_t j = 0; j < n; ++j) - { - float d = buffer[j] - buffer[j + tau]; - sum += d * d; - } - diff[tau] = sum; + float aRe = fftInput[static_cast(2 * k)]; + float aIm = fftInput[static_cast(2 * k + 1)]; + float bRe = fftOutput[static_cast(2 * k)]; + float bIm = fftOutput[static_cast(2 * k + 1)]; + fftInput[static_cast(2 * k)] = aRe * bRe + aIm * bIm; + fftInput[static_cast(2 * k + 1)] = aRe * bIm - aIm * bRe; + } + + fft->performRealOnlyInverseTransform(fftInput.data()); + + float powerTerm0 = 0.0f; + for (size_t j = 0; j < n; ++j) + powerTerm0 += linearBuffer[j] * linearBuffer[j]; + + float powerTermTau = powerTerm0; + + diff[0] = 0.0f; + for (size_t tau = 1; tau < n; ++tau) + { + powerTermTau += linearBuffer[n + tau - 1] * linearBuffer[n + tau - 1] + - linearBuffer[tau - 1] * linearBuffer[tau - 1]; + diff[tau] = powerTerm0 + powerTermTau - 2.0f * fftInput[tau]; } - // Step 2: Cumulative mean normalized difference function cmndf[0] = 1.0f; float runningSum = 0.0f; @@ -79,7 +108,6 @@ void YinPitchDetector::analyse() cmndf[tau] = 1.0f; } - // Step 3: Absolute threshold size_t tauEstimate = 0; for (size_t tau = 2; tau < n; ++tau) { @@ -98,7 +126,6 @@ void YinPitchDetector::analyse() return; } - // Step 4: Parabolic interpolation float betterTau = static_cast(tauEstimate); if (tauEstimate > 0 && tauEstimate < n - 1) diff --git a/source/dsp/YinPitchDetector.h b/source/dsp/YinPitchDetector.h index 1cb869e..a327fa7 100644 --- a/source/dsp/YinPitchDetector.h +++ b/source/dsp/YinPitchDetector.h @@ -1,5 +1,7 @@ #pragma once +#include +#include #include struct PitchResult @@ -19,16 +21,21 @@ class YinPitchDetector void analyse(); double analysisSR = 44100.0; - int windowSize = 2048; - int halfWindow = 1024; - - int decimation = 1; - int decimationCounter = 0; - float decimationAccum = 0.0f; + int windowSize = 0; + int halfWindow = 0; + int hopSize = 0; std::vector buffer; + std::vector linearBuffer; int writePos = 0; - bool bufferFull = false; + int hopCounter = 0; + bool windowFilled = false; + + std::unique_ptr fft; + int fftOrder = 0; + int fftSize = 0; + std::vector fftInput; + std::vector fftOutput; std::vector diff; std::vector cmndf; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 0579ecf..ce6bb6b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -8,15 +8,15 @@ FetchContent_Declare( ) FetchContent_MakeAvailable(Catch2) -add_executable(HdnRingmodTests +juce_add_console_app(HdnRingmodTests + PRODUCT_NAME "HDN Ring Modulator Tests" +) + +target_sources(HdnRingmodTests PRIVATE TestOscillator.cpp TestYinPitchDetector.cpp TestPitchSmoother.cpp TestParameters.cpp -) - -# PitchSmoother is header-only so doesn't need a source entry here -target_sources(HdnRingmodTests PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../source/dsp/Oscillator.cpp ${CMAKE_CURRENT_SOURCE_DIR}/../source/dsp/YinPitchDetector.cpp ) @@ -25,7 +25,15 @@ target_include_directories(HdnRingmodTests PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../source ) -target_link_libraries(HdnRingmodTests PRIVATE Catch2::Catch2WithMain) +target_compile_definitions(HdnRingmodTests PRIVATE + JUCE_WEB_BROWSER=0 + JUCE_USE_CURL=0 +) + +target_link_libraries(HdnRingmodTests PRIVATE + juce::juce_dsp + Catch2::Catch2WithMain +) include(CTest) include(Catch) diff --git a/tests/TestPitchSmoother.cpp b/tests/TestPitchSmoother.cpp index 133e66f..463a946 100644 --- a/tests/TestPitchSmoother.cpp +++ b/tests/TestPitchSmoother.cpp @@ -81,7 +81,7 @@ TEST_CASE("PitchSmoother: zero smoothing passes frequency through immediately") Catch::Matchers::WithinAbs(440.0, 0.1)); } -TEST_CASE("PitchSmoother: smoothing slows convergence") +TEST_CASE("PitchSmoother: smoothing slows convergence within a semitone") { PitchSmoother fast, slow; fast.prepare(44100.0); @@ -92,32 +92,82 @@ TEST_CASE("PitchSmoother: smoothing slows convergence") slow.setSmoothingAmount(0.9f); slow.setSensitivity(1.0f); - fast.process(220.0f, 1.0f); - slow.process(220.0f, 1.0f); + fast.process(440.0f, 1.0f); + slow.process(440.0f, 1.0f); float fastOut = 0.0f, slowOut = 0.0f; - for (int i = 0; i < 1000; ++i) + for (int i = 0; i < 100; ++i) { - fastOut = fast.process(440.0f, 1.0f); - slowOut = slow.process(440.0f, 1.0f); + fastOut = fast.process(443.0f, 1.0f); + slowOut = slow.process(443.0f, 1.0f); } - REQUIRE(std::abs(fastOut - 440.0f) < std::abs(slowOut - 440.0f)); + REQUIRE(std::abs(fastOut - 443.0f) < std::abs(slowOut - 443.0f)); } -TEST_CASE("PitchSmoother: works in log-frequency domain (octave-uniform)") +TEST_CASE("PitchSmoother: smoothing operates in log-frequency domain") { PitchSmoother smoother; smoother.prepare(44100.0); smoother.setSmoothingAmount(0.5f); smoother.setSensitivity(1.0f); - smoother.process(220.0f, 1.0f); + smoother.process(440.0f, 1.0f); float output = 0.0f; - for (int i = 0; i < 4410; ++i) - output = smoother.process(440.0f, 1.0f); + for (int i = 0; i < 100; ++i) + output = smoother.process(445.0f, 1.0f); + + REQUIRE(output > 440.0f); + REQUIRE(output < 445.0f); +} + +TEST_CASE("PitchSmoother: instant lock on note change") +{ + PitchSmoother smoother; + smoother.prepare(44100.0); + smoother.setSmoothingAmount(1.0f); + smoother.setSensitivity(1.0f); - REQUIRE(output > 220.0f); - REQUIRE(output < 440.0f); + smoother.process(440.0f, 1.0f); + + float output = smoother.process(880.0f, 1.0f); + REQUIRE_THAT(static_cast(output), + Catch::Matchers::WithinAbs(880.0, 0.1)); +} + +TEST_CASE("PitchSmoother: holds during alternating confidence gaps") +{ + PitchSmoother smoother; + smoother.prepare(44100.0); + smoother.setSmoothingAmount(0.0f); + smoother.setSensitivity(0.5f); + + smoother.process(440.0f, 1.0f); + + for (int i = 0; i < 10; ++i) + { + float held = smoother.process(0.0f, 0.0f); + REQUIRE_THAT(static_cast(held), + Catch::Matchers::WithinAbs(440.0, 0.1)); + + float active = smoother.process(440.0f, 1.0f); + REQUIRE_THAT(static_cast(active), + Catch::Matchers::WithinAbs(440.0, 0.1)); + } +} + +TEST_CASE("PitchSmoother: stable output on sustained note") +{ + PitchSmoother smoother; + smoother.prepare(44100.0); + smoother.setSmoothingAmount(0.5f); + smoother.setSensitivity(1.0f); + + for (int i = 0; i < 1000; ++i) + smoother.process(440.0f, 1.0f); + + float output = smoother.process(440.0f, 1.0f); + REQUIRE_THAT(static_cast(output), + Catch::Matchers::WithinAbs(440.0, 0.01)); } diff --git a/tests/TestYinPitchDetector.cpp b/tests/TestYinPitchDetector.cpp index 7c88007..fcd74fb 100644 --- a/tests/TestYinPitchDetector.cpp +++ b/tests/TestYinPitchDetector.cpp @@ -5,6 +5,17 @@ static constexpr double twoPi = 6.283185307179586476925; +static int computeWindowSize(double sampleRate) +{ + int halfWindow = static_cast(std::ceil(sampleRate / 70.0)); + return 2 * halfWindow; +} + +static int computeHopSize(double sampleRate) +{ + return static_cast(std::ceil(sampleRate * 0.003)); +} + static void feedSine(YinPitchDetector& yin, double sampleRate, float freq, int numSamples) { double phase = 0.0; @@ -18,6 +29,23 @@ static void feedSine(YinPitchDetector& yin, double sampleRate, float freq, int n } } +static int feedSineUntilDetection(YinPitchDetector& yin, double sampleRate, float freq, int maxSamples) +{ + double phase = 0.0; + double inc = twoPi * static_cast(freq) / sampleRate; + + for (int i = 0; i < maxSamples; ++i) + { + float sample = static_cast(std::sin(phase)); + yin.feedSample(sample); + phase += inc; + + if (yin.getResult().frequency > 0.0f) + return i + 1; + } + return maxSamples; +} + TEST_CASE("YIN: detects 440 Hz sine at 44100 Hz") { YinPitchDetector yin; @@ -46,6 +74,19 @@ TEST_CASE("YIN: detects 220 Hz sine at 44100 Hz") REQUIRE(result.confidence > 0.5f); } +TEST_CASE("YIN: detects E2 (82.4 Hz) at 44100 Hz") +{ + YinPitchDetector yin; + yin.prepare(44100.0); + + feedSine(yin, 44100.0, 82.4f, 44100); + + auto result = yin.getResult(); + REQUIRE(result.frequency > 0.0f); + REQUIRE_THAT(static_cast(result.frequency), + Catch::Matchers::WithinRel(82.4, 0.03)); +} + TEST_CASE("YIN: detects 880 Hz sine at 48000 Hz") { YinPitchDetector yin; @@ -120,3 +161,61 @@ TEST_CASE("YIN: confidence is clamped to [0, 1]") REQUIRE(result.confidence >= 0.0f); REQUIRE(result.confidence <= 1.0f); } + +TEST_CASE("YIN: first detection within one window fill") +{ + YinPitchDetector yin; + yin.prepare(44100.0); + + int winSize = computeWindowSize(44100.0); + int hop = computeHopSize(44100.0); + int samplesNeeded = feedSineUntilDetection(yin, 44100.0, 440.0f, 44100); + + REQUIRE(samplesNeeded <= winSize + hop); + REQUIRE(yin.getResult().frequency > 0.0f); +} + +TEST_CASE("YIN: detects pitch change after silence") +{ + YinPitchDetector yin; + yin.prepare(44100.0); + + int winSize = computeWindowSize(44100.0); + int hop = computeHopSize(44100.0); + + for (int i = 0; i < winSize * 2; ++i) + yin.feedSample(0.0f); + + REQUIRE(yin.getResult().frequency == 0.0f); + + int samplesNeeded = feedSineUntilDetection(yin, 44100.0, 440.0f, 44100); + + REQUIRE(samplesNeeded <= winSize + hop); + + auto result = yin.getResult(); + REQUIRE(result.frequency > 0.0f); + REQUIRE_THAT(static_cast(result.frequency), + Catch::Matchers::WithinRel(440.0, 0.03)); +} + +TEST_CASE("YIN: tracks frequency sweep across hop intervals") +{ + YinPitchDetector yin; + yin.prepare(44100.0); + + int winSize = computeWindowSize(44100.0); + + feedSine(yin, 44100.0, 440.0f, winSize + 1000); + + auto resultBefore = yin.getResult(); + REQUIRE(resultBefore.frequency > 0.0f); + REQUIRE_THAT(static_cast(resultBefore.frequency), + Catch::Matchers::WithinRel(440.0, 0.03)); + + feedSine(yin, 44100.0, 880.0f, winSize + 1000); + + auto resultAfter = yin.getResult(); + REQUIRE(resultAfter.frequency > 0.0f); + REQUIRE_THAT(static_cast(resultAfter.frequency), + Catch::Matchers::WithinRel(880.0, 0.05)); +}