Instantly share code, notes, and snippets.

Last active May 17, 2023 06:08
THD+N calculator

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.

``````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.

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 from __future__ import division import sys from scipy.signal import blackmanharris from numpy.fft import rfft, irfft from numpy import argmax, sqrt, mean, absolute, arange, log10 import numpy as np try: import soundfile as sf except ImportError: from scikits.audiolab import Sndfile def rms_flat(a): """ Return the root mean square of all the elements of *a*, flattened out. """ return sqrt(mean(absolute(a)**2)) def find_range(f, x): """ Find range between nearest local minima from peak at index x """ for i in arange(x+1, len(f)): if f[i+1] >= f[i]: uppermin = i break for i in arange(x-1, 0, -1): if f[i] <= f[i-1]: lowermin = i + 1 break return (lowermin, uppermin) def THDN(signal, sample_rate): """ Measure the THD+N for a signal and print the results Prints the estimated fundamental frequency and the measured THD+N. This is calculated from the ratio of the entire signal before and after notch-filtering. Currently this tries to find the "skirt" around the fundamental and notch out the entire thing. A fixed-width filter would probably be just as good, if not better. """ # Get rid of DC and window the signal # TODO: Do this in the frequency domain, and take any skirts with it? signal -= mean(signal) windowed = signal * blackmanharris(len(signal)) # TODO Kaiser? # Measure the total signal before filtering but after windowing total_rms = rms_flat(windowed) # Find the peak of the frequency spectrum (fundamental frequency), and # filter the signal by throwing away values between the nearest local # minima f = rfft(windowed) i = argmax(abs(f)) # Not exact print('Frequency: %f Hz' % (sample_rate * (i / len(windowed)))) lowermin, uppermin = find_range(abs(f), i) f[lowermin: uppermin] = 0 # Transform noise back into the signal domain and measure it # TODO: 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))) def load(filename): """ Load a wave file and return the signal, sample rate and number of channels. Can be any format that libsndfile supports, like .wav, .flac, etc. """ try: wave_file = sf.SoundFile(filename) signal = wave_file.read() except ImportError: wave_file = Sndfile(filename, 'r') signal = wave_file.read_frames(wave_file.nframes) channels = wave_file.channels sample_rate = wave_file.samplerate return signal, sample_rate, channels def analyze_channels(filename, function): """ Given a filename, run the given analyzer function on each channel of the file """ signal, sample_rate, channels = load(filename) print('Analyzing "' + filename + '"...') if channels == 1: # Monaural function(signal, sample_rate) elif channels == 2: # Stereo if np.array_equal(signal[:, 0], signal[:, 1]): print('-- Left and Right channels are identical --') function(signal[:, 0], sample_rate) else: print('-- Left channel --') function(signal[:, 0], sample_rate) print('-- Right channel --') function(signal[:, 1], sample_rate) else: # Multi-channel for ch_no, channel in enumerate(signal.transpose()): print('-- Channel %d --' % (ch_no + 1)) function(channel, sample_rate) files = sys.argv[1:] if files: for filename in files: try: analyze_channels(filename, THDN) except Exception as e: print('Couldn\'t analyze "' + filename + '"') print(e) print() else: sys.exit("You must provide at least one file to analyze")

### tmbouman commented Mar 14, 2021

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!

### 200502002 commented Mar 17, 2023

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?

### endolith commented Mar 17, 2023

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

### KiDo-Ruan commented May 10, 2023

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!

### endolith commented May 10, 2023

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

### KiDo-Ruan commented May 11, 2023

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.

### endolith commented May 11, 2023

@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

### KiDo-Ruan commented May 12, 2023

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.

### endolith commented May 15, 2023

@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?

### KiDo-Ruan commented May 17, 2023

It work for me. Thanks you!

### KiDo-Ruan commented May 17, 2023

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