hello
when you do modal analysis / FRF measurements with a hammer , you should split the record into individual hit input and output signals , then do the averaging of the FRFs (and coherence)
I may not directly give the code for your data , but the process is described here , based on the simulated response of a mass / spring / dash pot system. You can see the measured FRF pretty much overlays with the theoretical Bode plot of the system
%% Modal Analysis / hammer test / FRF acquisition and FRF computation tool
%% FFT processing parameters
Fs = 1024*4;
nfft = 1024*2;
channels = 1;
% dummy hammer test data - or real data if you have
Nhits = 8;
duration_hits = 1; % mean time duration between hits (in seconds)
% Mass / spring / dash pot model
M = 25; % kgs
fn = 100; % Hz (natural frequency)
K = M*(2*pi*fn)^2;
C = 0.05*sqrt(K*M);
% analog model acceleration vs force
num = [1 0 0];
den = [M C K];
% discretization
dt = 1/Fs;
samples_one_hit = duration_hits*Fs;
[numd,dend] = c2dm(num,den,1/Fs,'tustin');
% create some dummy hammer gaussian pulses with variable force amplitude
f0 = Fs/2.5; % pulse frequency
Var_Sq = 1; % pulse variance (squared)
t = dt*(0:samples_one_hit-1)';
t_pulse = max(t)/8;
offset = 0.5*duration_hits*(rand(Nhits,1)); % random time shifts
amplitude = 10*(1 + 1*(rand(Nhits,1))); % random amplitude
force_signal = [];
for k = 1:Nhits
pulse_signal = amplitude(k)*exp(-(((t-t_pulse-offset(k)).^2).*(f0^2))./Var_Sq);
force_signal = [force_signal;pulse_signal];
end
% now all samples are concatenated to make one signal with random delta t
% between pulses
samples = numel(force_signal);
time = dt*(0:samples-1)';
% apply model to force signal (input) to create acceleration signal (output)
accel_signal = filter(numd,dend,force_signal);
accel_signal = accel_signal + 0.001*randn(size(accel_signal));
% now let's trigger the data (using "zero crossing" on the force signal)
% this is basically the same way a FFT/network analyser operates
pre_trig_duration = 10e-3; % add some time to shift the pulse signal to the right
force_signal_level = 1;
[Zx] = find_zc(time,force_signal,force_signal_level);
Zx = Zx - pre_trig_duration;
%recommended max nfft
recommended_max_nfft = floor(min(Fs*diff(Zx)));
%% check valid segments (duration in samples must be above nfft
% pulses durations in samples :
pd = [Fs*diff(Zx); samples-Fs*Zx(end)];% add the last pulse duration :
data_valid = find(pd>=nfft);
figure(2),
subplot(2,1,1),plot(time,force_signal)
% warning message
if nfft>recommended_max_nfft
text(max(time)*0.1,max(force_signal)*0.9,['Selected nfft = ' num2str(nfft) ' is above recommended max value : ' num2str(recommended_max_nfft)]);
end
subplot(2,1,2),plot(time,accel_signal)
%% FFT processing
% windowing
windowx = ones(nfft,1); % no window on force signal
windowy = ones(nfft,1); % no window on response sensor signal
% windowy = exp(-3*(0:nfft-1)/nfft)'; % 0.1 exp window on response sensor signal
% initialisation
Pxx = zeros(nfft,1);
Pxy = zeros(nfft,channels);
Pyy = zeros(nfft,channels);
% check if data_valid is not empty (can happen if you choose nfft grossly
% oversized)
if isempty(data_valid)
error('No data or nfft too large - check data & nfft settings');
else % we have data and nfft is correctly defined
for ck = 1:numel(data_valid)
% do it only for valid segments
k = data_valid(ck);
% select time data
ind_start = find(time>=Zx(k),1,'first');
% check valid segments (duration in samples must be above nfft
ind = (ind_start :ind_start + nfft -1);
time_measure = time(ind);
time_measure = time_measure - time_measure(1);
force_signal_measure = force_signal(ind);
accel_signal_measure = accel_signal(ind);
figure(3), % just to show the individual gaussian pulses with different time offsets
subplot(2,1,1),plot(time_measure,force_signal_measure)
hold on
subplot(2,1,2),plot(time_measure,accel_signal_measure)
hold on
leg_str1{ck} = ['hit # ' num2str(k)];
leg_str2{ck} = ['record # ' num2str(k)];
% FFT processing
x = force_signal_measure;
y = accel_signal_measure;
xw = windowx.*x;
yw = (windowy*ones(1,channels)).*y;
Xx = fft(xw,nfft);
Yy = fft(yw,nfft,1);
Xx2 = abs(Xx).^2;
Xy2 = Yy.*(conj(Xx)*ones(1,channels));
Yy2 = abs(Yy).^2;
Pxx = Pxx + Xx2;
Pxy = Pxy + Xy2;
Pyy = Pyy + Yy2;
end
subplot(2,1,1),legend(leg_str1); % add legend of valid data
subplot(2,1,2),legend(leg_str2); % add legend of valid data
% Select first half
if ~any(any(imag([x y])~=0)) % if x and y are not complex
if rem(nfft,2) % nfft odd
select = [1:(nfft+1)/2];
else
select = [1:nfft/2+1]; % include DC AND Nyquist
end
Pxx = Pxx(select);
Pxy = Pxy(select,:);
Pyy = Pyy(select,:);
else
select = 1:nfft;
end
% set up output parameters
Txy = Pxy ./ (Pxx*ones(1,channels)); % transfer function estimate
Cxy = (abs(Pxy).^2)./((Pxx*ones(1,channels)).*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;
end
%% plot model vs identified FRF
flow = 10;
fhigh = Fs/2.56;
freq = logspace(log10(flow),log10(fhigh),300);
[gm,pm] = bode(num,den,2*pi*freq);
figure(4)
subplot(3,1,1),semilogx(freq,20*log10(gm),'b',f,20*log10(abs(Txy)),'r');
title(['FRF comparison / Valid records : ' num2str(data_valid')])
xlim([flow fhigh]);
ylabel('Magnitude (dB)');
legend('input model','measured');
subplot(3,1,2),semilogx(freq,pm,'b',f,180/pi*(angle(Txy)),'r');
xlim([flow fhigh]);
ylabel('Phase (°)');
legend('input model','measured');
subplot(3,1,3),semilogx(freq,ones(size(freq)),'b',f,Cxy,'r');
xlim([flow fhigh]);
xlabel('Frequency (Hz)');
ylabel('Coh');
legend('input model','measured');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)
% find time values of positive slope zero (or "threshold" value ) y
% crossing points
% put data in rows
x = x(:);
y = y(:);
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end
