Unfortunately, there are 2 versions of this. The other is here: https://github.com/endolith/waveform-analyzer I intend to either completely combine them or completely separate them, eventually.
Somewhat crude THD+N calculator in Python
Measures the total harmonic distortion plus noise (THD+N) for a given input signal, by guessing the fundamental frequency (finding the peak in the FFT), and notching it out in the frequency domain. This is a THDR measurement, meaning the denominator is the total distorted signal, not a bandpassed fundamental.
Depends on Audiolab and SciPy
Example of usage, with 997 Hz full-scale sine wave generated by Adobe Audition at 96 kHz, showing the 16-bit quantization distortion:
> python thdcalculator.py "perfect 997 Hz no dither.flac" Frequency: 997.000000 Hz THD+N: 0.0016% or -96.1 dB
(Is this right? Theoretical SNR of a full-scale sine is 1.761+6.02⋅16 = -98.09 dB. Close, at least.)
- THDF is the fundamental alone vs the harmonics alone
- THDR is the total distorted signal vs the harmonics alone
- THD+N is usually measured like THDR: the entire signal (not just the fundamental) vs the entire signal with the fundamental notched out (including noise and other components, not just the harmonics). (With low distortion figures, the difference between the entire signal and the fundamental is negligible.) This script is THD+NR (if that's how you write that).
The primary problem with the current script is that I don't know how much of the surrounding region of the peak to throw away. Probably should be related to the mainlobe width of the windowing function, rather than what it's currently doing. To match other test equipment would just use a fixed bandwidth, though:
width = 50 f[i-width: i+width+1] = 0
Instead of a fixed width, it currently just tries to find the nearest local minima and throw away everything between them. It works for almost all cases, but on peaks with wider "skirts", it gets stuck at any notches. Should this be considered part of the peak or part of the noise (jitter)?
By comparison, Audio Precision manual states "Bandreject Response typically –3 dB at 0.725 f0 & 1.38 f0", which is about 0.93 octaves.
Also it computes the FFT for the entire sample, which is a waste of time. Use short samples.
Adobe Audition with dither:
997 Hz 8-bit -49.8 997 Hz 16-bit -93.4 997 Hz 32-bit -143.9
Art Ludwig's Sound Files (http://members.cox.net/artludwig/):
File Claimed Measured (dB) Reference 0.0% 0.0022% -93.3 Single-ended triode 5.0% 5.06% -25.9 Solid state 0.5% 0.51% -45.8
Comparing a test device on an Audio Precision System One 22 kHz filtered vs recorded into my 96 kHz 24-bit sound card and measured with this script:
Frequency AP THD+N Script THD+N 40 1.00% 1.04% 100 0.15% 0.19% 100 0.15% 0.14% 140 0.15% 0.17% 440 0.056% 0.057% 961 0.062% 0.067% 1021 0.080% 0.082% 1440 0.042% 0.041% 1483 0.15% 0.15% 4440 0.048% 0.046% 9974 7.1% 7.8% 10036 0.051% 0.068% 10723 8.2% 9.3% 13640 12.2% 16.8% 19998 20.2% 56.3% (nasty intermodulation distortion) 20044 0.22% 0.30%
So it's mostly accurate. Mostly.
Thanks for writing this. I find for simple waveforms, a flattop window is much more accurate. So I:
changed line 3 to: from scipy.signal import flattop
changed line 52 to: windowed = signal * flattop(len(signal)) #
That drastically increased my accuracy, but flattops have wider band effects so if you have a lot of different frequency content, this may not be the best window. A hanning window is another good option for audio data. Just my 2 cents. Thanks!