Skip to content

Instantly share code, notes, and snippets.

@endolith

endolith/readme.md

Last active October 18, 2024 13:16
Download ZIP

Revisions

  1. endolith revised this gist Jun 28, 2019. 1 changed file with 18 additions and 12 deletions.
    30 changes: 18 additions & 12 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -46,7 +46,9 @@ def THDN(signal, sample_rate):
    if not better.
    """
    # Get rid of DC and window the signal
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it?

    # 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
    @@ -57,15 +59,18 @@ def THDN(signal, sample_rate):
    # minima
    f = rfft(windowed)
    i = argmax(abs(f))
    print 'Frequency: %f Hz' % (sample_rate * (i / len(windowed))) # Not exact

    # 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
    # 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))
    print("THD+N: %.4f%% or %.1f dB" % (THDN * 100, 20 * log10(THDN)))


    def load(filename):
    @@ -93,35 +98,36 @@ def analyze_channels(filename, function):
    file
    """
    signal, sample_rate, channels = load(filename)
    print 'Analyzing "' + 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 --'
    print('-- Left and Right channels are identical --')
    function(signal[:, 0], sample_rate)
    else:
    print '-- Left channel --'
    print('-- Left channel --')
    function(signal[:, 0], sample_rate)
    print '-- Right channel --'
    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)
    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 ''
    print('Couldn\'t analyze "' + filename + '"')
    print(e)
    print()
    else:
    sys.exit("You must provide at least one file to analyze")
  2. endolith revised this gist Jun 26, 2016. 1 changed file with 48 additions and 30 deletions.
    78 changes: 48 additions & 30 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -1,20 +1,26 @@
    from __future__ import division
    import sys
    from scikits.audiolab import Sndfile
    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 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
    """
    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]:
    @@ -26,17 +32,18 @@ def find_range(f, x):
    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
    """
    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,
    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
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it?
    @@ -45,11 +52,12 @@ def THDN(signal, sample_rate):
    # 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
    # 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))
    print 'Frequency: %f Hz' % (sample_rate * (i / len(windowed))) # Not exact
    print 'Frequency: %f Hz' % (sample_rate * (i / len(windowed))) # Not exact
    lowermin, uppermin = find_range(abs(f), i)
    f[lowermin: uppermin] = 0

    @@ -59,22 +67,31 @@ def THDN(signal, sample_rate):
    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.
    """
    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.
    """
    wave_file = Sndfile(filename, 'r')
    signal = wave_file.read_frames(wave_file.nframes)
    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
    """
    """
    Given a filename, run the given analyzer function on each channel of the
    file
    """
    signal, sample_rate, channels = load(filename)
    print 'Analyzing "' + filename + '"...'

    @@ -83,14 +100,14 @@ def analyze_channels(filename, function):
    function(signal, sample_rate)
    elif channels == 2:
    # Stereo
    if np.array_equal(signal[:,0],signal[:,1]):
    if np.array_equal(signal[:, 0], signal[:, 1]):
    print '-- Left and Right channels are identical --'
    function(signal[:,0], sample_rate)
    function(signal[:, 0], sample_rate)
    else:
    print '-- Left channel --'
    function(signal[:,0], sample_rate)
    function(signal[:, 0], sample_rate)
    print '-- Right channel --'
    function(signal[:,1], sample_rate)
    function(signal[:, 1], sample_rate)
    else:
    # Multi-channel
    for ch_no, channel in enumerate(signal.transpose()):
    @@ -102,8 +119,9 @@ def analyze_channels(filename, function):
    for filename in files:
    try:
    analyze_channels(filename, THDN)
    except:
    except Exception as e:
    print 'Couldn\'t analyze "' + filename + '"'
    print e
    print ''
    else:
    sys.exit("You must provide at least one file to analyze")
  3. endolith revised this gist May 30, 2014. 1 changed file with 5 additions and 3 deletions.
    8 changes: 5 additions & 3 deletions readme.md
    View file
    Original file line number Diff line number Diff line change
    @@ -26,11 +26,13 @@ at 96 kHz, showing the 16-bit quantization distortion:
    (Is this right? Theoretical SNR of a full-scale sine is 1.761+6.02⋅16 = -98.09 dB. Close, at least.)

    * THD<sub>F</sub> is the fundamental alone vs the harmonics alone
    * THD<sub>R</sub> is the total distorted signal vs the harmonics alone (this script)
    * THD<sub>R</sub> is the total distorted signal vs the harmonics alone
    * THD+N is usually measured like THD<sub>R</sub>: the entire signal (not just
    the fundamental) vs the entire signal with the fundamental notched out. (With
    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.)
    fundamental is negligible.) This script is THD+N<sub>R</sub> (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
  4. endolith revised this gist May 30, 2014. 1 changed file with 16 additions and 11 deletions.
    27 changes: 16 additions & 11 deletions readme.md
    View file
    Original file line number Diff line number Diff line change
    @@ -1,10 +1,15 @@
    Updated version here: https://github.com/endolith/waveform-analyzer
    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.
    and notching it out in the frequency domain. This is a THD<sub>R</sub>
    measurement, meaning the denominator is the total distorted signal, not a
    bandpassed fundamental.

    Depends on Audiolab and SciPy

    @@ -20,17 +25,17 @@ at 96 kHz, showing the 16-bit quantization distortion:

    (Is this right? Theoretical SNR of a full-scale sine is 1.761+6.02⋅16 = -98.09 dB. Close, at least.)

    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.)
    * THD<sub>F</sub> is the fundamental alone vs the harmonics alone
    * THD<sub>R</sub> is the total distorted signal vs the harmonics alone (this script)
    * THD+N is usually measured like THD<sub>R</sub>: the entire signal (not just
    the fundamental) vs the entire signal with the fundamental notched out. (With
    low distortion figures, the difference between the entire signal and the
    fundamental is negligible.)

    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 the way to match
    other test equipment is to just calculate the width of a certain bandwidth,
    but is that really ideal?
    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
  5. endolith revised this gist Oct 25, 2012. 1 changed file with 2 additions and 0 deletions.
    2 changes: 2 additions & 0 deletions readme.md
    View file
    Original file line number Diff line number Diff line change
    @@ -1,3 +1,5 @@
    Updated version here: https://github.com/endolith/waveform-analyzer

    Somewhat crude THD+N calculator in Python

    Measures the total harmonic distortion plus noise (THD+N) for a given input
  6. endolith revised this gist Mar 9, 2012. 1 changed file with 3 additions and 2 deletions.
    5 changes: 3 additions & 2 deletions readme.md
    View file
    Original file line number Diff line number Diff line change
    @@ -16,8 +16,7 @@ at 96 kHz, showing the 16-bit quantization distortion:
    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.)
    (Is this right? Theoretical SNR of a full-scale sine is 1.761+6.02⋅16 = -98.09 dB. Close, at least.)

    According to the never-wrong Wikipedia:

    @@ -46,11 +45,13 @@ 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
  7. endolith renamed this gist Mar 9, 2012. 1 changed file with 31 additions and 29 deletions.
    60 changes: 31 additions & 29 deletions readme.txt → readme.md
    View file
    Original file line number Diff line number Diff line change
    @@ -5,20 +5,22 @@ signal, by guessing the fundamental frequency (finding the peak in the FFT),
    and notching it out in the frequency domain.

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

    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
    @@ -29,8 +31,8 @@ the surrounding region of the peak to throw away. Probably the way to match
    other test equipment is to just calculate the width of a certain bandwidth,
    but is that really ideal?

    width = 50
    f[i-width: i+width+1] = 0
    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,
    @@ -44,35 +46,35 @@ 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
    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
    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%
    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.
  8. endolith revised this gist Mar 12, 2010. 1 changed file with 26 additions and 9 deletions.
    35 changes: 26 additions & 9 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -27,6 +27,17 @@ def find_range(f, x):
    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
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it?
    windowed = signal * blackmanharris(len(signal)) # TODO Kaiser?
    @@ -49,44 +60,50 @@ def THDN(signal, sample_rate):
    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.
    """
    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_file(filename):
    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
    THDN(signal, sample_rate)
    function(signal, sample_rate)
    elif channels == 2:
    # Stereo
    if np.array_equal(signal[:,0],signal[:,1]):
    print '-- Left and Right channels are identical --'
    THDN(signal[:,0], sample_rate)
    function(signal[:,0], sample_rate)
    else:
    print '-- Left channel --'
    THDN(signal[:,0], sample_rate)
    function(signal[:,0], sample_rate)
    print '-- Right channel --'
    THDN(signal[:,1], sample_rate)
    function(signal[:,1], sample_rate)
    else:
    # Multi-channel
    for ch_no, channel in enumerate(signal.transpose()):
    print '-- Channel %d --' % (ch_no + 1)
    THDN(channel, sample_rate)
    function(channel, sample_rate)

    files = sys.argv[1:]
    if files:
    for filename in files:
    try:
    analyze_file(filename)
    analyze_channels(filename, THDN)
    except:
    print 'Couldn\'t analyze "' + filename + '"'
    print ''
    else:
    sys.exit("You must provide at least one file to analyze")

  9. endolith revised this gist Mar 11, 2010. 1 changed file with 6 additions and 6 deletions.
    12 changes: 6 additions & 6 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -28,7 +28,7 @@ def find_range(f, x):

    def THDN(signal, sample_rate):
    # Get rid of DC and window the signal
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it?
    windowed = signal * blackmanharris(len(signal)) # TODO Kaiser?

    # Measure the total signal before filtering but after windowing
    @@ -78,15 +78,15 @@ def analyze_file(filename):
    for ch_no, channel in enumerate(signal.transpose()):
    print '-- Channel %d --' % (ch_no + 1)
    THDN(channel, sample_rate)

    print '\n'

    files = sys.argv[1:]
    if files:
    for filename in files:
    #try:
    analyze_file(filename)
    #print 'No URL found in file', filename
    try:
    analyze_file(filename)
    except:
    print 'Couldn\'t analyze "' + filename + '"'
    print ''
    else:
    sys.exit("You must provide at least one file to analyze")

  10. endolith revised this gist Mar 11, 2010. 1 changed file with 30 additions and 19 deletions.
    49 changes: 30 additions & 19 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -55,27 +55,38 @@ def load(filename):
    sample_rate = wave_file.samplerate
    return signal, sample_rate, channels

    filename = sys.argv[1]
    signal, sample_rate, channels = load(filename)
    def analyze_file(filename):
    signal, sample_rate, channels = load(filename)

    print 'Analyzing "' + filename + '"...'
    print 'Analyzing "' + filename + '"...'

    if channels == 1:
    # Monaural
    THDN(signal, sample_rate)
    elif channels == 2:
    # Stereo
    if np.array_equal(signal[:,0],signal[:,1]):
    print '--Left and Right channels are identical:--'
    THDN(signal[:,0], sample_rate)
    if channels == 1:
    # Monaural
    THDN(signal, sample_rate)
    elif channels == 2:
    # Stereo
    if np.array_equal(signal[:,0],signal[:,1]):
    print '-- Left and Right channels are identical --'
    THDN(signal[:,0], sample_rate)
    else:
    print '-- Left channel --'
    THDN(signal[:,0], sample_rate)
    print '-- Right channel --'
    THDN(signal[:,1], sample_rate)
    else:
    print '--Left channel:--'
    THDN(signal[:,0], sample_rate)
    print '--Right channel:--'
    THDN(signal[:,1], sample_rate)
    # Multi-channel
    for ch_no, channel in enumerate(signal.transpose()):
    print '-- Channel %d --' % (ch_no + 1)
    THDN(channel, sample_rate)

    print '\n'

    files = sys.argv[1:]
    if files:
    for filename in files:
    #try:
    analyze_file(filename)
    #print 'No URL found in file', filename
    else:
    # Multi-channel
    for ch_no, channel in enumerate(signal.transpose()):
    print '--Channel %d:--' % (ch_no + 1)
    THDN(channel, sample_rate)
    sys.exit("You must provide at least one file to analyze")

  11. endolith revised this gist Mar 11, 2010. 1 changed file with 27 additions and 5 deletions.
    32 changes: 27 additions & 5 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -4,6 +4,7 @@
    from scipy.signal import blackmanharris
    from numpy.fft import rfft, irfft
    from numpy import argmax, sqrt, mean, absolute, arange, log10
    import numpy as np

    def rms_flat(a):
    """Return the root mean square of all the elements of *a*, flattened out.
    @@ -42,18 +43,39 @@ def THDN(signal, sample_rate):
    f[lowermin: uppermin] = 0

    # Transform noise back into the signal domain and measure it
    # Could probably calculate the RMS directly in the frequency domain instead
    # 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):
    wave_file = Sndfile(filename, 'r')
    signal = wave_file.read_frames(wave_file.nframes)[:,1] # TODO: Handle each channel separately
    signal = wave_file.read_frames(wave_file.nframes)
    channels = wave_file.channels
    sample_rate = wave_file.samplerate
    return signal, sample_rate
    return signal, sample_rate, channels

    filename = sys.argv[1]
    signal, sample_rate = load(filename)
    THDN(signal, sample_rate)
    signal, sample_rate, channels = load(filename)

    print 'Analyzing "' + filename + '"...'

    if channels == 1:
    # Monaural
    THDN(signal, sample_rate)
    elif channels == 2:
    # Stereo
    if np.array_equal(signal[:,0],signal[:,1]):
    print '--Left and Right channels are identical:--'
    THDN(signal[:,0], sample_rate)
    else:
    print '--Left channel:--'
    THDN(signal[:,0], sample_rate)
    print '--Right channel:--'
    THDN(signal[:,1], sample_rate)
    else:
    # Multi-channel
    for ch_no, channel in enumerate(signal.transpose()):
    print '--Channel %d:--' % (ch_no + 1)
    THDN(channel, sample_rate)

  12. endolith revised this gist Mar 11, 2010. 1 changed file with 11 additions and 11 deletions.
    22 changes: 11 additions & 11 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -25,16 +25,10 @@ def find_range(f, x):
    break
    return (lowermin, uppermin)

    def load(filename):
    wave_file = Sndfile(filename, 'r')
    signal = wave_file.read_frames(wave_file.nframes)[:,1] # TODO: Handle each channel separately
    fs = wave_file.samplerate
    return signal, fs

    def THDN(signal, fs):
    def THDN(signal, sample_rate):
    # Get rid of DC and window the signal
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it
    windowed = signal * blackmanharris(len(signal))
    windowed = signal * blackmanharris(len(signal)) # TODO Kaiser?

    # Measure the total signal before filtering but after windowing
    total_rms = rms_flat(windowed)
    @@ -43,7 +37,7 @@ def THDN(signal, fs):
    # the signal by throwing away values between the nearest local minima
    f = rfft(windowed)
    i = argmax(abs(f))
    print 'Frequency: %f Hz' % (fs * (i / len(windowed))) # Not exact
    print 'Frequency: %f Hz' % (sample_rate * (i / len(windowed))) # Not exact
    lowermin, uppermin = find_range(abs(f), i)
    f[lowermin: uppermin] = 0

    @@ -52,8 +46,14 @@ def THDN(signal, fs):
    noise = irfft(f)
    THDN = rms_flat(noise) / total_rms
    print "THD+N: %.4f%% or %.1f dB" % (THDN * 100, 20 * log10(THDN))

    def load(filename):
    wave_file = Sndfile(filename, 'r')
    signal = wave_file.read_frames(wave_file.nframes)[:,1] # TODO: Handle each channel separately
    sample_rate = wave_file.samplerate
    return signal, sample_rate

    filename = sys.argv[1]
    signal, fs = load(filename)
    THDN(signal, fs)
    signal, sample_rate = load(filename)
    THDN(signal, sample_rate)

  13. endolith revised this gist Mar 11, 2010. 1 changed file with 31 additions and 24 deletions.
    55 changes: 31 additions & 24 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -25,28 +25,35 @@ def find_range(f, x):
    break
    return (lowermin, uppermin)

    def load(filename):
    wave_file = Sndfile(filename, 'r')
    signal = wave_file.read_frames(wave_file.nframes)[:,1] # TODO: Handle each channel separately
    fs = wave_file.samplerate
    return signal, fs

    def THDN(signal, fs):
    # Get rid of DC and window the signal
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it
    windowed = signal * blackmanharris(len(signal))

    # 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))
    print 'Frequency: %f Hz' % (fs * (i / len(windowed))) # Not exact
    lowermin, uppermin = find_range(abs(f), i)
    f[lowermin: uppermin] = 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))

    filename = sys.argv[1]
    wave_file = Sndfile(filename, 'r')
    signal = wave_file.read_frames(wave_file.nframes)[:,1] # TODO: Handle each channel separately
    fs = wave_file.samplerate

    # Get rid of DC and window the signal
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it
    windowed = signal * blackmanharris(len(signal))

    # 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))
    print 'Frequency: %f Hz' % (fs * (i / len(windowed))) # Not exact
    lowermin, uppermin = find_range(abs(f), i)
    f[lowermin: uppermin] = 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))
    signal, fs = load(filename)
    THDN(signal, fs)

  14. endolith revised this gist Mar 11, 2010. 2 changed files with 6 additions and 4 deletions.
    4 changes: 2 additions & 2 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -54,8 +54,8 @@ 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:
    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%
    6 changes: 4 additions & 2 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -1,6 +1,6 @@
    from __future__ import division
    import sys
    from scikits.audiolab import flacread
    from scikits.audiolab import Sndfile
    from scipy.signal import blackmanharris
    from numpy.fft import rfft, irfft
    from numpy import argmax, sqrt, mean, absolute, arange, log10
    @@ -26,7 +26,9 @@ def find_range(f, x):
    return (lowermin, uppermin)

    filename = sys.argv[1]
    signal, fs, enc = flacread(filename)
    wave_file = Sndfile(filename, 'r')
    signal = wave_file.read_frames(wave_file.nframes)[:,1] # TODO: Handle each channel separately
    fs = wave_file.samplerate

    # Get rid of DC and window the signal
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it
  15. endolith revised this gist Mar 11, 2010. 1 changed file with 3 additions and 2 deletions.
    5 changes: 3 additions & 2 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -40,7 +40,8 @@ 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.
    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
    @@ -74,4 +75,4 @@ Frequency AP THD+N Script THD+N
    19998 20.2% 56.3% (nasty intermodulation distortion)
    20044 0.22% 0.30%

    So it's mostly accurate. Mostly.
    So it's mostly accurate. Mostly.
  16. endolith revised this gist Mar 9, 2010. 1 changed file with 1 addition and 1 deletion.
    2 changes: 1 addition & 1 deletion thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -29,7 +29,7 @@ def find_range(f, x):
    signal, fs, enc = flacread(filename)

    # Get rid of DC and window the signal
    signal -= mean(signal)
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it
    windowed = signal * blackmanharris(len(signal))

    # Measure the total signal before filtering but after windowing
  17. endolith revised this gist Mar 9, 2010. 1 changed file with 2 additions and 1 deletion.
    3 changes: 2 additions & 1 deletion readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -37,7 +37,8 @@ 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.
    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.

  18. endolith revised this gist Mar 9, 2010. 1 changed file with 2 additions and 0 deletions.
    2 changes: 2 additions & 0 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -37,6 +37,8 @@ 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:
  19. endolith revised this gist Mar 2, 2010. 1 changed file with 1 addition and 1 deletion.
    2 changes: 1 addition & 1 deletion readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -50,7 +50,7 @@ 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 22 kHz filtered vs recorded
    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
  20. endolith revised this gist Mar 2, 2010. 1 changed file with 1 addition and 1 deletion.
    2 changes: 1 addition & 1 deletion readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -35,7 +35,7 @@ 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?
    be considered part of the peak or part of the noise (jitter)?

    Also it computes the FFT for the entire sample, which is a waste of time. Use short samples.

  21. endolith revised this gist Dec 6, 2009. 1 changed file with 1 addition and 1 deletion.
    2 changes: 1 addition & 1 deletion thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -39,7 +39,7 @@ def find_range(f, x):
    # the signal by throwing away values between the nearest local minima
    f = rfft(windowed)
    i = argmax(abs(f))
    print 'Frequency: %f Hz' % (fs * (i / len(windowed)))
    print 'Frequency: %f Hz' % (fs * (i / len(windowed))) # Not exact
    lowermin, uppermin = find_range(abs(f), i)
    f[lowermin: uppermin] = 0

  22. endolith revised this gist Dec 6, 2009. 1 changed file with 8 additions and 8 deletions.
    16 changes: 8 additions & 8 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -8,21 +8,21 @@ Depends on Audiolab and SciPy
    * http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/audiolab/
    * http://www.scipy.org/

    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.
    (Is this right? Theoretical SNR of a FS sine is 1.761+6.02*16 = -98.09 dB.
    Close, at least.)

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

    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 the way to match
  23. endolith revised this gist Dec 6, 2009. 2 changed files with 4 additions and 5 deletions.
    6 changes: 3 additions & 3 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -50,8 +50,8 @@ 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:
    Comparing a test device on an Audio Precision 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%
    @@ -68,7 +68,7 @@ Frequency AP THD+N Script THD+N
    10036 0.051% 0.068%
    10723 8.2% 9.3%
    13640 12.2% 16.8%
    19998 20.2% 56.3%
    19998 20.2% 56.3% (nasty intermodulation distortion)
    20044 0.22% 0.30%

    So it's mostly accurate. Mostly.
    3 changes: 1 addition & 2 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -6,8 +6,7 @@
    from numpy import argmax, sqrt, mean, absolute, arange, log10

    def rms_flat(a):
    """
    Return the root mean square of all the elements of *a*, flattened out.
    """Return the root mean square of all the elements of *a*, flattened out.
    """
    return sqrt(mean(absolute(a)**2))
  24. endolith revised this gist Dec 6, 2009. 1 changed file with 2 additions and 0 deletions.
    2 changes: 2 additions & 0 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -37,6 +37,8 @@ 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?

    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
  25. endolith revised this gist Dec 6, 2009. 2 changed files with 46 additions and 48 deletions.
    34 changes: 26 additions & 8 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -4,8 +4,11 @@ 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:
    Depends on Audiolab and SciPy
    * http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/audiolab/
    * http://www.scipy.org/

    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
    @@ -21,6 +24,19 @@ 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.

    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 the way to match
    other test equipment is to just calculate the width of a certain bandwidth,
    but is that really ideal?

    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?

    Adobe Audition with dither:
    997 Hz 8-bit -49.8
    997 Hz 16-bit -93.4
    @@ -40,15 +56,17 @@ Frequency AP THD+N Script THD+N
    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%
    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.05% 0.05%
    4440 0.048% 0.046%
    9974 7.1% 7.8%
    10036 0.05% 0.12%
    10036 0.051% 0.068%
    10723 8.2% 9.3%
    20044 0.22% 0.60%
    13640 12.2% 16.8%
    19998 20.2% 56.3%
    20044 0.22% 0.30%

    So it's mostly accurate. Mostly.
    60 changes: 20 additions & 40 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -1,11 +1,9 @@
    from __future__ import division
    import sys
    from scikits.audiolab import flacread
    from scipy.signal import blackmanharris
    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

    from numpy import argmax, sqrt, mean, absolute, arange, log10

    def rms_flat(a):
    """
    @@ -14,55 +12,37 @@ def rms_flat(a):
    """
    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)
    def find_range(f, x):
    """Find range between nearest local minima from peak at index x
    """
    # 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)

    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)

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

    # Get rid of DC
    # Get rid of DC and window the signal
    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
    # 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)

    # 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
    i = argmax(abs(f))
    print 'Frequency: %f Hz' % (fs * (i / len(windowed)))
    lowermin, uppermin = find_range(abs(f), i)
    f[lowermin: uppermin] = 0

    # Transform noise back into the signal domain and measure it
    # Could probably calculate the RMS directly in the frequency domain instead
  26. endolith revised this gist Dec 5, 2009. 1 changed file with 12 additions and 5 deletions.
    17 changes: 12 additions & 5 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -1,19 +1,25 @@
    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.
    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.)
    * 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:
    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.
    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
    @@ -26,7 +32,8 @@ 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:
    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%
  27. endolith revised this gist Dec 5, 2009. 1 changed file with 1 addition and 0 deletions.
    1 change: 1 addition & 0 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -15,6 +15,7 @@ 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
  28. endolith revised this gist Dec 5, 2009. 1 changed file with 21 additions and 0 deletions.
    21 changes: 21 additions & 0 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -15,6 +15,10 @@ 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.

    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
    @@ -23,3 +27,20 @@ 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%
  29. endolith revised this gist Dec 5, 2009. 2 changed files with 11 additions and 12 deletions.
    20 changes: 9 additions & 11 deletions readme.txt
    View file
    Original file line number Diff line number Diff line change
    @@ -7,21 +7,19 @@ 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.)

    Comparing a test signal on an Audio Precision vs recorded into my 24-bit sound card with this script:

    Claimed freq Script freq AP THD+N Script THD+N
    440.1 440.1 0.06% 0.06%
    10003.9 10003.0 0.05% 0.07%
    15000.0 14999.1 0.13%
    19987.8 19986.0 0.22% 0.21%
    10000.0 9999.1 8.16% 7.85%
    997.0 996.9 0.06% 0.08%
    100.0 100.0 0.15% 0.14%

    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.

    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:

    3 changes: 2 additions & 1 deletion thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -61,7 +61,8 @@ def parabolic(f, x):
    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.
    f[i-10: i+11] = 0
    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
  30. endolith revised this gist Dec 4, 2009. 1 changed file with 2 additions and 2 deletions.
    4 changes: 2 additions & 2 deletions thdncalculator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -58,7 +58,7 @@ def parabolic(f, x):
    # throwing away the region around it
    i = argmax(abs(f)) # absolute?
    true_i = parabolic(abs(f), i)[0]
    print 'Frequency:\t%f Hz' % (fs * true_i / len(windowed))
    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.
    f[i-10: i+11] = 0
    @@ -67,4 +67,4 @@ def parabolic(f, x):
    # Could probably calculate the RMS directly in the frequency domain instead
    noise = irfft(f)
    THDN = rms_flat(noise) / total_rms
    print "THD+N: \t%.4f%% or %.1f dB" % (THDN * 100, 20 * log10(THDN))
    print "THD+N: %.4f%% or %.1f dB" % (THDN * 100, 20 * log10(THDN))