FFT plotting trouble: Concatenated inputs

I am currently working on tap testing and am running into an issues with my FFT now that I am trying to Concatenate them to improve the "resolution". The plots max values fit what I expect however i need it to be a line and not have all the points underneath. I'm uncertian if my issue is in my measurements or in my code.

clear all
close all
clc
% each channel is notated with _n where n is the channel
% Hammer is on channel 1
% Accelerometer is on channel 2

% dataset_plot(11,'y')
% dataset_plot(12,'y')
% dataset_plot(13,'y')
%
% dataset_plot(21,'y')
dataset_plot(22,'y')

loading 22_y_1_1.csv
loading 22_y_2_1.csv
loading 22_y_3_1.csv
loading 22_y_4_1.csv
loading 22_y_5_1.csv

Warning: Integer operands are required for colon operator when used as index.

Warning: Integer operands are required for colon operator when used as index.

The goal is to have a single clear line so that I can create the transfer function to describe the system.

function dataset_plot(n,direction)
%% data Loading
out_2 = [];
out_1 = [];
for f_n = 1:5
    fprintf('loading '+ string(n) + '_' + direction + '_' + string(f_n) + '_1.csv'+'\n')
    raw_1 = table2array(readtable(string(n)+'_' + direction + '_' + string(f_n) + '_1.csv'));% read in hammer
    raw_cut_1 = raw_1(3:height(raw_1),:); %remove header

    if f_n == 1 %if its the first loop find signal length for window
        window = hanning(length(raw_cut_1));
        sample_rate  = mean(diff(raw_cut_1(:,1))); %time between samples
        Fs = 1/sample_rate; %Frequency in HZ
    end

    % Zero correction
    Hammer_Zeroing = raw_cut_1(1:(length(raw_cut_1(:,2))*0.2)-1);%fix the offset
    Hammer_adj = mean(Hammer_Zeroing);
    raw_cut_1(:,2) = raw_cut_1(:,2) - (Hammer_adj);
    %concatenate

    out_1 = [out_1; raw_cut_1(:,2).*hanning(length(raw_cut_1))]; 
    
    raw_2 = table2array(readtable(string(n)+'_' + direction + '_' + string(f_n) + '_2.csv'));% read in accel
    raw_cut_2 = raw_2(3:height(raw_1),:); %remove header
    
    %Zero correction
    Accel_Zeroing = raw_cut_2(1:(length(raw_cut_2(:,2))*0.2)-1);%fix the offset
    Accel_adj = mean(Accel_Zeroing);
    raw_cut_2(:,2) = raw_cut_2(:,2) - (Accel_adj);
    
    %concatenate
    out_2 = [out_2; raw_cut_2(:,2).*hanning(length(raw_cut_2))];
    clear raw_cut_2 raw_cut_1 raw_1 raw_2 padding Hammer_adj Hammer_Zeroing Accel_adj Accel_Zeroing
end
%% create X axis
Hammer = out_1;%(:,2);
Accel = out_2;%(:,2);
x_axis = linspace(0,sample_rate*length(Accel),length(Accel));
Accel = vertcat(Accel,ones(length(x_axis)-length(Accel),1).*mean(diff(Accel)));
meas_time_2 = sample_rate*length(x_axis);

%% run FFT functions
[H_f,H_P1] = FFT_ploting(smoothdata(Hammer,"gaussian",20),meas_time_2);
[A_f,A_P1] = FFT_ploting(smoothdata(Accel,"gaussian",20),meas_time_2);
%% fix units
H_P1 = H_P1./0.00225; %convert from V to N @ 2.25 mV/N
A_P1 = A_P1./0.0102; %convert from V to m/s^2 10.2 mV/(m/s²)

%% plot
close all

figure('Name','raw data check')
hold on
plot(Accel)
plot(smoothdata(Accel,"gaussian",10))

figure('Name','FFT Plots')
hold on
yyaxis left
plot(H_f,H_P1);
ylabel 'lbf'
yyaxis right
plot(A_f,A_P1);
ylabel 'm/s^2'
grid on
xlim([5 3000]); %accelerometer limit of 3000Hz
%we are focusing on 5 to 500 Hz
title('Accelerometer & Hammer FFT')
legend('Hammer FFT','Acceleromter FFT')
xlabel('Hz');
end

function [f,P1] = FFT_ploting(fun_data,meas_time)
% count = count + 1; %FFT itteration count
Fs = (meas_time/length(fun_data))^-1; %frequecy in Hz
T = 1/Fs; %meas time between samples (inverse)
L = length(fun_data);
t = (0:L-1)*T;
f = Fs*(0:(L-1)/2)/L;
fft_data = fft(fun_data,L);
P2 = (abs(fft_data/L));
P1 = P2(1:(L+1)/2);
P1(2:end) = 2*P1(2:end);
end

Answer 1

[Mathieu NOE]

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