the 997hz testing sin signal

[szuer2018]

hi, recently, I have been trying the THD and THDN, could you offer the testing wav files?

[endolith]

https://drive.google.com/drive/folders/1mvybLleH5gvskzoBYj1TlgxTHTeebf_9?usp=sharing

[szuer2018]

hi, endolith. Thanks for sharing the 997hz tesing sin signals to me. I tried out the codes on your git, but got THDN=0.002146% or -53.37dB (THDN: 0.001404% or -97.05dB when setting weight='A') when setting weight=None in THDN function. Is that correct? Here is the codes:

if name == "main":

test_wav_path = r'.\Audition 997 Hz 16 bit.flac'

analyze_channels(test_wav_path, THDN)

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 + '"...')
signal = signal/max(abs(signal))
if channels == 1:
# Monaural
function(signal, sample_rate)
elif channels == 2:
# Stereo
if 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)

def THDN(signal, fs, weight=None):
"""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.

This notch-filters by nulling out the frequency coefficients ±10% of the
fundamental

TODO: Make R vs F reference a parameter (currently is R)
TODO: Or report all of the above in a dictionary?

"""
# Get rid of DC and window the signal
signal = np.asarray(signal) + 0.0  # Float-like array
# TODO: Do this in the frequency domain, and take any skirts with it?
signal -= mean(signal)

window = general_cosine(len(signal), flattops['HFT248D'])
windowed = signal * window
del signal

# Zero pad to nearest power of two
new_len = next_fast_len(len(windowed))
windowed = concatenate((windowed, zeros(new_len - len(windowed))))

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

# Find the peak of the frequency spectrum (fundamental frequency)
f = rfft(windowed)
i = argmax(abs(f))
true_i = parabolic(log(abs(f)), i)[0]
frequency = fs * (true_i / len(windowed))

# Filter out fundamental by throwing away values ±10%
lowermin = int(true_i * 0.9)
uppermin = int(true_i * 1.1)
f[lowermin: uppermin] = 0
# TODO: Zeroing FFT bins is bad

# Transform noise back into the time domain and measure it
noise = irfft(f)
# TODO: RMS and A-weighting in frequency domain?  Parseval?

if weight is None:
    # TODO: Return a dict or list of frequency, THD+N?
    THDN = rms_flat(noise) / total_rms
    print(f'THD: {THDN * 100}% or {dB(THDN*100)}dB')
    return THDN
elif weight == 'A':
    # Apply A-weighting to residual noise (Not normally used for
    # distortion, but used to measure dynamic range with -60 dBFS signal,
    # for instance)
    noise = A_weight(noise, fs)
    # TODO: filtfilt? tail end of filter?
    # TODO: Return a dict or list of frequency, THD+N?
    THDN = rms_flat(noise) / total_rms
    print(f'THD: {THDN * 100}% or {dB(THDN)}dB')
    return THDN
else:
    raise ValueError('Weighting not understood')