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

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?

Answer 1

[Mathieu NOE]

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')

[Christina Reid]

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]

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]

ahh! Thank you! I was trying different things and wanted I think I wanted to see how the system would respond from 0-25 Hz and I never changed it. I am going to see if I can incorporate your code and my code to produce the same graphs with the given input/output signal example

[Christina Reid]

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

[Mathieu NOE]

hello Christina

ilaplace works on symbolic transfer functions, not numerical data (is what tfestimate generates)

[Christina Reid]

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]

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]

also , I believe there was a minor error in case c)

[Mathieu NOE]

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]

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