2
$\begingroup$

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 a much lower distortion value for the unclipped signal. Would these be considered valid results?

Unclipped time domain measurement - impulse response - window: unclipped

Clipped time domain measurement - 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)

    @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()
$\endgroup$
3
  • $\begingroup$ Another idea. Why not compare it with values returned by REW? $\endgroup$
    – jojeck
    Commented Sep 5, 2020 at 8:06
  • $\begingroup$ I am having the same problem. Did you find the answer? Thanks, mm $\endgroup$ Commented Aug 18, 2022 at 14:23
  • $\begingroup$ @MohammadMohammadi Welcome to SE.SP! Please do not post an "answer" that is not an answer to the question asked. If you have a related question, please ask a new question and refer to this question... and state specifically why the answer here doesn't help you (so we can help!). $\endgroup$
    – Peter K.
    Commented Aug 18, 2022 at 15:27

2 Answers 2

1
$\begingroup$

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

$\endgroup$
4
  • 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$ Commented Sep 4, 2020 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
    Commented Sep 4, 2020 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$ Commented Sep 4, 2020 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
    Commented Sep 4, 2020 at 18:05
0
$\begingroup$

One of the problems with logsweeps is that the first and last cycles have a discontinuity. If you insert zeros before and after the signal, then at the first sample it goes from zero to sin(wt), so the first derivative goes from zero to wcos(wt), which means it has a step.

The same happens at the end of the signal. This is like multiplying an infinite length signal with a rectangular window, so in the frequency domain your spectrum will be convolved with the spectrum of the rectangular window, and that has a lot of skirts which will pollute the results and make them unusable.

If you don't insert zeros before and after the signal... FFT still considers your signal as periodic, so it has a discontinuity between the last sample and the first sample. This spreads lots of garbage all over the spectrum.

One solution is to apply "fade in" and "fade out", for example by multiplying the beginning and end of the signal with your favorite window function, sliced in half so it makes a nice smooth fade.

If you intend to go up to Fs/2 then this means you'll lose a bit of frequency response at the top. It is not a problem at the bottom, since you can just use a lower starting frequency.

In addition, if your DUT has any kind of impulse response, you have to acquire the output signal after the sweep ends in order to measure it, so the resulting input signal will have the sweep in the middle surrounded by zeros.

If your DUT is AC-coupled (or you use an AC-coupled soundcard) you must be very careful about the first cycles, because the first cycle of the waveform (going up) will bias the AC coupling caps and result in a rather long settling time. This is another way to say a gated sinewave contains frequencies both above and below its oscillation frequency, so you will see the settling time of the analog highpass filter in the result. A solution is to carefully adjust the fade-in rate to minimize the amount of this "DC shift" at the beginning of the waveform.

The length of the FFT should be a power of two for speed, so you can adjust how much sweep you put into it depending on the expected length of the DUT impulse response.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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