# Estimation of the frequency response data using Matlab command invfreqz?

I would like to determine frequency response and then impulse response of the displacement equation (eq. 1 please see screen shot of the task below). In this example we study a response of the finite beam to a harmonic force at certain position.

The description of my work is summarized here:

References 1 L. Meirovitch, Fundamentals of Vibrations. McGraw-Hill, Boston, 2001. 2 C.R. Fuller, S.J. Elliott and P.A. Nelson, Active Control of Vibration. Academic Press, London, 1996.

Consider a simply supported beam with parameters given in MATLAB function inputparameter_beam()

Here is my approach: The first step is to compute frequency response of a vibrating beam. The total response of the beam (see eq. 1) includes the harmonic time component. However to compute the FR I have omitted this component. Subsequently, the summation I have performed in terms of changed angular frequency ‘omega’. For this computation please refer to function pressuresignal_beam(). The vector of angular frequency ‘womega’ is chosen arbitrarily, my choice is based on half of sampling frequency ½ fs=2000/2 =1000 Hz. The impulse response ‘Sw{i}’ calculated using numerator ‘B{i}’ and denominator ‘A{i}’ of invfreqz matlab command. To obtain Stable estimation of the system I have put large degree of polynomial in nominator as well as denominator. This matlab command invfreqz is a “discrete filter least squares fit to frequency response data”.

The matlab script ‘SysID_vibrating_beam’ is divided into separate function for more clarification. The function are: inputparameter_beam; pressuresignal_beam; You can find useful plots in “SysID_vibrating_beam” which are currently comment out: Impulse response plot; Pressure registered at the sensor; Frequency response function plot. By compiling whole script automatically plot show up: z-plane poles and zeros and bode plot;

Any ideas how to obtain stable/proper estimation of frequency response data and then an impulse responses?

%% SysID_vibrating_beam- Vibration of simply supported beam
clc; clear all; close all;

% Parameter to the function invfreqz, which can be found in function
'pressuresignal_beam'
NB=8;  % poles
NA=3;  % zeros
womega=0:1000;
%---------------------------------------------------------
% Initial Parameters declaration
[fs,t,wn,kn,xa,xs,etha,f_hat,m,M,N,
whitenoise,Nsamples,Pos]=imputparameter_beam();
% ---------------------------------------------------------
% Signals p on N microphone positions at the constant spacing
[Sw, p,w_disp,hh,freq] = pressuresignal_beam(whitenoise,M,N,kn,xa,xs,etha,wn,f_hat,m,fs,Pos,1,NB,NA,womega);

%% #1 Impulse response plot
figure()
% stem(Sw{Pos});
plot(Sw{Pos})
grid on;
title(['Impulse response']);
xlabel('Sample number');
ylabel('Amplitude');
xlim([0 20000])

%% #2 Pressure registered at the sensor introducing white noise excitation
figure()
plot(t,p{Pos})
xlabel('time [s]');
ylabel('Signal');
title('Pressure registered at the sensor position 7');
axis tight
grid on

%% #3 frequency response function of a simply supported beam
figure()
plot(womega, w_disp{Pos})
xlabel('Frequency [Hz]');
ylabel('Displacement [m]');
title('Frequency response function');
axis tight
grid on
set(gca,'XMinorGrid','on')

function [fs,t,wn,kn,xa,xs,etha,f_hat,m,M,N,
whitenoise,Nsamples,Pos]=imputparameter_beam();
Nsamples = 2e4;         % number of samples for disturbance signal d
fs = 2000;              % sampling frequency in Hz
t = (0:Nsamples-1)/fs;  % time vector
l=0.38;                 % beam length in m
h=0.002;                % beam height in m
b=0.04;                 % beam width in m
E=2.1e11;               % elasticity modulus in Pa
rho=8800;               % density of material in kg/m3
xa=0.038;               % position of complex force f in m
etha=0.03;              % loss factor
f_hat=1;                % complex amlitude in N
m=rho*l*b*h;            % mass
m_prime=rho*h;          % mass per unit area
I=b*h^3/12;             % moment of inertia
Pos=7;                  % choice of sensor position to observe
N=30;                   % summation range
M=20;                   % number of sensors distributed along beam
%---------------------------------------------------------
for j=1:M-1             % calculate sensor position
xs(j+1)=j/(M-1)*l;
end
%%---------------------------------------------------------
j=0;
for a=1:N
kn(a)=a*pi()/l;                            % characteristic wavenumbers
wn(a)=kn(a).^2*sqrt((E*I)/(b*m_prime));    % angular resonance frequencies
end
whitenoise = wgn(Nsamples,1,0);             % white noise generation
end

function [Sw, p,w_disp,hh,freq] = pressuresignal_beam(whitenoise,M,N,kn,xa,xs,etha,wn,f_hat,m,fs,Pos,vall,NB,NA,womega);
%% The total response of the beam without the harmonic time component 'exp(-1j*wom*t)'
womega=womega';     % vector declaration of angular frequency in rad/s (selected freely/predicted range)
% used to calculate displacement w_disp
for i=1:M
w=0;
for a=1:N
w=w + (  sin(kn(a)*xa) *sin(kn(a)*xs(i))  )./(         womega.^2+2*1i*etha*womega*wn(a)-wn(a)^2 );    % Eq 1.
end
w_disp{i}=-f_hat/m*w;
end
%% Signals p on N microphone positions at the constant spacing
wnorm = linspace(0, pi(), length(w_disp{1}));          % contains the     normalized frequency values within the interval [0,Pi]
clear A B                                             % used in invfreqz
for i=1:M
[B{i},A{i}]=invfreqs(w_disp{i},wnorm,NB,NA);      % Discrete filter least     squares fit to frequency response data
[Sw{i},ts{i}]=impz(B{i},A{i},20001,fs);                     % Impulse response from numerator and denominator
[hh{i},freq{i}] = freqz(B{i},A{i},'whole',fs);              % Frequency     response of digital filter
p{i} = filter(Sw{i},1,whitenoise);                          % Expected pressure on the sensore
end

if vall==1;
figure()           % Z-plane zero-pole plot.
zplane(B{Pos},A{Pos})

figure()           % Bode frequency response of dynamic systems
bode(B{Pos},A{Pos})
else
end
end