# Identification of Nonlinear Systems in Series

# Measurement on an Acoustical wave-guide

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

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.

## MATLAB® Code

The whole matlab code with all the functions and the measured_data.mat file can be downloaded here: Nonlinear_System_Identification_Waveguide.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

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

where is the rate of frequency increase defined as

## 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 figure; semilogx(20*log10(abs(M1.HHFRF(1:3,:).')./2e-5),'linewidth',2); 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; figure; semilogx(20*log10(abs(M2.HHFRF(1:3,:).')./2e-5),'linewidth',2); 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); end 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); end % estimate the filters G(:,frs+1) = pinv(XX)*M2.HHFRF(:,frs+1); end

## 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); end % plot results figure; for k=1:Nbb semilogx(f_axis(k*f1u:f2u),20*log10(abs(synth_HHFRF(k,k*f1u:f2u))./2e-5),'linewidth',2); hold all end hold off hold on for k=1:Nbb semilogx(k*fx,20*log10(abs(theo_HHFRF(k,:))./2e-5),'*'); hold all end 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 end 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); end % plot results figure; co = 342.1; % speed of sound s = round(dist/co*96000); % delay due to the propagation pts = round(3*fs/fo); % points to plot plot((0:pts-1)./fs*1000,y((1:pts)+s),'linewidth',2); hold all plot((0:pts-1)./fs*1000,y_burgers(1:pts),'r--','linewidth',2); 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;