System response from input and output signals

My input is chirp signal .How could i get the frequency response of the sysytem from the output . Using this i want to create a database

Answer 1

[Star Strider]

First, use the System Identification Toolbox (iddata then ssest and compare then lsim) or the Signal Processing Toolbox (invfreqz then filtfilt) to estimate the system and use your chirp function as the input for the simulation.

[suniya vs]

My system is non linear. There is some harmonics in the output signal. Using the above method i cant replicate my output signal. Simulated signal does not contain harmonics.

[suniya vs]

clear all;close all;clc;

[inpSig,fs]=audioread('varFsin_1K.wav');

[outSig, Fs] = audioread('AliELENvarFsin_1K261104b500_8_23.wav');

windowLength = 256;

win = hamming(windowLength,"periodic");

overlap = round(0.75 * windowLength);

ffTLength = windowLength;

Ts=1;

data = iddata(outSig,inpSig,Ts)

nx = 10;

sys = ssest(data,nx);

compare(data,sys)

dt = 1/Fs;

t=0:dt:(length(inpSig)*dt)-dt;

z=lsim(sys,inpSig,t)

spectrogram(z,win,overlap,ffTLength,Fs,'yaxis');

[Star Strider]

It should be possible to approximate a linear system with the input and output signals. It may be necessary to increase the system order.

I need the signals themselves to experiment with them in order to identify the system. I do not see them posted.

[suniya vs]

Thank you for the reply.i am not able to attach .wav files here

[suniya vs]

Thank you for the reply.i am not able to attach .wav files here

[Star Strider]

As always, my pleasure!

You can use the zip function to zip them and then upload the .zip file.

[suniya vs]

Pls see the attache files

[Star Strider]

Looking at the input and output signals and their spectral analyses, I am surprised that the results are as good as they are. You can continue increasing the system order, however I doubt the results will improve significantly. I doubt that the system you are modeling is actually linear.

My analysis —

Uz1 = unzip('AliELENvarFsin...b500_8_23.zip')

Uz1 = 1×1 cell array
    {'AliELENvarFsin_1K261104b500_8_23.wav'}

Uz2 = unzip('varFsin_1K.zip')

Uz2 = 1×1 cell array
    {'varFsin_1K.wav'}

[inpSig,fs]=audioread(Uz2{1});
size_inp = size(inpSig)

size_inp = 1×2

20001 1

fs

fs = 2000

[outSig, Fs] = audioread(Uz1{:});
size_out = size(outSig)

size_out = 1×2

20001 1

Fs

Fs = 2000

windowLength = 256;
win = hamming(windowLength,"periodic");
overlap = round(0.75 * windowLength);
ffTLength = windowLength;
Ts=1/Fs;
data = iddata(outSig,inpSig,Ts)

data =

Time domain data set with 20001 samples.
Sample time: 0.0005 seconds              
                                         
Outputs      Unit (if specified)         
   y1                                    
                                         
Inputs       Unit (if specified)         
   u1                                    
                                         

nx = 10;
sys = ssest(data,nx);
compare(data,sys)

dt = 1/Fs;
t=0:dt:(length(inpSig)*dt)-dt;
z=lsim(sys,inpSig,t)

z = 20001×1

1.0e+00 * 0 0 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

spectrogram(z,win,overlap,ffTLength,Fs,'yaxis');

opts = bodeoptions;
opts.FreqUnits = 'Hz';

figure
bodeplot(sys, opts)
grid

[p1,f1] = pspectrum(inpSig,fs);
[p2,f2] = pspectrum(outSig,Fs);

figure
yyaxis left
plot(f1, p1)
ylabel('Input')
yyaxis right
plot(f2, p2)
ylabel('Output')
grid
xlabel('Frequency (Hz)')

figure
pspectrum(inpSig,fs,'spectrogram')
colormap(turbo)

figure
pspectrum(outSig,Fs,'spectrogram')
colormap(turbo)

L = numel(inpSig);
t = linspace(0, L-1, L)/fs;

[FTs1,Fv] = FFT1(inpSig,t);
[FTs2,Fv] = FFT1(outSig,t);
transfcn = FTs2 ./ FTs1;

figure
plot(Fv, mag2db(abs(transfcn)))
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
title('Transfer Function')

function [FTs1,Fv] = FFT1(s,t)
    s = s(:);
    t = t(:);
    L = numel(t);
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
    Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv);
end

.

[suniya vs]

Sir, but there is a lot of harmonics in the output signal 'outsig'. It is not in the output of simulated system 'z'.

[Star Strider]

A linear system may not have the ability to model all of the harmonics. I usually prefer ssest for systems, however here tfest may be more appropriate, although the system wil have a very high order (32 is as high as I can get here without the code timing out) so experiment with higher orders —

Uz1 = unzip('AliELENvarFsin...b500_8_23.zip')

Uz1 = 1×1 cell array
    {'AliELENvarFsin_1K261104b500_8_23.wav'}

Uz2 = unzip('varFsin_1K.zip')

Uz2 = 1×1 cell array
    {'varFsin_1K.wav'}

[inpSig,fs]=audioread(Uz2{1});
size_inp = size(inpSig)

size_inp = 1×2

20001 1

fs

fs = 2000

[outSig, Fs] = audioread(Uz1{:});
size_out = size(outSig)

size_out = 1×2

20001 1

Fs

Fs = 2000

windowLength = 256;
win = hamming(windowLength,"periodic");
overlap = round(0.75 * windowLength);
ffTLength = windowLength;
Ts=1/Fs;
data = iddata(outSig,inpSig,Ts)

data =

Time domain data set with 20001 samples.
Sample time: 0.0005 seconds              
                                         
Outputs      Unit (if specified)         
   y1                                    
                                         
Inputs       Unit (if specified)         
   u1                                    
                                         

nx = 32;
sys = tfest(data,nx);

figure
compare(data,sys)

dt = 1/Fs;
t=0:dt:(length(inpSig)*dt)-dt;
z=lsim(sys,inpSig,t)

z = 20001×1

1.0e+00 * 0 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000

spectrogram(z,win,overlap,ffTLength,Fs,'yaxis');

opts = bodeoptions;
opts.FreqUnits = 'Hz';

figure
bodeplot(sys, opts)
grid

[p1,f1] = pspectrum(inpSig,fs);
[p2,f2] = pspectrum(outSig,Fs);

figure
yyaxis left
plot(f1, p1)
ylabel('Input')
yyaxis right
plot(f2, p2)
ylabel('Output')
grid
xlabel('Frequency (Hz)')

figure
pspectrum(inpSig,fs,'spectrogram')
colormap(turbo)

figure
pspectrum(outSig,Fs,'spectrogram')
colormap(turbo)

L = numel(inpSig);
t = linspace(0, L-1, L)/fs;

[FTs1,Fv] = FFT1(inpSig,t);
[FTs2,Fv] = FFT1(outSig,t);
transfcn = FTs2 ./ FTs1;

figure
plot(Fv, mag2db(abs(transfcn)))
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
title('Transfer Function')

function [FTs1,Fv] = FFT1(s,t)
    s = s(:);
    t = t(:);
    L = numel(t);
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
    Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv);
end

.

[suniya vs]

Sir ,There is some harmonics in the system output but it isnit in the simulated output. i want a system that gives similar output.

[Star Strider]

The harmonics are probably there. You will need to increase the system order beyond 32 to get an accurate representation, although modeling all of them may not be possible. (It turns out that tfest models this system better than ssest, so use tfest instead.)

EDIT — (25 Aug 2023 at 14:37)

Another approach (signal processing) is to use the firls function to approximate the transfer function —

Uz1 = unzip('AliELENvarFsin...b500_8_23.zip')

Uz1 = 1×1 cell array
    {'AliELENvarFsin_1K261104b500_8_23.wav'}

Uz2 = unzip('varFsin_1K.zip')

Uz2 = 1×1 cell array
    {'varFsin_1K.wav'}

[inpSig,fs]=audioread(Uz2{1});
size_inp = size(inpSig)

size_inp = 1×2

20001 1

fs

fs = 2000

[outSig, Fs] = audioread(Uz1{:});
size_out = size(outSig)

size_out = 1×2

20001 1

Fs

Fs = 2000

L = numel(inpSig);
t = linspace(0, L-1, L)/fs;

[FTs1,Fv] = FFT1(inpSig,t);
[FTs2,Fv] = FFT1(outSig,t);
transfcn = FTs2 ./ FTs1;

figure
plot(Fv, mag2db(abs(transfcn)))
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
title('Transfer Function')

n = 350;
f = Fv(1:end-1);
a = abs(transfcn(1:end-1));
b = firls(n, f/max(f), a);

[hz,fz] = freqz(b, 1, 2^16, Fs);

figure
semilogy(f, a, 'DisplayName','Transfer Function')
hold on
plot(fz, abs(hz), 'DisplayName',["FIR Filter (Order "+string(n)+")"])
hold off
grid
legend('Location','best')

function [FTs1,Fv] = FFT1(s,t)
    s = s(:);
    t = t(:);
    L = numel(t);
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
    Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv);
end

Use filtfilt to implement the filter.

.