Improve DSP measurement accuracy

README.md

@@ -9,7 +9,7 @@
 [![GitHub](https://img.shields.io/badge/GitHub-CryptoJones%2FVibeComposing--Analyzer-181717?logo=github&logoColor=white)](https://github.com/CryptoJones/VibeComposing-Analyzer)
 [![C99](https://img.shields.io/badge/C-C99-A8B9CC?logo=c&logoColor=white)](waveform.c)
 [![VLC 3.0.x](https://img.shields.io/badge/VLC-3.0.x-FF8800?logo=vlcmediaplayer&logoColor=white)](https://www.videolan.org/vlc/)
-[![Version](https://img.shields.io/badge/version-3.6.0-orange)](https://github.com/CryptoJones/VibeComposing-Analyzer)
+[![Version](https://img.shields.io/badge/version-3.6.1-orange)](https://github.com/CryptoJones/VibeComposing-Analyzer)
 
 > **Sister project:** [**FL-Studio-MCP-Server**](https://github.com/CryptoJones/FL-Studio-MCP-Server)
 > — an MCP server that lets Claude Code drive FL Studio and generate editable `.flp` projects. Where this

dix_pitch.h

@@ -11,10 +11,14 @@
 #include <math.h>
 
 #define DIX_A440 440.0f
+#define DIX_PITCH_MAX_HARMONICS 6
 
 static const char *const NOTE_NAMES[12] =
 { "C", "C#", "D", "D#", "E", "F", "F#", "G", "G#", "A", "A#", "B" };
 
+static const float DIX_PITCH_HARMONIC_WEIGHTS[DIX_PITCH_MAX_HARMONICS] =
+{ 1.00f, 0.65f, 0.45f, 0.33f, 0.25f, 0.20f };
+
 /* Krumhansl-Kessler key profiles (major / natural minor tonal hierarchies). */
 static const float DIX_KRUMHANSL_MAJ[12] =
 { 6.35f, 2.23f, 3.48f, 2.33f, 4.38f, 4.09f, 2.52f, 5.19f, 2.39f, 3.66f, 2.29f, 2.88f };
@@ -41,6 +45,124 @@ static inline float dix_cents_off( float hz )
     return ( midi - (float)lroundf( midi ) ) * 100.f;
 }
 
+static inline float dix_spectrum_mag_at_bin( const float *lin_mag, int half_bins, float bin )
+{
+    if( lin_mag == 0 || half_bins <= 2 )
+        return 0.f;
+    if( bin <= 1.f )
+        return lin_mag[1];
+    if( bin >= (float)( half_bins - 1 ) )
+        return lin_mag[half_bins - 1];
+
+    int i = (int)floorf( bin );
+    float frac = bin - (float)i;
+    return lin_mag[i] * ( 1.f - frac ) + lin_mag[i + 1] * frac;
+}
+
+static inline float dix_spectrum_local_peak_bin( const float *lin_mag, int half_bins,
+                                                 float center_bin )
+{
+    if( lin_mag == 0 || half_bins <= 2 )
+        return 0.f;
+    int i = (int)lroundf( center_bin );
+    if( i < 1 ) i = 1;
+    if( i > half_bins - 2 ) i = half_bins - 2;
+    float a = lin_mag[i - 1], b = lin_mag[i], c = lin_mag[i + 1];
+    float denom = a - 2.f * b + c;
+    float delta = denom != 0.f ? 0.5f * ( a - c ) / denom : 0.f;
+    if( delta < -0.5f ) delta = -0.5f;
+    if( delta > 0.5f ) delta = 0.5f;
+    return (float)i + delta;
+}
+
+/* Harmonic pitch score for an FFT-bin candidate. Single pure tones still win at
+ * their actual bin, while a missing-fundamental stack can beat one loud overtone
+ * when several upper partials line up at integer multiples. */
+static inline float dix_harmonic_pitch_score( const float *lin_mag, int half_bins,
+                                              float candidate_bin )
+{
+    if( lin_mag == 0 || half_bins <= 2 || candidate_bin < 1.f )
+        return 0.f;
+
+    float peak = 0.f;
+    for( int h = 1; h <= DIX_PITCH_MAX_HARMONICS; h++ )
+    {
+        float hb = candidate_bin * (float)h;
+        if( hb >= (float)( half_bins - 1 ) )
+            break;
+        float m = dix_spectrum_mag_at_bin( lin_mag, half_bins, hb );
+        if( m > peak ) peak = m;
+    }
+    if( peak <= 1e-12f )
+        return 0.f;
+
+    float acc = 0.f;
+    float wsum = 0.f;
+    int hits = 0;
+    int fundamental_present = 0;
+    for( int h = 1; h <= DIX_PITCH_MAX_HARMONICS; h++ )
+    {
+        float hb = candidate_bin * (float)h;
+        if( hb >= (float)( half_bins - 1 ) )
+            break;
+        float w = DIX_PITCH_HARMONIC_WEIGHTS[h - 1];
+        float m = dix_spectrum_mag_at_bin( lin_mag, half_bins, hb );
+        acc += w * m;
+        wsum += w;
+        if( m >= peak * 0.08f )
+        {
+            hits++;
+            if( h == 1 )
+                fundamental_present = 1;
+        }
+    }
+    if( wsum <= 0.f || hits <= 0 )
+        return 0.f;
+    float continuity = sqrtf( (float)hits ) * sqrtf( sqrtf( (float)hits ) );
+    float score = ( acc / sqrtf( wsum ) ) * continuity;
+    if( !fundamental_present && hits < 3 )
+        score *= 0.35f;
+    return score;
+}
+
+static inline float dix_harmonic_pitch_refine_hz( const float *lin_mag, int half_bins,
+                                                  int fft_size, float sample_rate,
+                                                  float candidate_bin )
+{
+    if( lin_mag == 0 || half_bins <= 2 || fft_size <= 0 || sample_rate <= 0.f )
+        return 0.f;
+
+    float peak = 0.f;
+    for( int h = 1; h <= DIX_PITCH_MAX_HARMONICS; h++ )
+    {
+        float hb = candidate_bin * (float)h;
+        if( hb >= (float)( half_bins - 1 ) )
+            break;
+        float pk = dix_spectrum_local_peak_bin( lin_mag, half_bins, hb );
+        float m = dix_spectrum_mag_at_bin( lin_mag, half_bins, pk );
+        if( m > peak ) peak = m;
+    }
+    if( peak <= 1e-12f )
+        return candidate_bin * sample_rate / (float)fft_size;
+
+    float hz_acc = 0.f;
+    float hz_w = 0.f;
+    for( int h = 1; h <= DIX_PITCH_MAX_HARMONICS; h++ )
+    {
+        float hb = candidate_bin * (float)h;
+        if( hb >= (float)( half_bins - 1 ) )
+            break;
+        float pk = dix_spectrum_local_peak_bin( lin_mag, half_bins, hb );
+        float mag = dix_spectrum_mag_at_bin( lin_mag, half_bins, pk );
+        if( mag < peak * 0.08f )
+            continue;
+        float w = DIX_PITCH_HARMONIC_WEIGHTS[h - 1] * mag;
+        hz_acc += ( pk * sample_rate / (float)fft_size ) / (float)h * w;
+        hz_w += w;
+    }
+    return hz_w > 1e-12f ? hz_acc / hz_w : candidate_bin * sample_rate / (float)fft_size;
+}
+
 /* Centered (Pearson) correlation of two 12-length vectors. */
 static inline float dix_corr12( const float *a, const float *b )
 {

dix_prod.h

@@ -7,6 +7,7 @@
 #include <math.h>
 
 #define DIX_PROD_BAND_COUNT 8
+#define DIX_PROD_PI 3.14159265358979323846f
 
 enum dix_prod_band
 {
@@ -27,6 +28,16 @@ typedef struct
     float hi_hz;
 } dix_prod_band_t;
 
+typedef struct
+{
+    float b0, b1, b2, a1, a2;
+} dix_prod_biquad_t;
+
+typedef struct
+{
+    float x1, x2, y1, y2;
+} dix_prod_biquad_state_t;
+
 static const dix_prod_band_t DIX_PROD_BANDS[DIX_PROD_BAND_COUNT] =
 {
     { "SUB",   20.f,    60.f },
@@ -51,30 +62,138 @@ static inline float dix_prod_rel_db( float numerator, float denominator )
     return 10.f * log10f( ( numerator + 1e-20f ) / ( denominator + 1e-20f ) );
 }
 
-/* Mean squared FFT magnitude in a frequency range. Averaging by bin count keeps
- * wide bands from looking louder only because they contain more bins. */
-static inline float dix_prod_band_power( const float *lin_mag, int half_bins,
-                                         float sample_rate, int fft_size,
-                                         float lo_hz, float hi_hz )
+static inline dix_prod_biquad_t dix_prod_lowpass2( float sample_rate,
+                                                   float cutoff_hz, float q )
 {
-    if( lin_mag == 0 || half_bins <= 1 || sample_rate <= 0.f || fft_size <= 0 )
+    dix_prod_biquad_t c = { 1.f, 0.f, 0.f, 0.f, 0.f };
+    if( sample_rate <= 0.f || cutoff_hz <= 0.f || q <= 0.f )
+        return c;
+    if( cutoff_hz > sample_rate * 0.45f )
+        cutoff_hz = sample_rate * 0.45f;
+
+    float w0 = 2.f * DIX_PROD_PI * cutoff_hz / sample_rate;
+    float cw = cosf( w0 );
+    float alpha = sinf( w0 ) / ( 2.f * q );
+    float a0 = 1.f + alpha;
+    c.b0 = ( 1.f - cw ) * 0.5f / a0;
+    c.b1 = ( 1.f - cw ) / a0;
+    c.b2 = c.b0;
+    c.a1 = ( -2.f * cw ) / a0;
+    c.a2 = ( 1.f - alpha ) / a0;
+    return c;
+}
+
+static inline float dix_prod_biquad_process( const dix_prod_biquad_t *c,
+                                             dix_prod_biquad_state_t *s,
+                                             float x )
+{
+    float y = c->b0 * x + c->b1 * s->x1 + c->b2 * s->x2
+              - c->a1 * s->y1 - c->a2 * s->y2;
+    s->x2 = s->x1; s->x1 = x;
+    s->y2 = s->y1; s->y1 = y;
+    return y;
+}
+
+/* Cheap 4x intersample peak estimate using Catmull-Rom cubic interpolation over
+ * one segment y0..y1 with neighbor samples ym1/y2. It is not a standards-grade
+ * true-peak meter, but it catches many intersample overs better than sample peak. */
+static inline float dix_prod_interp_peak4( float ym1, float y0, float y1, float y2 )
+{
+    float peak = fabsf( y0 );
+    float ay1 = fabsf( y1 );
+    if( ay1 > peak ) peak = ay1;
+    for( int i = 1; i < 4; i++ )
+    {
+        float t = (float)i * 0.25f;
+        float t2 = t * t;
+        float t3 = t2 * t;
+        float y = 0.5f * ( 2.f * y0 + ( -ym1 + y1 ) * t
+                 + ( 2.f * ym1 - 5.f * y0 + 4.f * y1 - y2 ) * t2
+                 + ( -ym1 + 3.f * y0 - 3.f * y1 + y2 ) * t3 );
+        float ay = fabsf( y );
+        if( ay > peak ) peak = ay;
+    }
+    return peak;
+}
+
+/* FFT-bin power integrated over a frequency range with fractional edge bins.
+ * Treat each bin as covering center +/- half a bin so low bands do not jump by
+ * one whole bin when a boundary lands between centers. */
+static inline float dix_prod_band_energy( const float *lin_mag, int half_bins,
+                                          float sample_rate, int fft_size,
+                                          float lo_hz, float hi_hz )
+{
+    if( lin_mag == 0 || half_bins <= 1 || sample_rate <= 0.f || fft_size <= 0 ||
+        hi_hz <= lo_hz )
         return 0.f;
 
-    int lo = (int)floorf( lo_hz * (float)fft_size / sample_rate );
-    int hi = (int)ceilf( hi_hz * (float)fft_size / sample_rate );
+    float df = sample_rate / (float)fft_size;
+    float nyq_edge = ( (float)half_bins - 0.5f ) * df;
+    if( lo_hz < 0.f ) lo_hz = 0.f;
+    if( hi_hz > nyq_edge ) hi_hz = nyq_edge;
+    if( hi_hz <= lo_hz ) return 0.f;
+
+    double acc = 0.0;
+    int lo = (int)floorf( lo_hz / df - 0.5f );
+    int hi = (int)ceilf( hi_hz / df + 0.5f );
     if( lo < 1 ) lo = 1;
     if( hi >= half_bins ) hi = half_bins - 1;
-    if( hi < lo ) return 0.f;
 
-    double acc = 0.0;
-    int count = 0;
     for( int i = lo; i <= hi; i++ )
     {
+        float bin_lo = ( (float)i - 0.5f ) * df;
+        float bin_hi = ( (float)i + 0.5f ) * df;
+        float ov_lo = bin_lo > lo_hz ? bin_lo : lo_hz;
+        float ov_hi = bin_hi < hi_hz ? bin_hi : hi_hz;
+        float weight = ( ov_hi - ov_lo ) / df;
+        if( weight <= 0.f ) continue;
         double m = (double)lin_mag[i];
-        acc += m * m;
-        count++;
+        acc += m * m * (double)weight;
     }
-    return count > 0 ? (float)( acc / (double)count ) : 0.f;
+    return (float)acc;
+}
+
+/* Mean squared FFT magnitude in a frequency range. Averaging by effective bin
+ * count keeps wide bands from looking louder only because they contain more
+ * bins; fractional edges keep narrow low-frequency ranges stable. */
+static inline float dix_prod_band_power( const float *lin_mag, int half_bins,
+                                         float sample_rate, int fft_size,
+                                         float lo_hz, float hi_hz )
+{
+    if( lin_mag == 0 || half_bins <= 1 || sample_rate <= 0.f || fft_size <= 0 ||
+        hi_hz <= lo_hz )
+        return 0.f;
+
+    float df = sample_rate / (float)fft_size;
+    float nyq_edge = ( (float)half_bins - 0.5f ) * df;
+    float clipped_lo = lo_hz < 0.f ? 0.f : lo_hz;
+    float clipped_hi = hi_hz > nyq_edge ? nyq_edge : hi_hz;
+    if( clipped_hi <= clipped_lo ) return 0.f;
+
+    float bins = ( clipped_hi - clipped_lo ) / df;
+    if( bins <= 1e-9f ) return 0.f;
+    return dix_prod_band_energy( lin_mag, half_bins, sample_rate, fft_size,
+                                 clipped_lo, clipped_hi ) / bins;
+}
+
+/* Energy per octave. This is the useful comparator for production balance
+ * because the playbook's pink-noise target is roughly equal energy per octave,
+ * not equal energy per linear-Hz bin. */
+static inline float dix_prod_band_octave_density( const float *lin_mag, int half_bins,
+                                                  float sample_rate, int fft_size,
+                                                  float lo_hz, float hi_hz )
+{
+    if( lin_mag == 0 || half_bins <= 1 || sample_rate <= 0.f || fft_size <= 0 )
+        return 0.f;
+    float df = sample_rate / (float)fft_size;
+    float nyq_edge = ( (float)half_bins - 0.5f ) * df;
+    if( lo_hz < 1.f ) lo_hz = 1.f;
+    if( hi_hz > nyq_edge ) hi_hz = nyq_edge;
+    if( hi_hz <= lo_hz ) return 0.f;
+    float octaves = log2f( hi_hz / lo_hz );
+    if( octaves <= 1e-9f ) return 0.f;
+    return dix_prod_band_energy( lin_mag, half_bins, sample_rate, fft_size,
+                                 lo_hz, hi_hz ) / octaves;
 }
 
 static inline float dix_prod_corr( double sum_lr, double sum_ll, double sum_rr )

tests/test_pitch.c

@@ -39,6 +39,26 @@ int main( void )
            "cents_off(444) ~ +15.7 (sharp)" );
     check( dix_cents_off( 437.f ) < 0.f, "cents_off(437) is flat (negative)" );
 
+    /* Harmonic pitch scoring */
+    float spec[4096] = { 0.f };
+    spec[440] = 1.f;
+    check( dix_harmonic_pitch_score( spec, 4096, 440.f ) >
+           dix_harmonic_pitch_score( spec, 4096, 220.f ),
+           "harmonic_score pure sine stays at its actual bin" );
+    check( nearf( dix_harmonic_pitch_refine_hz( spec, 4096, 8192, 8192.f, 440.f ),
+                  440.f, 0.01f ),
+           "harmonic_refine returns pure sine frequency" );
+
+    memset( spec, 0, sizeof( spec ) );
+    spec[880] = 1.0f; spec[1320] = 0.9f; spec[1760] = 0.8f;
+    spec[2200] = 0.6f; spec[2640] = 0.5f;
+    check( dix_harmonic_pitch_score( spec, 4096, 440.f ) >
+           dix_harmonic_pitch_score( spec, 4096, 880.f ),
+           "harmonic_score can recover a missing fundamental from aligned partials" );
+    check( nearf( dix_harmonic_pitch_refine_hz( spec, 4096, 8192, 8192.f, 440.f ),
+                  440.f, 0.01f ),
+           "harmonic_refine estimates missing fundamental frequency" );
+
     /* 12-vector correlation */
     float a[12], b[12];
     for( int i = 0; i < 12; i++ ) { a[i] = (float)i; b[i] = (float)i; }

tests/test_prod.c

@@ -37,6 +37,33 @@ int main( void )
     check( nearf_( mud, 1.f, 0.0001f ), "band_power isolates the requested band" );
     check( nearf_( dix_prod_rel_db( sub, mud ), 6.0206f, 0.01f ), "rel_db returns 6.02 dB for 4x power" );
 
+    dix_prod_biquad_t lp = dix_prod_lowpass2( 48000.f, 200.f, 0.70710678f );
+    float dc_gain = ( lp.b0 + lp.b1 + lp.b2 ) / ( 1.f + lp.a1 + lp.a2 );
+    check( nearf_( dc_gain, 1.f, 0.001f ), "lowpass2 has unity DC gain" );
+    dix_prod_biquad_state_t st = { 0.f, 0.f, 0.f, 0.f };
+    float y = 0.f;
+    for( int i = 0; i < 4000; i++ )
+        y = dix_prod_biquad_process( &lp, &st, 1.f );
+    check( y > 0.99f && y < 1.01f, "lowpass2 passes steady DC" );
+
+    check( nearf_( dix_prod_interp_peak4( 0.5f, 0.5f, 0.5f, 0.5f ), 0.5f, 0.0001f ),
+           "interp_peak4 keeps flat signal peak" );
+    check( dix_prod_interp_peak4( 0.f, 0.8660254f, 0.8660254f, 0.f ) > 0.90f,
+           "interp_peak4 catches between-sample peak" );
+
+    float edge[512] = { 0.f };
+    edge[10] = 2.f;
+    check( nearf_( dix_prod_band_power( edge, 512, 1024.f, 1024, 9.75f, 10.25f ), 4.f, 0.0001f ),
+           "band_power handles fractional bin edges" );
+
+    float pinkish[512] = { 0.f };
+    for( int i = 1; i < 512; i++ )
+        pinkish[i] = 1.f / sqrtf( (float)i );
+    float oct1 = dix_prod_band_octave_density( pinkish, 512, 1024.f, 1024, 20.f, 40.f );
+    float oct2 = dix_prod_band_octave_density( pinkish, 512, 1024.f, 1024, 40.f, 80.f );
+    check( fabsf( dix_prod_rel_db( oct1, oct2 ) ) < 0.25f,
+           "octave_density treats pink-like energy as balanced across octaves" );
+
     check( nearf_( dix_prod_corr( 4.0, 4.0, 4.0 ), 1.f, 0.0001f ), "corr identical channels = 1" );
     check( nearf_( dix_prod_corr( -4.0, 4.0, 4.0 ), -1.f, 0.0001f ), "corr inverted channels = -1" );
     check( nearf_( dix_prod_width( 9.0, 0.0 ), 0.f, 0.0001f ), "width mono center = 0" );