Exponential Swept Sine Distortion

Question

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:

Clipped time domain meas\impulse response\window:

Frequency domain data:

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)

    @property
    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

    @property
    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(
        self,
        points,
        signal_index=None,
        start_time: float = -0.05,
        stop_time: float = 0.1,
        window="hann",
        start_percent=10,
        end_percent=10,
    ) -> 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.grid()
            plt.plot(meas_points, meas)
            plt.subplot(2, 1, 2)
            plt.grid()
            plt.plot(ir_points, impulse_response)
            plt.twinx()
        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()
    fig.suptitle("unclipped")
    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()
    fig.suptitle("clipped")
    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.figure()
    plt.grid()
    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])
    plt.show()

Answer 1

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