I was able to get passable results on fundamental measurements using and exponential swept sine. Now I am trying to get distortion information from the same measurement but am puzzled by the results. I expected much a much lower distortion value for the unclipped signal. Would these be considered valid results?

Unclipped time domain meas\impulse response\window: unclipped

Clipped time domain meas\impulse response\window: clipped

Frequency domain data: frequency domain

Here is the code to get these graphs:

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

class SweptSineMeas(object):
    def __init__(self, duration, sample_rate, freq_start, freq_stop):
        self.duration = duration
        self.sample_rate = sample_rate
        self.freq_start = freq_start
        self.freq_stop = freq_stop
        self.sample_points = np.arange(0, self.duration, 1 / self.sample_rate)
        self.sweep_rate = np.log(self.freq_stop / self.freq_start)

    def stimulus(self):
        log_swept_sine = np.sin(
            (2 * np.pi * self.freq_start * self.duration / self.sweep_rate)
            * (np.exp(self.sample_points * self.sweep_rate / self.duration) - 1)
        return log_swept_sine

    def inverse_filter(self):
        decay_map = np.exp(self.sample_points * self.sweep_rate / self.duration) * 10
        inverse_filter = self.stimulus[::-1] / decay_map
        return inverse_filter

    def _impulse_reponse(self, meas, inverse_filter):
        z = np.zeros((meas.size - inverse_filter.size))
        inverse_filter = np.concatenate((inverse_filter, z))
        impulse_response = sig.fftconvolve(meas, inverse_filter, mode="same")
        return impulse_response

    def _window(
        start_time: float = -0.05,
        stop_time: float = 0.1,
    ) -> np.array:
        if signal_index is None:
            signal_index = int(points / 2)
        start_skirt_points = abs(int(start_time / (1 / self.sample_rate)))
        end_skirt_points = int(stop_time / (1 / self.sample_rate))
        window_points = start_skirt_points + end_skirt_points

        start_skirt = np.zeros(signal_index - start_skirt_points)
        start_window_points = int(window_points * (start_percent / 100))
        start_window = sig.windows.get_window(window, start_window_points * 2)
        start_window = start_window[:start_window_points]

        end_skirt = np.zeros(points - signal_index - end_skirt_points)
        end_window_points = int(window_points * (end_percent / 100))
        end_window = sig.windows.get_window(window, end_window_points * 2)
        end_window = end_window[end_window_points - 1 :: -1]

        middle_window = np.ones(window_points - (start_window.size + end_window.size))
        return np.concatenate((start_skirt, start_window, middle_window, end_window, end_skirt))

    def spectrum_mag(self, meas, window_start, window_stop, plot=False):
        impulse_response = self._impulse_reponse(meas, self.inverse_filter)
        meas_points = np.arange(0, meas.size / self.sample_rate, 1 / self.sample_rate)
        ir_points = np.arange(0, impulse_response.size / self.sample_rate, 1 / self.sample_rate)
        window = self._window(impulse_response.size, start_time=window_start, stop_time=window_stop)

        if plot is True:
            plt.subplot(2, 1, 1)
            plt.plot(meas_points, meas)
            plt.subplot(2, 1, 2)
            plt.plot(ir_points, impulse_response)
        plt.plot(ir_points, window)

        windowed_meas = impulse_response * window
        mag = np.fft.rfft(windowed_meas)
        freq = np.fft.rfftfreq(windowed_meas.size, 1 / self.sample_rate)

        return freq, 20 * np.log10(np.abs(mag))

if __name__ == "__main__":
    fund_window_start = -0.05
    fund_window_stop = 0.3
    dst_window_start = -0.4
    dst_window_stop = -0.05

    ssm = SweptSineMeas(1, 48000, 10, 10000)
    stim = ssm.stimulus
    meas = stim
    fig = plt.figure()
    freq, fnd_raw = ssm.spectrum_mag(meas, fund_window_start, fund_window_stop, plot=True)
    freq, dst_raw = ssm.spectrum_mag(meas, dst_window_start, dst_window_stop)
    meas = np.clip(stim, -0.5, 0.5)
    fig = plt.figure()
    freq, fnd_clipped = ssm.spectrum_mag(meas, fund_window_start, fund_window_stop, plot=True)
    freq, dst_clipped = ssm.spectrum_mag(meas, dst_window_start, dst_window_stop)
    plt.semilogx(freq, fnd_raw, "-r")
    plt.semilogx(freq, dst_raw, "--r")
    plt.semilogx(freq, fnd_clipped, "-g")
    plt.semilogx(freq, dst_clipped, "--g")
    plt.ylim([-18, 66])
  • $\begingroup$ Another idea. Why not compare it with values returned by REW? $\endgroup$ – jojek Sep 5 '20 at 8:06

Hello NatanBackwards and welcome to DSP SE.

Your results seem valid to me but I would like to mention a couple of things here.

First of all, you have to keep in mind that due to the way "log-sine sweeps" work the part before the impulse response is comprised of the distortion components. Judging from the way you have handled the measurements I believe you already know that, but you have to keep in mind that according to this article, there can be some "leakage" of the low frequency distortion components (odd ordered) in the actual impulse response. This may cause your results to be somewhat off from the expected ones. Additionally, you have to keep in mind that the method introduces some pre-ringing (as well as post-ringing) in most situations. More information about that and possible (partial) "remedies" can be found in this article by Angelo Farina.

The main thing to mention here is that you should be very careful of your impulse response window width. A possibly "long-enough" distortion impulse response (the impulse response of the first distortion component) may leak from the "anti-causal" to the "causal" impulse response window (the terms causal and anti-causal are used after Farina who is also the inventor of the method and are not necessarily related to the beginning of the time axis). Similarly, the pre-ringing of the impulse response may leak into the distortion impulse response window, which would result in energy from the impulse response being considered as distortion.

Other than that, two short comments are:

  1. Without further knowledge of the Device-Under-Test (DUT) we can't really conclude on the "validity" of the results you provide.
  2. A visibly clean measured response does not in any situation guarantee a distortion free transfer function. The limitation of the graphing tool/device as well as the inability of the user to notice a slight difference in a graphed function (in no way I am trying to imply that you have some kind of disability, I am just referring to the subjective ability of humans to visually distinguish small differences here) are sources of possible errors in the expectation. You can try to plot a sine-wave next to the measured spectrum and by increasing the amplitude of the sine try to see at which point you will be able to distinguish the time-domain signal is not a sine-wave anymore. You may be surprised to see how much distortion you have to introduce to visually realize the change in the input function (contrary to the hearing mechanism which may be able to distinguish the distortion earlier on).

Finally, I would like to mention that your results seem quite valid in the general sense and that by changing some things in the details (such as taking better care of the window width or position, the possible fade-in/fade-out of the input signal as well as the bandwidth) may or may not provide any considerable change in the results. I would be glad to read an update from you if you happen to reach any educated conclusions about your "problem" (and it would also benefit the community).

  • 1
    $\begingroup$ Thank you for the thorough response. A few points: 1. I am not clear from your wording but there is no DUT in my data so far this is just a simulation so far where I manually try clipped and unclipped data hoping to find meaningful response and distortion data. 2. My main issue is the fact that the unclipped distortion measurement is so high. It appears like the distortion is only 30dB down from the fundimental at 1khz. Is this correct? I would expect a value greater than 90dB on purely mathematical\simulation data. $\endgroup$ – NatanBackwards Sep 4 '20 at 15:23
  • $\begingroup$ Hhhhmmm, apologies I didn't realize you didn't use the code to measure. Well, for purely theoretical calculations the numbers are definitely "not-so-valid". One thing you chould try is to start your IR window at the first peak of the "causal" IR part and reject everything else (this will most probably introduce some phase errors though since you don't measure the initial delay of the system). Additionally, you should make sure you have correct normalizations to your calculations (summing up more points for the distortion than the IR may result in higher amplitude values). $\endgroup$ – ZaellixA Sep 4 '20 at 15:32
  • $\begingroup$ I updated the code for easy window configuration if you want to play with it. Increasing the fundamental measurement window size to be identical to the harmonic measurement didn't make my data any better. Reducing the distortion measurement window changes the measurement but doesn't make it "good". $\endgroup$ – NatanBackwards Sep 4 '20 at 16:30
  • $\begingroup$ Hhhhmmm, it seems that you may be using larger window sizes for the distortion components, which may result in different magnitudes in the frequency domain. This is implementation dependent and I will have to check it. I will make sure to keep you posted but you could try to match the sizes and see if you get any "better" results. $\endgroup$ – ZaellixA Sep 4 '20 at 18:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.