endolith / waveform-analysis

Functions and scripts for analyzing waveforms, primarily audio. This is currently somewhat disorganized and unfinished.
MIT License
259 stars 82 forks source link

the 997hz testing sin signal #21

Open szuer2018 opened 2 years ago

szuer2018 commented 2 years ago

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

endolith commented 2 years ago

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

szuer2018 commented 2 years ago

hi, endolith. Thanks for sharing the 997hz tesing sin signals to me. I tried out the codes on your git, but got

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')
endolith commented 3 months ago

I should add dynamic range THDN() tests for these Audition files and for the 997 Hz from 32-bit to 16-bit no dither.wav mentioned in the documentation, etc.