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.
In order to separate them from each other, we must first calculate the time-moments $\Delta t_n$ in which these higher harmonic impulse responses appear using Eq. (52). In order to obtain directly the number of samples by which the higher harmonic impulse responses are shifted, we multiply the time delays $\Delta t_n$ by the sampling frequency $f_s$.

	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.
Thus, we first separate the higher harmonic impulse responses at the integer shift calculated as the rounded version of the exact non-integer number of samples and, after the separation, we implement the non-integer delay in frequency domain using the time shift property of the Fourier transform $$ \mathcal{F} \left[ f(t-t_0) \right] = \exp (-j \omega t_0) \mathcal{F} \left[ f(t) \right]. $$ Finally, the HHFRs $H_n(f)$ are calculated as the Fourier transform of the higher harmonic impulse responses $h_n(t)$.

	% 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}
}