Frequency Response Function and FFT for Modal Analysis

I am experimentally testing a material to retrieve its natural frequencies through modal analysis. The type of testing is based around impulse response. My sampling frequency was 10kHz. My acccelerometer's sensisitivty was 100mV/g and my impact hammer's sensitivity was 2.25mV/N. My aim is to obtain an accurate FRF graph. My entire code:

%reading excel file

dataset = xlsread('BrassFRF.xlsx','Sheet1','A1:C50000');

%%Creating and filling data variables

a=dataset(:,2); 
n=dataset(:,3);

%%Converting acceleraltion voltages into acceleration (mm/sec^2)

a=a.*((1/0.1)*9.81);

%%Converting force voltages into Newton (N)

n=n.*(1/0.00225);

%%Creating FFT using acceleration

L=length(a);
NFFT=2.^nextpow2(L);
V=fft(a,NFFT)/L;
Fs=10000;
freq=Fs/2*linspace(0,1,NFFT/2+1);
subplot(3,1,1)
plot(freq,2*abs(V(1:NFFT/2+1)))
xlabel('Frequency (Hz)');
title('Velocity FFT ')

%%Creating FFT for using force

F=fft(n,NFFT)/L;
subplot(3,1,2)
plot(freq,2*abs(F(1:NFFT/2+1)))
xlabel('Frequency (Hz)');
title('Impulse FFT ')

%%FRF

FRF=V./F;
subplot(3,1,3)
plot(abs(FRF))
xlabel('Frequency (Hz)');
title('FRF')

Please tell me if you find any errors that can be the cause of an incorrect FRF. My suspect is the impulse FFT.

[Lingyuan Kong]

The last section is wrong, it should be:

FRF=V./F;

subplot(3,1,3)

plot(freq, abs(FRF(1:NFFT/2+1)))

xlabel('Frequency (Hz)');

title('FRF')

Answer 1

[durukan dilek]

did you obtain frf of acceleration/force? why did you write V/F? why did you use velocity?

L=length(a); NFFT=2.^nextpow2(L); V=fft(a,NFFT)/L; Fs=10000; freq=Fs/2*linspace(0,1,NFFT/2+1); subplot(3,1,1) plot(freq,2*abs(V(1:NFFT/2+1))) xlabel('Frequency (Hz)'); title('Velocity FFT ') set(gca, 'YScale','log')

Answer 2

[Deniz Kny]

how did you know your sampling frequency was 10khz? I have a accelerometer data and, trying to obtain fft like you. i need a Fs but i dont know how can i get this info.

[imthiyas manarikkal]

You can actually check your resonance frequency or natural frequency. The sampling frequency must be twice of your natural frequency (Nyquiste Theorem)

Answer 3

[Jose Ordaz]

Hello,

i'm working in the same way, trying to obtain the natural frequency and damping of some materiales using the information from a hammer and a accelerometer. Can you please help me with this?, how you solve your problem with the FRF?.

thanks

Answer 4

[Mathieu NOE]

hello

yes I a coming very late in the show but maybe there is still a need for a better code

you can improve the FRF estimation if you have multiple hits in the same record - with enough time in between hits to let the response decay to zero - see my other answer here : How can I obtain accurate and high-quality FRF results? - MATLAB Answers - MATLAB Central

results :

T = Mode Frequency Damping

____ _________ _________

1 1845.7 0.0072161

2 2336.8 0.0060413

3 2895.7 0.0080089

4 3515.7 0.0070846

code :

%% Load data
dataset = xlsread('BrassFRF.xls','Sheet1','A1:C50000');
%%Creating and filling data variables
accel_signal=dataset(:,2); 
force_signal=dataset(:,3);

%%Converting acceleraltion voltages into acceleration (mm/sec^2)
accel_signal=accel_signal.*((1/0.1)*9.81);
accel_unit = 'mm/sec^2';

%%Converting force voltages into Newton (N)
force_signal=force_signal.*(1/0.00225);
force_unit = 'N';

%% Extract force and acceleration data
time = dataset(:, 1);
dt = mean(diff(time));
Fs = 1/dt;


samples = numel(force_signal);
channels = size(accel_signal, 2);
time = (0:samples-1)*dt; % we need to re-create the time signal 

%% FFT processing parameters
nfft = 1024*2;

% Trigger the data (using "zero crossing" on the force signal) 
pre_trig_duration = 10e-3; % (seconds) - add some time to shift the pulse signal to the right
force_signal_level = 0.5*max(force_signal);

[Zx] = find_zc(time, force_signal, force_signal_level);
Zx = Zx - pre_trig_duration;

% Recommended max nfft
if numel(Zx) > 1
    recommended_max_nfft = floor(min(Fs * diff(Zx)));
else
    recommended_max_nfft = numel(time);
end

%% Check valid segments (duration in samples must be above nfft)
pd = [Fs * diff(Zx); samples - Fs * Zx(end)]; % Pulse durations in samples
data_valid = find(pd >= nfft);

figure(1),
subplot(2,1,1), plot(time, force_signal)
title('Force Signal (raw)');
xlabel('Time (s)');
ylabel(force_unit);

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)
title('Acceleration Signal (raw)');
xlabel('Time (s)');
ylabel(accel_unit);

%% FFT Processing - Combined Section
windowx = ones(nfft, 1); % No window on force signal (default)
% windowx(round(0.2*nfft:end)) = 0 ; % "force" window on force signal ( window = 1 during first 20% of time, then 0)
% windowx(round(0.1*nfft:end)) = 0 ; % "force" window on force signal ( window = 1 during first 10% of time, then 0)

windowy = ones(nfft, 1); % No window on acceleration signal
% windowy = exp(-2.3*(0:nfft-1)/nfft)'; % 0.1 exp window on response sensor signal 
% windowy = exp(-4.6*(0:nfft-1)/nfft)'; % 0.01 exp window on response sensor signal 

Pxx = zeros(nfft, 1);
Pxy = zeros(nfft, channels);
Pyy = zeros(nfft, channels);

% Check if data_valid is not empty
if isempty(data_valid)
    error('No data or nfft too large - check data & nfft settings');
else
    for ck = 1:numel(data_valid)
        k = data_valid(ck);
        % Select time data 
        ind_start = find(time >= Zx(k), 1, 'first');
        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(2), % just to show the individual gaussian pulses with different time offsets 
        subplot(2,1,1),plot(time_measure, windowx.*force_signal_measure)
        xlabel('Time (s)');
        ylabel(force_unit);
        hold on 
        subplot(2,1,2),plot(time_measure,(windowy*ones(1,channels)).*accel_signal_measure)
        xlabel('Time (s)');
        ylabel(accel_unit);
        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 FRF and Coherence
flow = 10;
fhigh = Fs / 2.56;


%% Modal Analysis (Natural Frequencies and Damping Ratios)
N = 4; % Number of dominant modes to identify
% Only use valid data (coherence above 0.9)
ind = Cxy > 0.9;

FR = [flow fhigh];
[fn, dr] = modalfit(Txy(ind), f(ind), Fs, N, 'FitMethod', 'pp', 'FreqRange', FR);
T = table((1:N)', fn, dr, 'VariableNames', {'Mode', 'Frequency', 'Damping'});

%%
figure(3)
subplot(3, 1, 1), semilogy(f, abs(Txy), 'r',fn,interp1(f,abs(Txy),fn),'db');

title(['FRF Estimation based on Valid Data Chunks']);
xlim([flow fhigh]);
ylabel(['Magnitude (' accel_unit ' / ' force_unit ')']);

subplot(3, 1, 2), plot(f, 180/pi * (angle(Txy)), 'r');
xlim([flow fhigh]);
ylabel('Phase (°)');

subplot(3, 1, 3), plot(f, Cxy, 'r');
xlim([flow fhigh]);
xlabel('Frequency (Hz)');
ylabel('Coherence');
ylim([0 1.1]);



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