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





Identification of Nonlinear Systems in Series

Simulation

(The whole matlab code can be downloaded here: Nonlinear_System_Identification_Simulation.zip)

Published in: A. Novak, B. Maillou, P. Lotton & L. Simon (2014) "Nonparametric Identification of Nonlinear Systems in Series", Instrumentation and Measurement, IEEE Transactions on. Vol. 63(8), pp. 2044-2051. (download)



The following MATLAB® code is a tool for the identification of a nonlinear system, which is preceded by another nonlinear system. The goal is to identify separately the second nonlinear system.

The method is based on the estimation of the Higher Harmonic Frequency Responses (HHFRFs) of distorted signals "u" and "y". The second nonlinear system is then modelled by non-parametric generalized Hammerstein model made up of power series associated with linear filters.



MATLAB® Code

The whole matlab code with all the functions can be downloaded here: Nonlinear_System_Identification_Simulation.zip

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

Synchronized Swept-Sine generation

swept_sine.f1 = 20; % initial frequency [Hz]
swept_sine.f2 = 5000; % final frequency [Hz]
swept_sine.fs = fs; % sampling frequency [Hz]
swept_sine.T_ = 15; % approximative time duration [s]
swept_sine.T_fade = [0.3000 0.1000]; % fade-in, fade-out [s]


swept_sine = swept_gene(swept_sine); % swept-sine signal generation

definition of the fisrt nonlinear system (NL1)

A1 = [1 -1.7 0.8];
B1 = [0.13 -0.1 0.03];

[H_1,W_1] = freqz(B1,A1);
u = atan(3*filter(B1,A1,swept_sine.s));


% plot the first nonlinear system characteristics
figure('Position',[243 474 1082 504]);

ax1 = subplot(121);
set(ax1,'FontSize',12,'FontName','Century','FontWeight','Bold');
plot(W_1/pi*fs/2/1000,20*log10(abs(H_1)),'linewidth',2); grid on
set(gca,'XLim',[swept_sine.f1 swept_sine.f2]./1000);
title('Frequency Response Function of the 1st NL system');
xlabel('Frequency [kHz]');
ylabel('Amplitude [dB]');

ax2 = subplot(122);
set(ax2,'FontSize',12,'FontName','Century','FontWeight','Bold');
plot(-1:0.01:1,atan(3*[-1:0.01:1]),'linewidth',2); grid on
title('Static nonlinearity of the 1st NL system');
xlabel('Input');
ylabel('Output');


% add a noise at the output of the 1st nonlinear system
SNRdB = 50;
A_nois = sqrt(mean(u.^2))/(10^(SNRdB/20));
u = u + A_nois.*randn(size(u,1),size(u,2));

definition of the second nonlinear system (NL2)

order = 3;

B = [1    0.1  0.5;
    -1.9 -0.3 -0.9;
     0.9  0.1  0.4];

A = [1    1    1;
    -1.9 -1.6 -1.9;
     0.95 0.8 0.9];

y = 0;
for k=1:order
    y = y + filter(B(:,k),A(:,k),u.^k);
end

% add a noise at the output of the second nonlinear system
SNRdB2 = 50;

A_nois = sqrt(mean(y.^2))/(10^(SNRdB2/20));
y = y + A_nois.*randn(size(y,1),size(y,2));

IDENTIFICATION PROCEDURE (of the second nonlinear system)

N = 9; % number of HHFRs
Nbb = 3; % number of branches of the Generalized Hammerstein model

% HHFRs of the powers of the signal x (after the NL1)
for k=1:Nbb
    [h, dt] = nonlin_conv(u.^k,swept_sine,N);
    u_HHFRF{k} =separation(h,swept_sine,N,dt);
end

% HHFRs of the signal y (after the NL2)
[h, dt] = nonlin_conv(y,swept_sine,N);
y_HHFRF = separation(h,swept_sine,N,dt);



f_axis = 0:fs-1; % frequency axis

% estimate the filters G(f) of the Generalized Hammerstein model
% (each frequency individually)
G = zeros(Nbb,length(f_axis));
for frs = swept_sine.f1:Nbb*swept_sine.f2
    % prepare the harmonics for each branch
    for n=1:Nbb
        XX(:,n) = u_HHFRF{n}(:,frs+1);
    end
    % estimate the filters
    G(:,frs+1) = pinv(XX)*y_HHFRF(:,frs+1);
end


for k=1:Nbb
    [H_2,W_2] = freqz(B(:,k),A(:,k),f_axis/fs*pi*2);

    figure;
    subplot(211)
    plot(f_axis(Nbb*swept_sine.f1:k*swept_sine.f2)./1000,...
        20*log10(abs(G(k,Nbb*swept_sine.f1:k*swept_sine.f2).')),...
        'linewidth',2); hold on
    plot(W_2(Nbb*swept_sine.f1:k*swept_sine.f2)/pi*fs/2/1000,...
        20*log10(abs(H_2(Nbb*swept_sine.f1:k*swept_sine.f2))),...
        'r--','linewidth',2); hold off; grid on
    legend('Estimated','Theoretical','location','best');
    xlim([swept_sine.f1 swept_sine.f2]./1000);
    xlabel('Frequency [kHz]');
    ylabel('Amplitude [dB]');
    title(['G_' num2str(k)]);

    subplot(212)
    plot(f_axis(Nbb*swept_sine.f1:k*swept_sine.f2)./1000,...
        angle(G(k,Nbb*swept_sine.f1:k*swept_sine.f2).'),...
        'linewidth',2); hold on
    plot(W_2(Nbb*swept_sine.f1:k*swept_sine.f2)/pi*fs/2/1000,...
        angle(H_2(Nbb*swept_sine.f1:k*swept_sine.f2)),...
        'r--','linewidth',2); hold off; grid on
    xlim([swept_sine.f1 swept_sine.f2]./1000);
    legend('Estimated','Theoretical','location','best');
    xlabel('Frequency [kHz]');
    ylabel('Angle [rad]');
end



© 2014 Antonin Novak - homepage (All Rights Reserved)