THD+N calculator

readme.txt

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.

According to the never-wrong Wikipedia:

* THD is the fundamental alone vs the harmonics alone
* THD+N is the entire signal (not just the fundamental) vs the entire signal 
with the fundamental notched out.  (For low distortion, the difference between 
the entire signal and the fundamental is negligible.)

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 FS sine is 1.761+6.02*16 = -98.09 dB.  
Close, at least.

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 vs recorded into my 24-bit sound 
card 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.06%       0.06%
961         0.06%       0.07%
1021        0.08%       0.09%
1440        0.04%       0.05%
1483        0.15%       0.15%
4440        0.05%       0.05%
9974        7.1%        7.8%
10036       0.05%       0.12%
10723       8.2%        9.3%
20044       0.22%       0.60%
13640       12.2%       16.8%
19998       20.2%       56.3%

thdncalculator.py

from __future__ import division
import sys
from scikits.audiolab import flacread
from numpy.fft import rfft, irfft
from numpy import argmax, sqrt, mean, absolute, linspace, log10, logical_and, average, diff, correlate
from matplotlib.mlab import find
from scipy.signal import blackmanharris, fftconvolve


def rms_flat(a):
    """
    Return the root mean square of all the elements of *a*, flattened out.
    
    """
    return sqrt(mean(absolute(a)**2))

def parabolic(f, x):
    """Quadratic interpolation for estimating the true position of an 
    inter-sample maximum when nearby samples are known.
    
    f is a vector and x is an index for that vector.
    
    Returns (vx, vy), the coordinates of the vertex of a parabola that goes
    through point x and its two neighbors.
    
    Example:
    Defining a vector f with a local maximum at index 3 (= 6), find local
    maximum if points 2, 3, and 4 actually defined a parabola.
    
    In [3]: f = [2, 3, 1, 6, 4, 2, 3, 1]
    
    In [4]: parabolic(f, argmax(f))
    Out[4]: (3.2142857142857144, 6.1607142857142856)
    
    """
    # Requires real division.  Insert float() somewhere to force it?
    xv = 1/2 * (f[x-1] - f[x+1]) / (f[x-1] - 2 * f[x] + f[x+1]) + x
    yv = f[x] - 1/4 * (f[x-1] - f[x+1]) * (xv - x)
    return (xv, yv)


filename = sys.argv[1]
signal, fs, enc = flacread(filename)

# Get rid of DC
signal -= mean(signal)

# Window it
windowed = signal * blackmanharris(len(signal))

# Measure the total signal before filtering but after windowing
total_rms = rms_flat(windowed)

# Fourier transform of windowed signal
f = rfft(windowed)

# Find the peak of the frequency spectrum (fundamental frequency) and filter by
# throwing away the region around it
i = argmax(abs(f))  # absolute?
true_i = parabolic(abs(f), i)[0]
print 'Frequency: %f Hz' % (fs * true_i / len(windowed))

# TODO: What's the right number of coefficients to use?  Probably depends on sample length, frequency? windowing etc.
width = 10
f[i-width: i+width+1] = 0

# Transform noise back into the signal domain and measure it
# Could probably calculate the RMS directly in the frequency domain instead
noise = irfft(f)
THDN = rms_flat(noise) / total_rms
print "THD+N:     %.4f%% or %.1f dB" % (THDN * 100, 20 * log10(THDN))

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!

@tmbouman Yeah the newer version uses flattop HFT248D: https://github.com/endolith/waveform_analysis/blob/master/waveform_analysis/thd.py#L63

This is elaborate, thanks!
I am looking to do something similar but in a WinPE environment on a windows on arm device. Numpy and Scipy aren't available on windows on ARM64 yet. Is there a way to do this differently without them as dependency?

@ 200502002 I don't know about that platform, but you could re-implement this in another language.

I'm a newer in audio analyze and I'm trying to understand your code
In the line 51: Can you explain for me why need to do "signal -= mean(signal)", or did you have any document about your method?
Many thanks!

@KiDo-Ruan That's removing the DC offset, or bias, from the signal

I had use your code to calculated THD of some sample file and use APx500 machine calculated it and the result is the same (the AP use blackmanharris window too, probably :D). From your code, i want to develop it, if i want to plot the signal in frequency domain, Did you know how can i calculate amplitude of signal?
I had find out some example at stackoverflow but i don't know if they're doing it right.

@KiDo-Ruan What kind of plot are you trying to do? I have a stem plot for spectrum here https://gist.github.com/endolith/236567

FFT transforms signals from the time domain to the frequency domain, right?
I'm trying to plot time domain and frequency domain so that user can see the waveform before and after performing the FFT.
Thanks for your example, i had find out some example but most of them take fixed values ​​such as f, N. My problem is, how can i make it with a random wav file? like as analyze a 20Hz to 20kHz wav file.

@KiDo-Ruan

"see the waveform before and after performing the FFT"

Do you mean plot the waveform and the spectrum?

Have you seen https://gist.github.com/endolith/98acf9824dbf10a01795?

Have you tried asking ChatGPT to explain it and help you with the code?

It work for me. Thanks you!

I'm in Viet Nam, ChatGPT is not enabled in my country :|