Antonin Novak - homepage

Postdoctoral Research Fellow

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

Email Address: ant.novak[at]

Antonin Novak

Identification of Nonlinear Systems in Series

Measurement on an Acoustical wave-guide

(The whole matlab code can be downloaded here:

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 (acoustical wave-guide) under test, which is preceded by another nonlinear system (amplifier + loudspeaker). The whole measured system consists thus of two nonlinear systems in series. The first one is the source of the acoustical signal (amplifier + loudspeaker) which is supposed to be linear, but (as shown below) is nonlinear and contaminates the measurements of the second systems' non-linearities. The goal of the method is to identify separately the second nonlinear system (acoustical wave-guide).

The method is based on the estimation of the Higher Harmonic Frequency Responses (HHFRFs) from the measurement of distorted input and output signals (microphone M1 and microphone M2). The second nonlinear system is then modeled by non-parametric generalized Hammerstein model made up of power series associated with linear filters.


The whole matlab code with all the functions and the measured_data.mat file can be downloaded here:

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

The synchronized swept-sine signal is genereted using the following formula:

$$ x(t) = \sin \left\{ 2 \pi f_1 L \left[ \exp \left( \frac{t}{L} \right)
-1 \right] \right\}, $$

where $L$ is the rate of frequency increase defined as

$$ L = \frac{1}{f_1} \texttt{Round} \left[ \frac{\hat{T} f_1}{\ln \left( \frac{f_2}{f_1} \right)} \right],$

Measurement on the acoustical wave-guide

the input signal "swept_sine" is send to the system under test (amplifier - loudspeaker - wave-guide). Two microphones M1 and M2 are flush mounted on the wave-guide wall. The first one is set at 0.2 m from the source, the second one is set at the distance of 4.2 m from the source. Their signals (responses to the swept-sine signal) are loaded in the following code as M1 and M2

load measured_data.mat M1 M2;

Higher Harmonic Freqency Responses (HHFRF)

N=9; % number of HHFRFs

% useful frequencies
f1u = 400;
f2u = 4.9e3;

% microphone 1 (M1) 20 cm from the source
[M1.h, M1.dt] = nonlin_conv(M1.signal,swept_sine,N);
M1.HHFRF = separation(M1.h,swept_sine,N,M1.dt,1,1800);

% microphone 2 (M2) 4.2 m from the source
[M2.h, M2.dt] = nonlin_conv(M2.signal,swept_sine,N);
M2.HHFRF = separation(M2.h,swept_sine,N,M2.dt,1000,2900);

% plot HHFRFs
xlim([f1u f2u]); ylim([100 160]);
title([{'HHFRFs at the first microphone (M1),'};...
       {'at the distance of 0.2 m from the source'}]);
xlabel('frequency [Hz]');
ylabel('acoustic pressure [dB SPL]');
legend({'1st harmonic','2nd harmonic','3rd harmonic'},'location','southwest');
grid on;

xlim([f1u f2u]); ylim([100 160]);
title([{'HHFRFs at the second microphone (M2),'};...
       {'at the distance of 4.2 m from the source.'}]);
xlabel('frequency [Hz]');
ylabel('acoustic pressure [dB SPL]');
legend({'1st harmonic','2nd harmonic','3rd harmonic'},'location','southwest');
grid on;

Identification of the nonlinear system NL2 between M1 and M2 (the wave-guide alone)

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

% HHFRFs of the powers of the signal at M1
for n=1:Nbb
    [h, dt] = nonlin_conv(M1.signal.^n,swept_sine,N);
    M1.HHFRF_NL{n} = separation(h,swept_sine,N,dt,1,1800);

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 = f1u:f2u
    % prepare the harmonics for each branch
    for n=1:Nbb
        XX(:,n) = M1.HHFRF_NL{n}(:,frs+1);
    % estimate the filters
    G(:,frs+1) = pinv(XX)*M2.HHFRF(:,frs+1);

Comparison with the Burgers' theory in frequency domain

% # Synthesized HHFRFs of the wave-guide (estimated NL model) for a given input pressure
p0 = 650; % input acoustic pressure [Pa]
C = powersin_ampl(Nbb,p0); % coefficients of the matrix fot power-to-harmonics conversion
synth_HHFRF = C.'*G; % synthesized HHFRFs

% # HHFRF based on the Burgers Equation
dist = M2.distance - M1.distance ; % distance between the fisrt and the second micropone
R=0.057/2 ; % radius of the wave-guide
fx = round(logspace(log10(f1u),log10(f2u),20)); % frequencies
for k = 1:length(fx)
    theo_HHFRF(:,k) = propag_burgers(Nbb, fx(k), p0, dist, R, 0.1);

% plot results
for k=1:Nbb
    hold all
hold off
hold on
for k=1:Nbb
    hold all
xlim([f1u f2u]); ylim([100 160]);
title([{'Synthesized HHFRFs at the second microphone (M2) compared with the theory.'};...
       {'The input signal at first microphone (M1) is a sine-wave (650 Pa).'}]);
xlabel('frequency [Hz]');
ylabel('acoustic pressure [dB SPL]');
legend({'1st harmonic (estimated)','2nd harmonic (estimated)','3rd harmonic (estimated)',...
        '1st harmonic (theory)','2nd harmonic (theory)','3rd harmonic(theory)'},'location','southwest');
grid on;

Comparison with the Burgers' theory in time domain

% Synthesis: harmoic sound pressure at M1
fo = 1400;               % inpit frequency
t = 0:1/fs:1-1/fs;       % time axis (1 second)
x = p0*sin(2*pi*fo*t);   % input signal

Y = 0;                   % prepare the output spectra
for n=1:Nbb              % for each branche of the Generalized Hammerstein model
    Y = Y + fft(x.^n).*G(n,:); % calculate the output spectra
y = ifft(Y,'symmetric'); % output signal (distorted)

% Burgers equation solution
YB = propag_burgers(Nbb, fo, p0, dist, R, 0.1); % fourier coefficinets
y_burgers = 0;  % prepare the output signal
for k=1:Nbb
    y_burgers = y_burgers + real(YB(k))*sin(k*2*pi*fo*t) + imag(YB(k))*cos(k*2*pi*fo*t);

% plot results
co = 342.1; % speed of sound
s = round(dist/co*96000);  % delay due to the propagation
pts = round(3*fs/fo); % points to plot
hold all
xlim([0 (pts-1)/fs*1000]);
title('Time waveforms of the acoustical pressure at M2 for an input signal at 1400 Hz');
xlabel('time [ms]');
ylabel('acoustic pressure [Pa]');
legend({'re-synthesized waveform from the estimated model',...
        'theoretical waveform calculated using the Burger’s theory'},'location','southwest');
grid on;

© 2014 Antonin Novak - homepage (All Rights Reserved)