Fix IndexError in THD() for low/fractional fundamental frequencies

THD() computes the harmonic count from the continuous fundamental
(num_harmonics = int((fs/2)/frequency)) but indexes each harmonic with the
rounded fundamental bin i = int(round(frequency*N/fs)). When the fractional bin
rounds up, i * num_harmonics overshoots the rfft length len(f) = N//2+1, so the
last harmonic's bin lands past Nyquist and f[i*h] is out of bounds. For a low
fundamental (e.g. THD(sig, 12000, freq=5.994): i=6, num_harmonics=1001,
i*1001=6006 > 6001) this raises IndexError.

Guard the index and break the loop once i*h reaches len(f). A harmonic at bin
>= len(f) is above the Nyquist frequency and cannot exist in the sampled
spectrum, so it must be excluded; break (not continue) is correct because h
increases monotonically. Harmonics in range are summed exactly as before, so
results for normal signals are unchanged.

Fixes #38.

Author: youdie006

tests/test_thd.py

@@ -149,6 +149,18 @@ def test_freq_parameter(self):
         explicit_thd = THD(signal, fs, freq=f)
         assert explicit_thd == pytest.approx(auto_thd)
 
+    def test_low_fundamental_no_indexerror(self):
+        # Regression test for issue #38.  num_harmonics is derived from the
+        # continuous fundamental frequency, but harmonics are indexed by the
+        # rounded fundamental bin i.  For a low/fractional fundamental that
+        # rounds up, i * num_harmonics overshot the rfft length and raised
+        # IndexError instead of returning a THD value.
+        fs = 12000  # Hz
+        signal = sine_wave(6, fs)  # 12000 samples, fundamental ~6 Hz
+        result = THD(signal, fs, freq=5.994)  # rounds to bin 6
+        assert np.isfinite(result)
+        assert result >= 0
+
 
 if __name__ == '__main__':
     pytest.main([__file__, "--capture=sys"])

waveform_analysis/thd.py

@@ -232,6 +232,13 @@ def THD(signal, fs, *, freq=None, ref='f', verbose=False):
     num_harmonics = int((fs/2)/frequency)
     harmonic_amplitudes = []
     for h in range(2, num_harmonics + 1):
+        # num_harmonics is derived from the continuous fundamental frequency,
+        # but harmonics are indexed by the rounded fundamental bin i.  When i
+        # rounds up, i * h can overshoot the rfft length, landing above the
+        # Nyquist frequency where no component can exist.  h increases
+        # monotonically, so break once we run past the spectrum (issue #38).
+        if i * h >= len(f):
+            break
         freq = frequency * h
         ampl = abs(f[i * h])
         harmonic_amplitudes.append(ampl)