JAES: Synchronized Swept-Sine

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