System response from input and output signals

Question

suniya vs 2023-8-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2010877-system-response-from-input-and-output-signals

编辑： Star Strider 2023-8-25

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-8-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2010877-system-response-from-input-and-output-signals#answer_1290457

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.

12 个评论
显示 10更早的评论隐藏 10更早的评论

Star Strider 2023-8-24

在 MATLAB Online 中打开

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

.

Star Strider 2023-8-24

编辑：Star Strider 2023-8-24

在 MATLAB Online 中打开

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 2023-8-25

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 2023-8-25

编辑：Star Strider 2023-8-25

在 MATLAB Online 中打开

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.

.

请先登录，再进行评论。

System response from input and output signals

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

12 个评论
显示 10更早的评论隐藏 10更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

System response from input and output signals

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

12 个评论 显示 10更早的评论隐藏 10更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论