Stable estimation of the frequency response data using Matlab command invfreqz?

I had a hard time 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?

 % %Vibration of simply supported beam
clc; clear all; close all;
% Parameter to the function invfreqz, which can be found in function 'pressuresignal_beam'
NB=12;  % poles
NA=32;  % 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});
grid on;
title(['Impulse response']);
xlabel('Sample number');
ylabel('Amplitude');
xlim([0 1000])
%%#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 [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, 2*pi(), length(womega));          % contains the normalized frequency values within the interval [0,2 Pi]
clear A B                                             % used in invfreqz  
for i=1:M
[B{i},A{i}]=invfreqz(w_disp{i},wnorm,'complex',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

% 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