Synchronized Swept Sine
The synchronized-swept-sine method is a nonlinear system identification method based on "nonlinear convolution" presented by Angelo Farina in AES 108th convention in Paris in 2000. The method can analyze a nonlinear system in terms of higher order impulse responses using a single swept sine signal with user-defined frequencies and duration. The synchronization of the swept sine, described in [Novak et al. (2015)], allows the phase synchronization of the higher order impulse responses. The higher order impulse responses can be next used, in several ways, to nonlinear modeling and synthesis. The most common way is the estimation of Diagonal Volterra Kernels that are equivalent to filters of Generalized Hammerstein Model (Parallel Cascaded Hammerstein Model).
[Novak et al. (2015)] Novak, A., Lotton, P., & Simon, L. (2015). Synchronized swept-sine: Theory, application, and implementation. Journal of the Audio Engineering Society, 63(10), 786-798.
Introduction
The following figure shows the basic principle of the Higher Harmonic Frequency Responses (HHFRs). The Synchronized Swept-Sine method can estimate these HHFRs within one short measurement.
Higher Harmonic Frequency Responses (HHFRs)
Before entering the details of the Synchronized Swept-Sine method, we must first understand what the Higher Harmonic Frequency Responses (HHFRs) are. One of the artifacts of a nonlinear system is a creation of higher harmonics when the system is excited with a harmonic signal (sine or cosine signal). These higher harmonics are or can be frequency and amplitude dependent.
Sine step-by-step
The following animation shows an example of the creation of higher harmonics. A harmonic signal with frequency $f_0$ (the power spectra is depicted on the left) is sent to a nonlinear system. This output signal is distorted in a nonlinear manner which is confirmed by the presence of higher harmonics (the figure in the middle). The first harmonic (the fundamental one) with frequency $f_0$ is depicted in blue, the second harmonic has the frequency $2 f_0$ and is painted in green, and the third harmonic with frequency $3 f_0$ is depicted in red.
The values of the harmonics, relative to the value of the fundamental harmonic of the input signal, are then the HHFRs (the figure on the right).
The animation shows one way of measuring these HHFRs (or in other words the dependency of the higher harmonics with frequency), that is using a stepped sine, changing the frequency of the sine and evaluating the output spectra. This procedure can be very time-consuming.
Swept-Sine
A technique based on a swept-sine signal can estimate the HHFRs much faster than a stepped sine. The following figure shows the spectrograms of swept-sine with exponentially varying frequency. The input spectrogram is on the right, the output (distorted one) is in the middle. The higher harmonics, whose frequency components are at the integer multiples of the fundamental frequency are visible in the output spectrogram. These higher harmonics can be understood as another swept-sine signal (negatively time-delayed). This property is used in the Synchronized Swept-Sine method, and the HHFRs can be estimated (figure on the right). The details are provided below.
Swept-Sine generation
The synchronized swept-sine is generated using \begin{equation} x(t) = \sin \left[ 2 \pi f_1 L \exp \left( \dfrac{t}{L} \right) \right], \end{equation} where \begin{equation} L = \dfrac{T}{\ln \left( \frac{f_2}{f_1} \right)}, \end{equation} and where $f_1$ and $f_2$ are the initial and the final frequency and $T$ is the duration of the swept-sine.
(click to show/hide the Matlab code header ...)
%% Synchronized Swept-Sine method for analysis and identification of Nonlinear Systems
%
% The example presented in this code shows the synchronized swept-sine
% method applied to a simple non-linear system. All the steps: signal
% generation, estimation of Higher Harmonic Frequency Responses, and
% estimation of the filters of the Generalized Hammerstein model (Diagonal
% Volterra Kernels) are shown.
% [Novak et al. (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.
%
% Antonin Novak (antonin.novak@univ-lemans.fr), 26/10/2018
% Laboratoire d'Acoustique de l'Université du Mans (LAUM, UMR CNRS 6613),
% 72085 Le Mans, France
%
% https://ant-novak.com
%
%%
clc; close all; clearvars;
fs = 192000; % sampling frequency
%% -- generation of the swept-sine signal
f1 = 1e3; % start frequency
f2 = 20e3; % end frequency
T = 10; % time duration
L = T/log(f2/f1);
t = (0:1/fs:T-1/fs);
x = sin(2*pi*f1*L*(exp(t/L)));
% -- end of signal generation
Measurement / Simulation
The input signal $x(t)$ is then used as an excitation signal of the nonlinear system under test and the output $y(t)$ is recorded.
This can be done either as part of a measurement or by a simulation. The following Matlab code shows a case of a simulated nonlinear system.
%% Nonlinear system
A = [1.0,-1.93,0.94; 1.0,-1.88,0.92; 1.0,-1.83,0.95];
B = [0.2,-0.38,0.18; 0.1,-0.19,0.09; 0.02,-0.041,0.02];
NL_system = @(x) filter(B(1,:),A(1,:),x) + filter(B(2,:),A(2,:),x.^2) + filter(B(3,:),A(3,:),x.^3);
y = NL_system(x);
% -- end of nonlinear system
Swept-Sine (non-linear) deconvolution
Linear time-invariant (LTI) system refresh
The behavior of a linear time-invariant (LTI) system with input signal $x(t)$ and output signal $y(t)$ is described by the well-known convolution integral
$$
y(t) = \int_{-\infty}^{\infty} x(t-\tau) h(\tau) d\tau,
$$
where $h(t)$ is the impulse response of the system.
The Frequency Response Function (FRF) $H(f)$ is defined as
$$
H(f) = \dfrac{Y(f)}{X(f)},
$$
where $X(f) = \text{FT} \Big\{ x(t) \Big\}$ and $Y(f) = \text{FT} \Big\{ y(t) \Big\} $, FT representing the Fourier Transform. The Impulse Response (IR) $h(t)$ of the system under test is next calculated using the inverse Fourier Transform of the FRF as
$$
h(t) = \text{FT}^{-1} \Big\{ H(f) \Big\}
\label{eq:TP_signal_h(t)}
$$
In system identifiction we often call this process a deconvolution since from the input $x(t)$ and output $y(t)$ signals we can obtained back the impulse response $h(t)$ of the LTI system under test.
While it is, in general, not allowed to use the linear time-invariant (LTI) theory on non-linear systems (e.g., the terms like impulse response, frequency response function, ...) it is, thanks to a signal processing magic, used for the Synchronized Swept Sine method. The method estimates the impulse response $h(t)$ (as if it was a linear system) as \begin{equation} h(t) = \text{FT}^{-1} \Big\{ \dfrac{Y(f)}{X(f)} \Big\}, \label{eq:h(t)} \end{equation} where $X(f)$ is the spectra of input signal $x(t)$ (the swept-sine signal) and $Y(f)$ is the spectra of the output signal $y(t)$ (the response of the system to the swept-sine signal). The spectra of the swept-sine signal $x(t)$ can be calculated in two ways. Either using the Fourier transform of the signal $x(t)$, or using directly the analytical expression \begin{equation} X(f) = \frac{1}{2} \sqrt{\frac{L}{f}} \exp \left\{j 2 \pi f L \left[1-\ln \left(\frac{f}{f_1} \right) \right] - j\frac{\pi}{4} \right\}, \label{eq:X(f)} \end{equation}
%% Nonlinear (de)convolution
% Nonlinear convolution in frequency domain:
Y = fft(y)./fs; % FFT of the output signal
len = length(y);
% frequency axis
f_ax = linspace(0,fs,len+1).'; f_ax(end) = [];
% analytical expression of the spectra of the synchronized swept-sine
% [Novak et al. (2015)] "Synchronized swept-sine: Theory, application,
% and implementation. JAES 63(10), pp.786-798., Eq. (42)
X = 1/2*sqrt(L./f_ax).*exp(1i*2*pi*f_ax*L.*(1 - log(f_ax/f1)) - 1i*pi/4);
% Deconvolution
H = Y./X; H(1) = 0; % avoid Inf at DC
h = ifft(H,'symmetric'); % higher order impulse responses
% -- end of the nonlinear (de)convoltuion
(click to show/hide the Matlab code for plot of the impulse response ...)
%% -- plot the Impulse Response
figure;
plot((-length(h):length(h)-1)./fs,[h; h], 'linewidth',2);
xlim([-9 2]);
xlabel('Time [s]');
title('Impulse Response');
grid on; box on;
% -- end of plot the Impulse Response
Separation of Higher Harmonic Impulse Responses (HHIRs)
After de-convolving the output $y(t)$ with the swept-sine signal $x(t)$ using Eqs. \eqref{eq:h(t)} and \eqref{eq:X(f)} the impulse response $h(t)$ (depicted in the previous figure) seems to contain several impulse responses $h_m(t)$ delayed in time by $\Delta t_m$, mathematically written as \begin{equation} h(t) = \sum_m h_m(t + \Delta t_m), \end{equation} where \begin{equation} \Delta t_n = L \ln(n), \label{eq:Delta_t_n} \end{equation}
%% -- separation of higher harmonic impulse responses
% positions of higher harmonic IRs
N = 3; % number of harmonics to study
dt = L.*log(1:N).*fs;
% rounded positions [samples]
dt_ = round(dt);
% non-integer difference
dt_rem = dt - dt_;
% circular periodisation of IR
len_IR = 2^12;
pre_IR = len_IR / 2;
h_pos = [h; h(1:len_IR)];
% frequency axis definition
axe_w = linspace(0,2*pi,len_IR+1).'; axe_w(end) = [];
% space localisation for IRs
hm = zeros(len_IR,N);
st0 = length(h_pos);
for n = 1:N
% start position of n-th IR
st = length(h) - dt_(n) - pre_IR;
% end position of n-th IR
ed = min(st + len_IR, st0);
% separation of IRs
hm(1:ed-st,n) = h_pos(st+1:ed);
% Non-integer sample delay correction
Hx = fft(hm(:,n)) .* exp(-1i*dt_rem(n)*axe_w);
hm(:,n) = ifft(Hx,'symmetric');
% last sample poition for the next IR
st0 = st - 1;
end
% -- end of the separation of higher harmonic impulse responses
Higher Harmonic Frequency Responses (HHFRs)
The so-called impulse responses $h_m(t)$ can be separated from each other. After separation and application of Fourier Transform on each of them, we get the HHFRF $H_m(f)$ \begin{equation} H_m(f) = \text{FT} \left\{ h_m(t) \right\}. \end{equation}
%% Higher Harmonic Frequency Responses (HHFRs)
Hm = fft(hm);
% -- end of Higher Harmonic Frequency Responses (HHFRs)
(click to show/hide the Matlab code for plot of the HHFRs ...)
%% -- plot the HHFRs
figure_line_colors = [0.0, 0.0, 1.0; 0.0, 0.5, 0.0; 1.0, 0.0, 0.0];
set(groot,'defaultAxesColorOrder',figure_line_colors);
figure;
semilogx(axe_w*fs/2/pi/1000,20*log10(abs(Hm)),'linewidth',2);
xlim([f1 N*f2]./1000); ylim([-60 0]);
set(gca,'FontSize',18,'FontName','Helvetica');
title('HHFRs','FontSize',16,'FontName','Courier New');
xlabel('Frequency [kHz]','FontSize',16,'FontName','Helvetica');
ylabel('Amplitude [dB]','FontSize',16,'FontName','Helvetica');
legend({'H_1(f)','H_2(f)','H_3(f)'});
set(gca,'XTick',[1 3 5 7 10 20 40 60]);
grid on; box on;
% -- end of plot the HHFRs
Generalized Hammerstein model (Diagonal Volterra Kernels)
The estimated Higher Harmonics Frequency Responses $H_m(f)$ can be, if needed, used for the calculation of the filters of the Generalized Hammerstein model (or in other words the Diagonal Volterra Kernels) $G_n(f)$.
This is realized using a simple matrix transform from $H_m(f)$ to $G_n(f)$ as
\begin{equation} H_m(f) = \sum_{n=1}^N A_{n,m} G_n(f) \end{equation} \begin{equation} A_{n,m}=\left\{ \begin{array}{l' s} \dfrac{(-1)^{2n+\frac{1-m}{2}}}{2^{n-1}} \binom{n}{\frac{n-m}{2}}, & \text{for}~n~\geq~m~\text{and}~(n+m)~\text{is even}, \\ 0,& \text{else}. \end{array} \right. \label{eq:AMatrix} \end{equation}
%% Estimation of filters Gm of Generalized Hammerstein model (Diagonal Volterra Kernels)
% Coefficients A of the transform between the Hm and Gm
A(N,N) = 0;
for n = 1:N
for m = 1:N
if ( (n>=m) && (mod(n+m,2)==0) )
% Eq. (48) of [Novak et al. (2010)]
A(n,m) = (((-1)^(2*n+(1-m)/2))/(2^(n-1)))*nchoosek(n,(n-m)/2) ;
end
end
end
% tranform the HHFR Hm to filters Gm
Gm = Hm/A;
% -- of the Estimation of filters Gm of Generalized Hammerstein model
(click to show/hide the Matlab code for plot of the $G_m(f)$ ...)
%% -- plot the filters G_m of the Generalized Hammerstein model
% calculate the theoretical FRF of the filters to compare them
G_theo = zeros(size(Hm));
for n=1:N
G_theo(:,n) = freqz(b_coeff(n,:), a_coeff(n,:), size(Gm,1), 'whole');
end
% and compare with the estimated filters Gm with the theoretical ones
f_ax = linspace(0,fs,length(Gm)+1); f_ax(end) = [];
YLimits = ([-25 0; -25 0; -35 -5]);
for n=1:N
figure;
axes;
semilogx(f_ax/1000,20*log10(abs(Gm(:,n))),'Color',figure_line_colors(n,:),'linewidth',2); hold all;
semilogx(f_ax/1000,20*log10(abs(G_theo(:,n))),'k--','linewidth',2); hold off;
xlim([N*f1 n*f2]./1000); ylim(YLimits(n,:))
set(gca,'FontSize',18,'FontName','Helvetica');
title(['Filter G_' num2str(n) '(f)'],'FontSize',16);
xlabel('Frequency [kHz]','FontSize',16,'FontName','Helvetica');
ylabel('Amplitude [dB]','FontSize',16,'FontName','Helvetica');
legend({'estimated','theoretical'});
set(gca,'XTick',[3 5 7 10 20 40 60]);
grid on; box on;
set(gca,'OuterPosition',[0 0 1/2 1])
axes;
semilogx(f_ax/1000,angle(Gm(:,n).*exp(-1i*axe_w*pre_IR)),'Color',figure_line_colors(n,:),'linewidth',2); hold all;
semilogx(f_ax/1000,angle(G_theo(:,n)),'k--','linewidth',2); hold off;
xlim([N*f1 n*f2]./1000);
set(gca,'FontSize',18,'FontName','Helvetica');
title(['Filter G_' num2str(n) '(f)'],'FontSize',16);
xlabel('Frequency [kHz]','FontSize',16,'FontName','Helvetica');
ylabel('Phase [rad]','FontSize',16,'FontName','Helvetica');
legend({'estimated','theoretical'});
set(gca,'XTick',[3 5 7 10 20 40 60]);
grid on; box on;
set(gca,'OuterPosition',[1/2 0 1/2 1])
end
% -- end of plot the filters G_m of the Generalized Hammerstein model
The synchronized swept-sine is generated using $$ x(t) = \sin \left[ 2 \pi f_1 L \exp \left( \dfrac{t}{L} \right) \right], $$ where $$ L = \dfrac{T}{\ln \left( \frac{f_2}{f_1} \right)}, $$ and where $f_1$ and $f_2$ are the initial and the final frequency and $T$ is the duration of the swept-sine.
The input signal $x(t)$ is then used as an excitation signal of the nonlinear system under test. The output $y(t)$ is acquired and the deconvolution process, to obtain the so-called impulse response $h(t)$ is next performed in the frequency domain as $$ h(t) = \text{FT}^{-1} \Big\{ \dfrac{\text{FT}[y(t)]}{X(f)} \Big\}, $$ where $X(f)$ is defined as $$ X(f) = \frac{1}{2} \sqrt{\frac{L}{f}} \exp \left\{j 2 \pi f L \left[1-\ln \left(\frac{f}{f_1} \right) \right] - j\frac{\pi}{4} \right\}. $$ The impulse response $h(t)$ consist of time-delayed higher harmonic impulse responses $h_m(t)$ $$ h(t) = \sum_m h_m(t+\Delta_m),$$ separated by time delays $$ \Delta t_m = L \ln(m).$$
The separated impulse responses $h_m(t)$ (each impulse response corresponds to $m$-th harmonic) are then Fourier Transformed $$ H_m(\omega) = \text{FT}[ h_m(t) ]$$ to obtain the Higher Harmonic Frequency Responses (HHFRs) $H_m(f)$. These responses can be used further for the estimation of filters of the Generalized Hammrestein System (also known as Diagonal Volterra Kernels) $G_n(f)$. $$ H_m(f) = \sum_{n=1}^N A_{n,m} G_n(f) $$ $$ A_{n,m}=\left\{ \begin{array}{l' s} \dfrac{(-1)^{2n+\frac{1-m}{2}}}{2^{n-1}} \binom{n}{\frac{n-m}{2}}, & \text{for}~n~\geq~m~\text{and}~(n+m)~\text{is even}, \\ 0,& \text{else}. \end{array} \right. $$
The following figure shows an example of the estimation of a nonlinear filter (4th order Generalized Hammrestein System) using the Synchronized Swept Sine. The Matlab and Python codes are available in the Tabs above.
Download the Matlab file: Matlab_Synchronized_Swept_Sine_all_in_one.m
%% Synchronized Swept-Sine method for analysis and identification of Nonlinear Systems
%
% The example presented in this code shows the synchronized swept-sine
% method applied on a simple non-linear system. All the steps: signal
% genartion, estimation of Higher Harmonic Frequency Responses, and
% estimation of the filters of Generalized Hammerstein model (Diagonal
% Volterra Kernels) are shown.
% [Novak et al. (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.
%
% Antonin Novak (antonin.novak@univ-lemans.fr), 26/10/2018
% Laboratoire d'Acoustique de l'Université du Mans (LAUM, UMR CNRS 6613),
% 72085 Le Mans, France
%
% https://ant-novak.com
%
%%
clc; close all; clearvars;
fs = 192000; % sampling frequency
%% -- generation of the swept-sine signal
f1 = 1e3; % start frequency
f2 = 20e3; % end frequency
T = 10; % time duration
L = T/log(f2/f1);
t = (0:1/fs:T-1/fs).';
x = sin(2*pi*f1*L*(exp(t/L)));
% -- end of signal generation
%% Nonlinear system
a_coeff = [1.0,-1.9,0.94; 1.0,-1.88,0.92; 1.0,-1.83,0.95];
b_coeff = [0.2,-0.38,0.18; 0.2,-0.38,0.18; 0.04,-0.082,0.04];
NL_system = @(x) filter(b_coeff(1,:),a_coeff(1,:),x) + filter(b_coeff(2,:),a_coeff(2,:),x.^2) + filter(b_coeff(3,:),a_coeff(3,:),x.^3);
y = NL_system(x);
% -- end of nonlinear system
%% Nonlinear (de)convolution
% Nonlinear convolution in frequency domain:
Y = fft(y)./fs; % FFT of the output signal
len = length(y);
% frequency axis
f_ax = linspace(0,fs,len+1).'; f_ax(end) = [];
% analytical expression of the spectra of the synchronized swept-sine
% [Novak et al. (2015)] "Synchronized swept-sine: Theory, application,
% and implementation. JAES 63(10), pp.786-798., Eq. (42)
X = 1/2*sqrt(L./f_ax).*exp(1i*2*pi*f_ax*L.*(1 - log(f_ax/f1)) - 1i*pi/4);
% Deconvolution
H = Y./X; H(1) = 0; % avoid Inf at DC
h = ifft(H,'symmetric'); % impulse response
% -- end of the nonlinear (de)convoltuion
%% -- plot the Impulse Response
figure;
plot((-length(h):length(h)-1)./fs,[h; h], 'linewidth',2);
xlim([-9 1]);
xlabel('Time [s]');
title('Impuls Response');
grid on; box on;
% -- end of plot the Impulse Response
%% -- separation of higher harmonic impulse responses
% positions of higher harmonic IRs
N = 3; % number of harmonics to study
dt = L.*log(1:N).*fs;
% rounded positions [samples]
dt_ = round(dt);
% non-integer difference
dt_rem = dt - dt_;
% circular periodisation of IR
len_IR = 2^12;
pre_IR = len_IR / 2;
h_pos = [h; h(1:len_IR)];
% frequency axis definition
axe_w = linspace(0,2*pi,len_IR+1).'; axe_w(end) = [];
% space localisation for IRs
hm = zeros(len_IR,N);
st0 = length(h_pos);
for n = 1:N
% start position of n-th IR
st = length(h) - dt_(n) - pre_IR;
% end position of n-th IR
ed = min(st + len_IR, st0);
% separation of IRs
hm(1:ed-st,n) = h_pos(st+1:ed);
% Non-integer sample delay correction
Hx = fft(hm(:,n)) .* exp(-1i*dt_rem(n)*axe_w);
hm(:,n) = ifft(Hx,'symmetric');
% last sample poition for the next IR
st0 = st - 1;
end
% -- end of the separation of higher harmonic impulse responses
%% Higher Harmonic Frequency Responses (HHFRs)
Hm = fft(hm);
% -- end of Higher Harmonic Frequency Responses (HHFRs)
%% -- plot the HHFRs
figure_line_colors = [0.0, 0.0, 1.0; 0.0, 0.5, 0.0; 1.0, 0.0, 0.0];
set(groot,'defaultAxesColorOrder',figure_line_colors);
figure;
semilogx(axe_w*fs/2/pi/1000,20*log10(abs(Hm)),'linewidth',2);
xlim([f1 N*f2]./1000); ylim([-60 0]);
set(gca,'FontSize',18,'FontName','Helvetica');
title('HHFRs','FontSize',16,'FontName','Courier New');
xlabel('Frequency [kHz]','FontSize',16,'FontName','Helvetica');
ylabel('Amplitude [dB]','FontSize',16,'FontName','Helvetica');
legend({'H_1(f)','H_2(f)','H_3(f)'});
set(gca,'XTick',[1 3 5 7 10 20 40 60]);
grid on; box on;
% -- end of plot the HHFRs
%% Estimation of filters Gm of Generalized Hammerstein model (Diagonal Volterra Kernels)
% Coefficients A of the transform between the Hm and Gm
A(N,N) = 0;
for n = 1:N
for m = 1:N
if ( (n>=m) && (mod(n+m,2)==0) )
% Eq. (48) of [Novak et al. (2010)]
A(n,m) = (((-1)^(2*n+(1-m)/2))/(2^(n-1)))*nchoosek(n,(n-m)/2) ;
end
end
end
% tranform the HHFR Hm to filters Gm
Gm = Hm/A;
% -- of the Estimation of filters Gm of Generalized Hammerstein model
%% -- plot the filters G_m of the Generalized Hammerstein model
% calculate the theoretical FRF of the filters to compare them
G_theo = zeros(size(Hm));
for n=1:N
G_theo(:,n) = freqz(b_coeff(n,:), a_coeff(n,:), size(Gm,1), 'whole');
end
% and compare with the estimated filters Gm with the theoretical ones
f_ax = linspace(0,fs,length(Gm)+1); f_ax(end) = [];
YLimits = ([-25 0; -25 0; -35 -5]);
for n=1:N
figure;
axes;
semilogx(f_ax/1000,20*log10(abs(Gm(:,n))),'Color',figure_line_colors(n,:),'linewidth',2); hold all;
semilogx(f_ax/1000,20*log10(abs(G_theo(:,n))),'k--','linewidth',2); hold off;
xlim([N*f1 n*f2]./1000); ylim(YLimits(n,:))
title(['Filter G_' num2str(n) '(f)'],'FontSize',14);
xlabel('Frequency [kHz]','FontSize',14,'FontName','Helvetica');
ylabel('Amplitude [dB]','FontSize',14,'FontName','Helvetica');
legend({'estimated','theoretical'});
set(gca,'XTick',[3 5 7 10 20 40 60]);
grid on; box on;
set(gca,'FontSize',14,'FontName','Helvetica');
set(gca,'OuterPosition',[0 0 1/2 1])
axes;
semilogx(f_ax/1000,angle(Gm(:,n).*exp(-1i*axe_w*pre_IR)),'Color',figure_line_colors(n,:),'linewidth',2); hold all;
semilogx(f_ax/1000,angle(G_theo(:,n)),'k--','linewidth',2); hold off;
xlim([N*f1 n*f2]./1000);
title(['Filter G_' num2str(n) '(f)'],'FontSize',14);
xlabel('Frequency [kHz]','FontSize',14,'FontName','Helvetica');
ylabel('Phase [rad]','FontSize',14,'FontName','Helvetica');
legend({'estimated','theoretical'});
set(gca,'XTick',[3 5 7 10 20 40 60]);
grid on; box on;
set(gca,'FontSize',14,'FontName','Helvetica');
set(gca,'OuterPosition',[1/2 0 1/2 1])
end
% -- end of plot the filters G_m of the Generalized Hammerstein model
Download all the Matlab files: Matlab_Synchronized_Swept_Sine_functions.zip
%% Synchronized Swept-Sine method for analysis and identification of Nonlinear Systems
%
% The example presented in this code shows the synchronized swept-sine
% method applied on a simple non-linear system. All the steps: signal
% genartion, estimation of Higher Harmonic Frequency Responses, and
% estimation of the filters of Generalized Hammerstein model (Diagonal
% Volterra Kernels) are provided by the external functions.
% [Novak et al. (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.
%
% Antonin Novak (antonin.novak@univ-lemans.fr), 29/10/2018
% Laboratoire d'Acoustique de l'Université du Mans (LAUM, UMR CNRS 6613),
% 72085 Le Mans, France
%
% https://ant-novak.com
%
%%
clc; close all; clearvars;
fs = 192000; % sampling frequency
%% -- generation of the swept-sine signal
f1 = 1e3; % start frequency
f2 = 20e3; % end frequency
T = 10; % time duration
[x,t,L] = synchronized_swept_sine(f1,f2,T,fs);
% x = fadeIn_fadeOut(x,0.5,0.05,fs);
%% Nonlinear system
a_coeff = [1.0,-1.9,0.94; 1.0,-1.88,0.92; 1.0,-1.83,0.95];
b_coeff = [0.2,-0.38,0.18; 0.2,-0.38,0.18; 0.04,-0.082,0.04];
NL_system = @(x) filter(b_coeff(1,:),a_coeff(1,:),x) + filter(b_coeff(2,:),a_coeff(2,:),x.^2) + filter(b_coeff(3,:),a_coeff(3,:),x.^3);
y = NL_system(x);
%% Nonlinear (de)convolution
% Nonlinear convolution in frequency domain:
Y = fft(y)./fs; % FFT of the output signal
% analytical expression of the spectra of the synchronized swept-sine
X = synchronized_swept_sine_spectra(f1,L,fs,length(y));
% Deconvolution
H = Y./X; H(1) = 0; % avoid Inf at DC
h = ifft(H,'symmetric'); % impulse response
%% -- separation of higher harmonic impulse responses
len_IR = 2^12;
pre_IR = len_IR / 2;
N = 3; % number of harmonics to study
dt = L.*log(1:N).*fs;
hm = synchronized_swept_sine_IR_separation(h,dt,len_IR,pre_IR);
% -- end of the separation of higher harmonic impulse responses
%% Higher Harmonic Frequency Responses (HHFRs)
Hm = fft(hm);
%% -- plot the HHFRs
axe_w = linspace(0,2*pi,len_IR+1).'; axe_w(end) = [];
figure_line_colors = [0.0, 0.0, 1.0; 0.0, 0.5, 0.0; 1.0, 0.0, 0.0];
set(groot,'defaultAxesColorOrder',figure_line_colors);
figure;
semilogx(axe_w*fs/2/pi/1000,20*log10(abs(Hm)),'linewidth',2);
xlim([f1 N*f2]./1000); ylim([-60 0]);
set(gca,'FontSize',18,'FontName','Helvetica');
title('HHFRs','FontSize',16,'FontName','Courier New');
xlabel('Frequency [kHz]','FontSize',16,'FontName','Helvetica');
ylabel('Amplitude [dB]','FontSize',16,'FontName','Helvetica');
legend({'H_1(f)','H_2(f)','H_3(f)'});
set(gca,'XTick',[1 3 5 7 10 20 40 60]);
grid on; box on;
% -- end of plot the HHFRs
%% Estimation of filters Gm of Generalized Hammerstein model (Diagonal Volterra Kernels)
% Coefficients A of the transform between the Hm and Gm
A = powersin(N);
% tranform the HHFR Hm to filters Gm
Gm = Hm/A;
%% -- plot the filters G_m of the Generalized Hammerstein model
% calculate the theoretical FRF of the filters to compare them
G_theo = zeros(size(Hm));
for n=1:N
G_theo(:,n) = freqz(b_coeff(n,:), a_coeff(n,:), size(Gm,1), 'whole');
end
% and compare with the estimated filters Gm with the theoretical ones
f_ax = linspace(0,fs,length(Gm)+1); f_ax(end) = [];
YLimits = ([-25 0; -25 0; -35 -5]);
for n=1:N
figure;
axes;
semilogx(f_ax/1000,20*log10(abs(Gm(:,n))),'Color',figure_line_colors(n,:),'linewidth',2); hold all;
semilogx(f_ax/1000,20*log10(abs(G_theo(:,n))),'k--','linewidth',2); hold off;
xlim([N*f1 n*f2]./1000); ylim(YLimits(n,:))
title(['Filter G_' num2str(n) '(f)'],'FontSize',14);
xlabel('Frequency [kHz]','FontSize',14,'FontName','Helvetica');
ylabel('Amplitude [dB]','FontSize',14,'FontName','Helvetica');
legend({'estimated','theoretical'});
set(gca,'XTick',[3 5 7 10 20 40 60]);
grid on; box on;
set(gca,'FontSize',14,'FontName','Helvetica');
set(gca,'OuterPosition',[0 0 1/2 1])
axes;
semilogx(f_ax/1000,angle(Gm(:,n).*exp(-1i*axe_w*pre_IR)),'Color',figure_line_colors(n,:),'linewidth',2); hold all;
semilogx(f_ax/1000,angle(G_theo(:,n)),'k--','linewidth',2); hold off;
xlim([N*f1 n*f2]./1000);
title(['Filter G_' num2str(n) '(f)'],'FontSize',14);
xlabel('Frequency [kHz]','FontSize',14,'FontName','Helvetica');
ylabel('Phase [rad]','FontSize',14,'FontName','Helvetica');
legend({'estimated','theoretical'});
set(gca,'XTick',[3 5 7 10 20 40 60]);
grid on; box on;
set(gca,'FontSize',14,'FontName','Helvetica');
set(gca,'OuterPosition',[1/2 0 1/2 1])
end