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

# Simulation

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 = 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]');