endolith/readme.md

## readme.md

      
    Raw
  

              readme.md
            
          
    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 THD_R
measurement, meaning the denominator is the total distorted signal, not a
bandpassed fundamental.
Depends on Audiolab and SciPy

http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/audiolab/
http://www.scipy.org/

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

THD_F is the fundamental alone vs the harmonics alone
THD_R is the total distorted signal vs the harmonics alone
THD+N is usually measured like THD_R: 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+N_R (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.

  
## thdncalculator.py
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")
	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")