JAES: Synchronized Swept-Sine
 October 30, 2015    Antonin Novak
research    research_nonlinear    A. Novak, P. Lotton & L. Simon (2015), “A. Novak, P. Lotton & L. Simon (2015), “Synchronized Swept-Sine: Theory, Application, and Implementation”, Journal of the Audio Engineering Society. Vol. 63(10), pp. 786-798.
Motivation
The paper entitled Synchronized SweptSine: Theory, Application, and Implementation deals with the design of an exponential swept-sine signal that is very often used to analyze nonlinear systems in terms of Higher Order Harmonic Frequency Responses (HHFRF). In the paper, the theory of exponential swept-sine measurements of nonlinear systems is reexamined. First, the synchronization procedure necessary for a proper analysis of higher harmonics is detailed leading to an improvement of the formula for the exponential swept-sine signal generation. Next, an analytical expression of spectra of the swept-sine signal is derived and used in the deconvolution of the impulse response. A Matlab code for generation of the synchronized swept-sine, deconvolution, and separation of the impulse responses is also given with discussion of some application issues.Swept-Sine Design
The exponential swept-sine signal is usually generated using a time-domain formula given by Angelo Farina as $$ \require{cancel} \cancel{x(t) = \sin \left\{ 2 \pi f_1 L \left[ \exp \left( \frac{t}{L} \right) -1 \right] \right\},} $$ $L$ being the rate of frequency increase. The synchronization procedure of the swept-sine, leading to a signal whose higher harmonics are perfectly synchronized with the original signal (see Figure below), is necessary for a correct phase estimation of HHFRF. $$ x(t) = \sin \left\{ 2 \pi f_1 L \left[ \exp \left( \frac{t}{L} \right) \right] \right\} $$ Contrary to the original definition the “-1” term is missing and no restriction on parameters $f_1$, $f_2$, or $T$ is necessary. As shown in the following Figure, the Synchronized swept-sine is synchronized with its harmonics.
Analytical Deconvolution
The exponential swept-sine signal $x(t)$ is then used as the input to the nonlinear system under test, and the output signal $y(t)$ is acquired. Then, the output is treated as being passed through a linear system. To obtain the impulse response $h(t)$, we simply deconvolve the output signal $y(t)$ with the original swept-sine signal $x(t)$. The deconvolution is usually made by inverting the signal $x(t)$ in time and applying an envelope adjustment. In this paper, an analytical formula for spectra of the signal $x(t)$ is derived, and the spectrum of the inverse filter $\tilde{X}(f)$ is then expressed as $$ \tilde{X}(f) = 2 \sqrt{\frac{f}{L}} \exp \left\{- j 2 \pi f L \left[1-\ln \left(\frac{f}{f_1} \right) \right] + j\frac{\pi}{4} \right\}, $$ that is used to obtain the deconvolved impulse response $h(t)$ as $$h(t) = \mathcal{F}^{-1} \left[ \mathcal{F}[y(t)] \tilde{X}(f)] \right].$$Synchronized Swept-Sine Signal Generation
The first step is the generation of the synchronized swept-sine signal. We first define the initial, final and sampling frequency and the time duration of the swept-sine signal.
f1 = 5; % initial frequency
f2 = 500; % final frequency
fs = 50000; % sampling frequency
T = 10; % atime duration
The choice of these parameters should be suited to the nonlinear system under test. In the example given in the paper, we analyzed the current distortion of a woofer using an acquisition system with sampling frequency $f_s$ = 50 kHz. Since the current distortion of a woofer is most important at very low frequencies we set the frequency range between $f_1 = 5$ Hz and $f_2 = 500$ Hz.
Next, we generate the signal using Eqs. (47) and (48)
L = round(f1/log(f2/f1)*T_)/f1;
T = L*log(f2/f1);
t = (0:ceil(fs*T)-1)./fs;
x = sin(2*pi*f1*L*(exp(t/L))).';
Windowing, fade-in and fade-out
For some applications it might be useful to apply a fade-in and/or fade-out window on the input signal. In the following code we propose a so-called raised cosine window, but any kind of window can be used.
% fade-in the input signal
fd1 = 4800; % number of samlpes
fade_in = (1-cos((0:fd1-1)/fd1*pi))/2;
index = 1:fd1;
x(index) = x(index).*fade_in.';
% fade-out the input signal
fd2 = 4800; % number of samlpes
fade_out = (1-cos((0:fd2-1)/fd2*pi))/2;
index = (1:fd2)-1;
x(end-index) = x(end-index).*fade_out.';
Note, that if we perform the deconvolution using (4) or (7) that applies the generated input signal data $x(t)$ to calculate the inverse filter, the fade-in fade-out windows are applied twice: once to the input signal $x(t)$ and once to the inverse filter. On the other hand, when the analytical formula for deconvolution (38) is used, the fading windows are applied only to the input signal $x(t)$.
Deconvolution
The deconvolution in frequency domain is implemented using Eq. (50) and the analytical expression of the inverse filter in frequency domain from Eq. (51).
% Output signal spectra
len = 2^ceil(log2(length(y)));
Y = fft(y,len)./fs; % FFT of y
% frequency axis
f_ax = linspace(0,fs,len+1).';
f_ax = f_ax(1:end-1);
% inverse filter in frequency domain
X_ = 2*sqrt(f_ax/L).*exp(-j*2*pi* ...
f_ax*L.*(1-log(f_ax/f1)) + j*pi/4);
% Deconvolution
H = Y.*X_;
H(1) = 0; % avoid Inf at DC
% higher harmonic impulse responses
h = ifft(H,'symmetric');
Separation of Higher Harmonic Impulse Responses
The deconvolution process results in the impulse response, that consists of time-delayed sum of higher harmonic impulse responses as depicted in the following Figure.
N = 3; % number of higher harmonics
% positions of higher harmonic IRs
dt = L.*log(1:N).*fs;
However, the number of samples dt defined as $\Delta t_n f_s$ is not necessarily an integer number leading to a non-integer sample shift to be realized. A schematic example of such a situation is depicted in below.

% Length of each IR (in samples)
len_IR = 2^12;
half_IR = round(len_IR/2);
% rounded positions [samples]
dt_ = round(dt);
% non-integer difference
dt_rem = dt - dt_;
% circular periodisation of IR
h_pos = [h; h(1:half_IR + len_IR)];
% frequency axis definition
axe_w = linspace(0,2*pi,len_IR+1).';
axe_w = axe_w(1:end-1);
% space localisation for IRs
hs = zeros(len_IR,N);
for n = 1:N
% start position of n-th IR
st = len-dt_(n)-half_IR;
% separation of IRs
hs(:,n) = h_pos(st:st+len_IR-1);
% apply windowing if needed
% Non-integer sample delay
Hx = fft(hs(:,n)).* ...
exp(-j*dt_rem(n)*axe_w);
hs(:,n) = ifft(Hx,'symmetric');
end
% Higher Harmonic Frequency Responses
Hs = fft(hs);
To smooth out the transitions between the impulse responses $h_n(t)$ and to control spectral ripples and time non-casual spreads, the impulse responses can be tapered prior to truncation using window functions such as the Hamming, the Bartlett, or the Kaizer windows.
Exponential swept-sine signals are very often used to analyze nonlinear audio systems. A reexamination of this methodology shows that a synchronization procedure is necessary for the proper analysis of higher harmonics. An analytical expression of spectra of the swept-sine signal is derived and used in the deconvolution of the impulse response. Matlab code for generation of the synchronized swept-sine, deconvolution, and separation of the impulse responses is given. This report provides a discussion of some application issues and an illustrative example of harmonic analysis of current distortion of a woofer. An analysis of the higher harmonics of the current distortion of a woofer is compared using both the synchronized and the non-synchronized swept-sine signals.
@article{novak2015synchronized, author={Novak, Antonin and Lotton, Pierrick and Simon, Laurent}, title={Synchronized Swept-Sine: Theory, Application, and Implementation}, journal={J. Audio Eng. Soc}, volume={63}, number={10}, pages={786-798}, year={2015}, publisher={Audio Engineering Society} }