Evaluate the system impose response function and the FRF with given input and output data

Question

Christina Reid 2021-1-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/707098-evaluate-the-system-impose-response-function-and-the-frf-with-given-input-and-output-data

评论： Christina Reid 2021-1-4

beam_experiment.mat

Hello!

I have the input and output data of a signal that is stored in a .mat file. I need help calulating the impulse response and the FRF. I have attached the .mat file. I know how to pull the info from the .mat file and the time length is 20 sec

data = importdata ('beam_experiment.mat');

x = transpose(data.x); %input

y = transpose(data.y); %output

fs = data.fs; % sampling frequency

T1 = 20;

I tried the follow:

data = iddata(y, x, 1/fs);

h = impulseest(data,2*max(time)*fs, min(time)*fs);

H = fft(h);

but get the following error message "Error using fft

Invalid data type. First argument must be double, single, int8, uint8, int16, uint16, int32,

uint32, or logical."

Could somone kindly help out?

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

请先登录，再进行评论。

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

Answer 1

Mathieu NOE 2021-1-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/707098-evaluate-the-system-impose-response-function-and-the-frf-with-given-input-and-output-data#answer_589403

在 MATLAB Online 中打开

HELLO

had to find another way than using the ID Toolbox (as I don't have it)

so I basically tried to fit an IIR model to the FRF and then you can create the impulse response

data = importdata ('beam_experiment.mat');
x = transpose(data.x); %input
y = transpose(data.y); %output
fs = data.fs;   % sampling frequency
T1 = 20;
NFFT = 2048;
NOVERLAP = 0.75*NFFT;
[Txy,F] = tfestimate(x,y,hanning(NFFT),NOVERLAP,NFFT,fs);
% IIR model
W = linspace(0,pi,length(F));
ind = find(F>5 & F <80); % frequency weighting ; data reliable between 5 and 80 Hz
WT = zeros(size(W));
WT(ind) = 1;
NB = 6;
NA = NB+2;
[B,A] = invfreqz(Txy,W,NB,NA,WT,1e4);
% check bode plots
[H,WW] = freqz(B,A,length(F));
figure(1),
subplot(2,1,1),plot(F,20*log10(abs(Txy)),'b',F,20*log10(abs(H)),'r');grid
subplot(2,1,2),plot(F,180/pi*(angle(Txy)),'b',F,180/pi*(angle(H)),'r');grid
% Impulse response 
[IR,X] = dimpulse(B,A);
samples = length(IR);
time = 1/fs*(1:samples);
figure(2),
plot(time,IR,'r');grid
title('Impulse response')
xlabel('Time (s)');
ylabel('Amplitude')

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

Christina Reid 2021-1-2

在 MATLAB Online 中打开

Screen Shot 2021-01-02 at 1.24.56 PM.png

Hi Mathieu,

Thank you so much. I was working on the following problem first ( see the 'Screen Shot 2021 - 01-02 at 1.24.56 PM.png' attachment) and the code below. With this problem, I am getting graphs that do not make any sense. From the problem, A1, A2, f1, and f2 are given and I need to find the correct fs value to correctly analyze the plots. I set fs= 100 ( which I thought would be correct since it is more than 2x the max frequency, but the plots are not looking correct.

I then tried to use this code to approach my first problem ( where instead of given the system paramaeters and impulse response function I had my input/ output data) I thought I could use impulseest function to find the impluse response, however I will use your code and see what I get.

In the meantime, would you mind looking at the following code and recommend a fs value?

% boot commands:
clear all; close all;clc 
% Program:
% Define the system parameters:
A1 = 400; 
A2 = 60; 
f1 = 25; 
f2 = 35; 
wn1 = 2*pi*f1; 
wn2 = 2*pi*f2; 
zeta1 = 0.001; 
zeta2 = 0.003;
wd1 = sqrt(1-zeta1^2)*wn1; 
wd2 = sqrt(1-zeta2^2)*wn2;
% Define the acquisition parameters:
fs = 100; 
T1 = 10; 
t1 = 0:1/fs:T1-1/fs;
% Evaluate the system impulse response function h(t) and the FRF H(f)
h = (A1/wd1)*exp(-zeta1*wn1*t1).*sin(wd1*t1)+(A2/wd2)*exp(-zeta2*wn2*t1).*sin(wd2*t1);
H = fft(h);
% Generate the input signal x(t), as white noise:
T = 10000;
x = randn(1,T*fs); 
% Generate the output signal y(t) as the system response to the input
y = filter(h,1,x); % we do not scale for convenience.
% Generate the measuremet noise acting on input and output for case (c)
nx = 0.5*randn(size(x)); % nx=0 for case (a)
ny = 0.5*randn(size(y)); % ny=0 for case (b)
% Calculate the (one-sided) spectral density functions using the segment
% averaging method with hanning window, 50% overlap and 1000 segments. 
% Then calculate the FRF and the coherence function.  
disp('Wait patiently...');
% Case a)
ya = y+ny;
[Gxx, f] = cpsd(x(1:T*fs),x(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(ya(1:T*fs),ya(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(x(1:T*fs),ya(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(ya(1:T*fs),x(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
H1a = Gxy./Gxx; 
H2a = Gyy./Gyx;
HTa = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammaa = abs(Gxy).^2./(Gxx.*Gyy);
disp('                ...still waiting?');
% case b)
xb = x+nx; 
[Gxx, f] = cpsd(xb(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(y(1:T*fs),y(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(y(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
H1b = Gxy./Gxx; 
H2b = Gyy./Gyx;
HTb = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammab = abs(Gxy).^2./(Gxx.*Gyy);
disp('                                  Time for a coffee?');
disp('                                                        S')
disp('                                                       __s__')
disp('                                                      |     |>')
disp('                                                      \_____/')
% case c)
% xc=xb; yc=ya;
[Gxx, f] = cpsd(xb(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(ya(1:T*fs),ya(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(ya(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
H1c = Gxy./Gxx; 
H2c = Gyy./Gyx;
HTc = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammac = abs(Gxy).^2./(Gxx.*Gyy);
clc;disp('Espresso');
% Plot of the results
figure
subplot(2,1,1),plot(f,20*log10(abs(H1a)), f,20*log10(abs(H2a)),f,20*log10(abs(HTa)),f,20*log10(abs(H(1:length(f)))),':')
title('case a: output noise only')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 25 -35 25])
subplot(2,1,2),plot(f,-unwrap(angle(H1a)), f,-unwrap(angle(H2a)),f,-unwrap(angle(HTa)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
axis([0 25 -4 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
subplot(2,1,1),plot(f,20*log10(abs(H1b)), f,20*log10(abs(H2b)),f,20*log10(abs(HTb)),f,20*log10(abs(H(1:length(f)))),':')
title('case b: input noise only')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 25 -35 25])
subplot(2,1,2),plot(f,-unwrap(angle(H1b)), f,-unwrap(angle(H2b)),f,-unwrap(angle(HTb)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
axis([0 25 -4 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
subplot(2,1,1),plot(f,20*log10(abs(H1c)), f,20*log10(abs(H2c)),f,20*log10(abs(HTc)),f,20*log10(abs(H(1:length(f)))),':')
title('case c: both input and output noise')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 25 -35 25])
subplot(2,1,2),plot(f,-unwrap(angle(H1c)), f,-unwrap(angle(H2c)),f,-unwrap(angle(HTc)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
axis([0 25 -4 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
plot(f, Gammaa,f,Gammab,f,Gammac); 
xlabel('Frequency (Hz)'); ylabel('Coherence function')
legend('case a','case b','case c')
axis([0 25 0 1])
clc
% END

Mathieu NOE 2021-1-2

在 MATLAB Online 中打开

hello again Christina

I simply commented the lines with caxis

because you force the x axis in the 0 to 25 Hz range, why ?

also , for my own interest, for case a , I plotted the "true" FRF , computed from the given imulse response

see

f,20*log10(abs(Txy)),'-.k',

maybe it can be of some interest for you ,

so the entire code below :

% boot commands:
clear all; close all;clc 
% Program:
% Define the system parameters:
A1 = 400; 
A2 = 60; 
f1 = 25; 
f2 = 35; 
wn1 = 2*pi*f1; 
wn2 = 2*pi*f2; 
zeta1 = 0.001; 
zeta2 = 0.003;
wd1 = sqrt(1-zeta1^2)*wn1; 
wd2 = sqrt(1-zeta2^2)*wn2;
% Define the acquisition parameters:
fs = 100; 
T1 = 10; 
t1 = 0:1/fs:T1-1/fs;
% Evaluate the system impulse response function h(t) and the FRF H(f)
h = (A1/wd1)*exp(-zeta1*wn1*t1).*sin(wd1*t1)+(A2/wd2)*exp(-zeta2*wn2*t1).*sin(wd2*t1);
H = fft(h);
% Generate the input signal x(t), as white noise:
T = 10000;
x = randn(1,T*fs); 
% Generate the output signal y(t) as the system response to the input
y = filter(h,1,x); % we do not scale for convenience.
% Generate the measuremet noise acting on input and output for case (c)
nx = 0.5*randn(size(x)); % nx=0 for case (a)
ny = 0.5*randn(size(y)); % ny=0 for case (b)
% Calculate the (one-sided) spectral density functions using the segment
% averaging method with hanning window, 50% overlap and 1000 segments. 
% Then calculate the FRF and the coherence function.  
disp('Wait patiently...');
% Case a)
ya = y+ny;
[Gxx, f] = cpsd(x(1:T*fs),x(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(ya(1:T*fs),ya(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(x(1:T*fs),ya(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(ya(1:T*fs),x(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
H1a = Gxy./Gxx; 
H2a = Gyy./Gyx;
HTa = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammaa = abs(Gxy).^2./(Gxx.*Gyy);
disp('                ...still waiting?');
% case b)
xb = x+nx; 
[Gxx, f] = cpsd(xb(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(y(1:T*fs),y(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(y(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
H1b = Gxy./Gxx; 
H2b = Gyy./Gyx;
HTb = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammab = abs(Gxy).^2./(Gxx.*Gyy);
disp('                                  Time for a coffee?');
disp('                                                        S')
disp('                                                       __s__')
disp('                                                      |     |>')
disp('                                                      \_____/')
% case c)
% xc=xb; yc=ya;
[Gxx, f] = cpsd(xb(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(ya(1:T*fs),ya(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(ya(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
H1c = Gxy./Gxx; 
H2c = Gyy./Gyx;
HTc = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammac = abs(Gxy).^2./(Gxx.*Gyy);
clc;disp('Espresso');
% "true" FRF
Txy = freqz(h,1,f,fs); % H = FREQZ(...,F,Fs) returns the complex frequency response 
% Plot of the results
figure
subplot(2,1,1),plot(f,20*log10(abs(Txy)),'-.k',f,20*log10(abs(H1a)), f,20*log10(abs(H2a)),f,20*log10(abs(HTa)),f,20*log10(abs(H(1:length(f)))),':')
title('case a: output noise only')
legend('true FRF','|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
% axis([0 25 -35 25])
subplot(2,1,2),plot(f,unwrap(angle(Txy)),'-.k',f,-unwrap(angle(H1a)), f,-unwrap(angle(H2a)),f,-unwrap(angle(HTa)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
% axis([0 25 -4 0.5])
legend('true FRF','\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
subplot(2,1,1),plot(f,20*log10(abs(H1b)), f,20*log10(abs(H2b)),f,20*log10(abs(HTb)),f,20*log10(abs(H(1:length(f)))),':')
title('case b: input noise only')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
% axis([0 25 -35 25])
subplot(2,1,2),plot(f,-unwrap(angle(H1b)), f,-unwrap(angle(H2b)),f,-unwrap(angle(HTb)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
% axis([0 25 -4 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
subplot(2,1,1),plot(f,20*log10(abs(H1c)), f,20*log10(abs(H2c)),f,20*log10(abs(HTc)),f,20*log10(abs(H(1:length(f)))),':')
title('case c: both input and output noise')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
% axis([0 25 -35 25])
subplot(2,1,2),plot(f,-unwrap(angle(H1c)), f,-unwrap(angle(H2c)),f,-unwrap(angle(HTc)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
% axis([0 25 -4 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
plot(f, Gammaa,f,Gammab,f,Gammac); 
xlabel('Frequency (Hz)'); ylabel('Coherence function')
legend('case a','case b','case c')
% axis([0 25 0 1])
clc
% END

Christina Reid 2021-1-4

在 MATLAB Online 中打开

Hello!

So I am actually trying to get the impulse response function, using, ilaplace, and getting the following error: 'Undefined function 'ilaplace' for input arguments of type 'double'

data = importdata ('beam_experiment.mat');
x = (data.x); %input
y = (data.y); %output
fs = data.fs;   % sampling frequency
T1 = 20;
NFFT = 2048;
NOVERLAP = 0.75*NFFT;
G = tfestimate(x,y);
h = ilaplace(Txy)

Above is what I have so far. Below is what I am trying to do. I have capalized a comment line right above where I am getting the error. I have decided to do it this way, because, ideally I would like my code to loook as such ( assuuming that h is properlly defined). The code below is a mix of my previous 2 questions.

data = importdata ('beam_experiment.mat');
x = (data.x); %input
y = (data.y); %output
fs = data.fs;   % sampling frequency
T1 = 20;
G = tfestimate(x,y);
%% HERE IS WHERE I AM TRYING TO DEFINE MY IMPULSE RESPONSE FUNCTION, BUT GETTING THE ERROR
h = ilaplace(Txy)
H = fft(h);
T = 10000;
% Generate the measuremet noise acting on input and output for case (c)
nx = 0.5*randn(size(x)); % nx=0 for case (a)
ny = 0.5*randn(size(y)); % ny=0 for case (b)
% Case a)
ya = y+ny;
[Gxx, f] = cpsd(x(1:T*fs),x(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(ya(1:T*fs),ya(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(x(1:T*fs),ya(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(ya(1:T*fs),x(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
H1a = Gxy./Gxx; 
H2a = Gyy./Gyx;
HTa = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammaa = abs(Gxy).^2./(Gxx.*Gyy);
disp('                ...still waiting?');
% case b)
xb = x+nx; 
[Gxx, f] = cpsd(xb(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(y(1:T*fs),y(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(y(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
H1b = Gxy./Gxx; 
H2b = Gyy./Gyx;
HTb = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammab = abs(Gxy).^2./(Gxx.*Gyy);
disp('                                  Time for a coffee?');
disp('                                                        S')
disp('                                                       __s__')
disp('                                                      |     |>')
disp('                                                      \_____/')
% case c)
% xc=xb; yc=ya;
[Gxx, f] = cpsd(xb(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(ya(1:T*fs),ya(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(ya(1:T*fs),xb(1:T*fs), hanning(T1*fs),T1*fs/2, T1*fs, fs);
H1c = Gxy./Gxx; 
H2c = Gyy./Gyx;
HTc = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammac = abs(Gxy).^2./(Gxx.*Gyy);
clc;disp('Espresso');
% Plot of the results
figure
subplot(2,1,1),plot(f,20*log10(abs(H1a)), f,20*log10(abs(H2a)),f,20*log10(abs(HTa)),f,20*log10(abs(H(1:length(f)))),':')
title('case a: output noise only')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
subplot(2,1,2),plot(f,-unwrap(angle(H1a)), f,-unwrap(angle(H2a)),f,-unwrap(angle(HTa)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
subplot(2,1,1),plot(f,20*log10(abs(H1b)), f,20*log10(abs(H2b)),f,20*log10(abs(HTb)),f,20*log10(abs(H(1:length(f)))),':')
title('case b: input noise only')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 25 -35 25])
subplot(2,1,2),plot(f,-unwrap(angle(H1b)), f,-unwrap(angle(H2b)),f,-unwrap(angle(HTb)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
subplot(2,1,1),plot(f,20*log10(abs(H1c)), f,20*log10(abs(H2c)),f,20*log10(abs(HTc)),f,20*log10(abs(H(1:length(f)))),':')
title('case c: both input and output noise')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 25 -35 25])
subplot(2,1,2),plot(f,-unwrap(angle(H1c)), f,-unwrap(angle(H2c)),f,-unwrap(angle(HTc)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)'))
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
plot(f, Gammaa,f,Gammab,f,Gammac); 
xlabel('Frequency (Hz)'); ylabel('Coherence function')
legend('case a','case b','case c')
clc
% END

Christina Reid 2021-1-4

在 MATLAB Online 中打开

Hi Mathieu,

Thank you for that. I Since I have the input/output data, I decided to find h by doing y./x and then using the reaming code to get the plots for case a, case b, case c, but there appears to be a lot of noise in the plots:

data = importdata ('beam_experiment.mat');
x = (data.x); %input
y = (data.y); %output
fs = data.fs;   % sampling frequency
T1 = 20;
t1 = 0:1/fs:T1-1/fs;
H = y./x;
%H = fft(h);
T = 20;
T2 = 10;
% Generate the measuremet noise acting on input and output for case (c)
nx = 0.5*randn(size(x)); % nx=0 for case (a)
ny = 0.5*randn(size(y)); % ny=0 for case (b)
% Case a)
ya = y+ny;
[Gxx, f] = cpsd(x(1:T*fs),x(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(ya(1:T*fs),ya(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(x(1:T*fs),ya(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(ya(1:T*fs),x(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
H1a = Gxy./Gxx; 
H2a = Gyy./Gyx;
HTa = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammaa = abs(Gxy).^2./(Gxx.*Gyy);
disp('                ...still waiting?');
% case b)
xb = x+nx; 
[Gxx, f] = cpsd(xb(1:T*fs),xb(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(y(1:T*fs),y(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(y(1:T*fs),xb(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
H1b = Gxy./Gxx; 
H2b = Gyy./Gyx;
HTb = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammab = abs(Gxy).^2./(Gxx.*Gyy);
disp('                                  Time for a coffee?');
disp('                                                        S')
disp('                                                       __s__')
disp('                                                      |     |>')
disp('                                                      \_____/')
% case c)
% xc=xb; yc=ya;
[Gxx, f] = cpsd(xb(1:T*fs),xb(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gyy, f] = cpsd(ya(1:T*fs),ya(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gyx, f] = cpsd(ya(1:T*fs),xb(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
H1c = Gxy./Gxx; 
H2c = Gyy./Gyx;
HTc = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammac = abs(Gxy).^2./(Gxx.*Gyy);
clc;disp('Espresso');
% Plot of the results
figure
subplot(2,1,1),plot(f,20*log10(abs(H1a)), f,20*log10(abs(H2a)),f,20*log10(abs(HTa)),f,20*log10(abs(H(1:length(f)))),':')
title('case a: output noise only')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 90 -80 40])
subplot(2,1,2),plot(f,-unwrap(angle(H1a)), f,-unwrap(angle(H2a)),f,-unwrap(angle(HTa)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
axis([0 90 -20 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
subplot(2,1,1),plot(f,20*log10(abs(H1b)), f,20*log10(abs(H2b)),f,20*log10(abs(HTb)),f,20*log10(abs(H(1:length(f)))),':')
title('case b: input noise only')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 90 -80 40])
subplot(2,1,2),plot(f,-unwrap(angle(H1b)), f,-unwrap(angle(H2b)),f,-unwrap(angle(HTb)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
axis([0 90 -20 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
subplot(2,1,1),plot(f,20*log10(abs(H1c)), f,20*log10(abs(H2c)),f,20*log10(abs(HTc)),f,20*log10(abs(H(1:length(f)))),':')
title('case c: both input and output noise')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 90 -80 40])
subplot(2,1,2),plot(f,-unwrap(angle(H1c)), f,-unwrap(angle(H2c)),f,-unwrap(angle(HTc)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
axis([0 90 -20 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
plot(f, Gammaa,f,Gammab,f,Gammac); 
xlabel('Frequency (Hz)'); ylabel('Coherence function')
legend('case a','case b','case c')
axis([0 90 0 1])
clc
% END

Mathieu NOE 2021-1-4

在 MATLAB Online 中打开

hello

I rewrote partially your function so it is more readable for me regarding the governing parameters (NFFT, NOVERLAP, etc...)

adding noises will degrade your estimation, and usually the longer the observation time (T) the better we can reject effects of noise. here we are time limited to 20 s.

I maximise the overlap so that also helps to reduce the noise effect , but to have more averaging I also reduced the NFFT value , so better noise rejection comes to the price of reduced frequency resolution

there is no free lunch in FFT analysis !

data = importdata ('beam_experiment.mat');
x = (data.x); %input
y = (data.y); %output
fs = data.fs;   % sampling frequency
samples = length(x);
T = 1/fs*samples;
NFFT = T*fs/4;
NOVERLAP = round(0.95*NFFT); % maximise overlap to reduce effect of exegeonous input / output noises
% Generate the measuremet noise acting on input and output for case (c)
nx = 0.5*randn(size(x)); % nx=0 for case (a)
ny = 0.5*randn(size(y)); % ny=0 for case (b)
% Case a)
ya = y+ny;
% [Gxx, f] = cpsd(x(1:T*fs),x(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);  %[Pxy,F] = CPSD(X,Y,WINDOW,NOVERLAP,NFFT,Fs)
% [Gyy, f] = cpsd(ya(1:T*fs),ya(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
% [Gxy, f] = cpsd(x(1:T*fs),ya(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
% [Gyx, f] = cpsd(ya(1:T*fs),x(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gxx, f] = cpsd(x,x,hanning(NFFT),NOVERLAP, NFFT, fs);  %[Pxy,F] = CPSD(X,Y,WINDOW,NOVERLAP,NFFT,Fs)
[Gyy, f] = cpsd(ya,ya,hanning(NFFT),NOVERLAP, NFFT, fs);
[Gxy, f] = cpsd(x,ya,hanning(NFFT),NOVERLAP, NFFT, fs);
[Gyx, f] = cpsd(ya,x,hanning(NFFT),NOVERLAP, NFFT, fs);
H1a = Gxy./Gxx; 
H2a = Gyy./Gyx;
HTa = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammaa = abs(Gxy).^2./(Gxx.*Gyy);
disp('                ...still waiting?');
% case b)
xb = x+nx; 
% [Gxx, f] = cpsd(xb(1:T*fs),xb(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
% [Gyy, f] = cpsd(y(1:T*fs),y(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
% [Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
% [Gyx, f] = cpsd(y(1:T*fs),xb(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gxx, f] = cpsd(xb,xb,hanning(NFFT),NOVERLAP, NFFT, fs);  %[Pxy,F] = CPSD(X,Y,WINDOW,NOVERLAP,NFFT,Fs)
[Gyy, f] = cpsd(y,y,hanning(NFFT),NOVERLAP, NFFT, fs);
[Gxy, f] = cpsd(xb,y,hanning(NFFT),NOVERLAP, NFFT, fs);
[Gyx, f] = cpsd(y,xb,hanning(NFFT),NOVERLAP, NFFT, fs);
H1b = Gxy./Gxx; 
H2b = Gyy./Gyx;
HTb = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammab = abs(Gxy).^2./(Gxx.*Gyy);
disp('                                  Time for a coffee?');
disp('                                                        S')
disp('                                                       __s__')
disp('                                                      |     |>')
disp('                                                      \_____/')
% case c)
% xc=xb; yc=ya;
% [Gxx, f] = cpsd(xb(1:T*fs),xb(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
% [Gyy, f] = cpsd(ya(1:T*fs),ya(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
% [Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);  % ? error here : ya instead of y ? 
% [Gyx, f] = cpsd(ya(1:T*fs),xb(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);
[Gxx, f] = cpsd(xb,xb,hanning(NFFT),NOVERLAP, NFFT, fs);  %[Pxy,F] = CPSD(X,Y,WINDOW,NOVERLAP,NFFT,Fs)
[Gyy, f] = cpsd(ya,ya,hanning(NFFT),NOVERLAP, NFFT, fs);
[Gxy, f] = cpsd(xb,ya,hanning(NFFT),NOVERLAP, NFFT, fs);
[Gyx, f] = cpsd(ya,xb,hanning(NFFT),NOVERLAP, NFFT, fs);
H1c = Gxy./Gxx; 
H2c = Gyy./Gyx;
HTc = (Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx);
Gammac = abs(Gxy).^2./(Gxx.*Gyy);
clc;disp('Espresso');
% Plot of the results
figure
subplot(2,1,1),plot(f,20*log10(abs(H1a)), f,20*log10(abs(H2a)),f,20*log10(abs(HTa)),f,20*log10(abs(H(1:length(f)))),':')
title('case a: output noise only')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 90 -80 40])
subplot(2,1,2),plot(f,-unwrap(angle(H1a)), f,-unwrap(angle(H2a)),f,-unwrap(angle(HTa)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
axis([0 90 -20 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
subplot(2,1,1),plot(f,20*log10(abs(H1b)), f,20*log10(abs(H2b)),f,20*log10(abs(HTb)),f,20*log10(abs(H(1:length(f)))),':')
title('case b: input noise only')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 90 -80 40])
subplot(2,1,2),plot(f,-unwrap(angle(H1b)), f,-unwrap(angle(H2b)),f,-unwrap(angle(HTb)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
axis([0 90 -20 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
subplot(2,1,1),plot(f,20*log10(abs(H1c)), f,20*log10(abs(H2c)),f,20*log10(abs(HTc)),f,20*log10(abs(H(1:length(f)))),':')
title('case c: both input and output noise')
legend('|\itH\rm_1|','|\itH\rm_2|','|\itH\rm_T|','|\itH\rm|')
xlabel('Frequency (Hz)'); ylabel('Modulus (dB)')
axis([0 90 -80 40])
subplot(2,1,2),plot(f,-unwrap(angle(H1c)), f,-unwrap(angle(H2c)),f,-unwrap(angle(HTc)),f,unwrap(angle(H(1:length(f)))),':')
xlabel('Frequency (Hz)'); ylabel('Phase (rad)')
axis([0 90 -20 0.5])
legend('\angle H_1','\angle H_2','\angle H_T','\angle H')
%
figure
plot(f, Gammaa,f,Gammab,f,Gammac); 
xlabel('Frequency (Hz)'); ylabel('Coherence function')
legend('case a','case b','case c')
axis([0 90 0 1])
clc
% END

Mathieu NOE 2021-1-4

在 MATLAB Online 中打开

see

% [Gxy, f] = cpsd(xb(1:T*fs),y(1:T*fs), hanning(T2*fs),T2*fs/2, T1*fs, fs);  % ? error here : ya instead of y ? 

I corrected it

Christina Reid 2021-1-4

aaaah! Bless your soul! I have been banging my head aganist the wall for this one!

请先登录，再进行评论。

Evaluate the system impose response function and the FRF with given input and output data

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

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Evaluate the system impose response function and the FRF with given input and output data

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

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

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