Antonin Novak - homepage


Postdoctoral Research Fellow

Laboratoire d'Acoustique de l'Université du Maine
UMR-CNRS 6613
Université du Maine
Av. Olivier Messiaen
72085 Le Mans Cedex 9

Email Address: ant.novak[at]gmail.com

Antonin Novak





Nonlinear System Identification using Synchronized Swept-Sine Method

(The whole matlab code can be downloaded here: swept_sine_demo.m)

Published in: A. Novak, L. Simon, F. Kadlec & P. Lotton (2010), "Nonlinear system identification using exponential swept-sine signal", Instrumentation and Measurement, IEEE Transactions on. Vol. 59(8), pp. 2220-2229. (download)



The following MATLAB® code is a tool for the identification of a nonlinear system using Synchronized Swept Sine method.

The method is based on the estimation of the Higher Harmonic Frequency Responses (HHFRFs) of distorted signal "y(t)". The nonlinear system under test is then modelled by non-parametric generalized Hammerstein model made up of power series associated with linear filters as depicted in the following Figure.



MATLAB® Code

The whole matlab code can be downloaded here: swept_sine_demo.m

clc; close all ; clear all;
fs = 96000;         % sampling frequency

f1 = 10;            % start frequency
f2 = 12000;         % final frequency
T_ = 10;            % approximative time duration

N = 4;              % order of nonlinearities

% Length of each impulse repsonse (in samples)
len_IR = 2^13;
pre_IR = len_IR/2;

Generation of the signal

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)-1));

% fade-in the input signal
fd1 = 480;  % 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 =  480;  % 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;


Nonlinear System definition

% the ARMA filters definition (ARMA order = 2, number of filters = N = 4)

Arma_A = [1.0000   -1.8996    0.9025 ;...
    1.0000   -1.9075    0.9409 ;...
    1.0000   -1.8471    0.8649 ;...
    1.0000   -1.7642    0.8464];

Arma_B = [1.0000   -1.9027    0.9409 ;...
    1.0000   -1.8959    0.9025 ;...
    0.700   -0.9176    0.4512 ;...
    0.300   -0.1836    0.0846];


% Computation of the kernels
g_theo = zeros(N,len_IR);
for k=1:N
    g_theo(k,:) = filter(Arma_B(k,:),Arma_A(k,:),[1 zeros(1,len_IR-1)]);
end
G_theo = fft(g_theo.').';

axe_f = ((1:len_IR)-1)./(len_IR-1).*fs;  % frequency axes

Nonlinear System response

y = 0 ;
for n = 1:N
    y = y + filter(Arma_B(n,:),Arma_A(n,:),[x.^n zeros(1,len_IR)]) ;
end

% figure; plot(y);
% title('Response of the Nonlinear System to the Sychronized Swept-Sine Signal');

Nonlinear convolution

y = y - mean(y);  % signal must not contain the mean value

% Nonlinear convolution in frequency domain:
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);

% definition of the inferse filter in spectral domain
% (A.Novak IEEE paper "Nonlinear System Identification Using Exponential
% Swept-Sine Signal",  Eq. (30) and (31) leading to Eq. (26) for the amplitude
% and (9) for the phase)
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

h = ifft(H,'symmetric');     % higher order impulse responses

Separation of higher impulse responses

(positions of higher harmonic IRs)

dt = L.*log(1:N).*fs;

% rounded positions [samples]
dt_ = round(dt);
% non-integer difference
dt_rem = dt - dt_;

% circular periodisation of IR
h_pos = [h h(1: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(N,len_IR);

st0 = length(h_pos);
for n = 1:N
    % start position of n-th IR
    st = len-dt_(n)-pre_IR;
    ed = min(st + len_IR-1, st0);
    st0 = st - 1;

    % separation of IRs
    hs(n,1:ed-st+1) = h_pos(st:ed);

    % Non-integer sample delay
    Hx = fft(hs(n,:)).* ...
        exp(-j*dt_rem(n)*axe_w);
    hs(n,:) = ifft(Hx,'symmetric');
end




% Higher Harmonics
Hs = fft(hs.').';

figure; axes('FontSize',12);
semilogx(axe_f,20*log10(abs(Hs)),'linewidth',2)
set(gca,'XLim',[n*f1 n*f2]);
grid on;

xlabel('Frequency [kHz]')
ylabel('Amplitude [-]')
title(['Higher Harmonics'])
ylim([-40 40]);

Conversion to Generalized Hammerstein model

% Matrix A generation
A(N,N) = 0;
for l = 1:N
    for k = 1:N
        if ((l>=k) && (mod(l+k,2)==0))
            A(k,l) = (-1)^(2*l+(1-k)/2)/ ...
                2^(l-1)*nchoosek(l,(l-k)/2);
        end
    end
end

% transformation
Gs = inv(A)*Hs;

% ifft with Hermitian property
gs = ifft(Gs.','symmetric').';
Gs = fft(gs.').';


Gs_roll = fft([gs(:,pre_IR+2:end) gs(:,1:pre_IR+1)].').';


for n = 1:N

    figure; axes('FontSize',12);
    semilogx(axe_f,abs(G_theo(n,:)),'-r','linewidth',2)
    hold on
    semilogx(axe_f,abs(Gs(n,:)),'-.b','linewidth',2)
    set(gca,'XLim',[N*f1 n*f2]);

    legend('Theoretical','Estimated')
    xlabel('Frequency [kHz]')
    ylabel('Amplitude [-]')
    title(['Hammerstein kernel (order ' num2str(n) ')'])
    grid on;

end


for n = 1:N

    figure; axes('FontSize',12);
    semilogx(axe_f,angle(G_theo(n,:)),'-r','linewidth',2)
    hold on
    semilogx(axe_f,angle(Gs_roll(n,:)),'-.b','linewidth',2)
    set(gca,'XLim',[N*f1 n*f2]);

    legend('Theoretical','Estimated')
    xlabel('Frequency [kHz]')
    ylabel('Phase [rad]')
    title(['Hammerstein kernel (order ' num2str(n) ')'])
    grid on;

end



© 2014 Antonin Novak - homepage (All Rights Reserved)