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:

enter image description 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()

enter image description here

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 
NB=8;  % poles
NA=3;  % zeros
% Initial Parameters declaration 
% ---------------------------------------------------------
% 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
% stem(Sw{Pos});
grid on;
title(['Impulse response']);
xlabel('Sample number');
 xlim([0 20000])

%% #2 Pressure registered at the sensor introducing white noise excitation 
xlabel('time [s]'); 
title('Pressure registered at the sensor position 7');
axis tight
grid on

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

function [fs,t,wn,kn,xa,xs,etha,f_hat,m,M,N, 
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
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
whitenoise = wgn(Nsamples,1,0);             % white noise generation 

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
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. 
%% 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

if vall==1;
figure()           % Z-plane zero-pole plot.

figure()           % Bode frequency response of dynamic systems

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.