Frequency Response Function and FFT for Modal Analysis

68 次查看（过去 30 天）
R SW
R SW 2016-12-12
回答： Mathieu NOE 2024-12-10，11:00
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.
  1 个评论
Lingyuan Kong
Lingyuan Kong 2021-1-22
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')

请先登录，再进行评论。

请先登录，再回答此问题。

回答（4 个）

durukan dilek
durukan dilek 2017-11-1
编辑：durukan dilek 2017-11-1
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')
Deniz Kny
Deniz Kny 2019-3-6
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.
  1 个评论
imthiyas manarikkal
imthiyas manarikkal 2019-11-8
You can actually check your resonance frequency or natural frequency. The sampling frequency must be twice of your natural frequency (Nyquiste Theorem)

请先登录，再进行评论。

Jose Ordaz
Jose Ordaz 2019-12-19
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
Mathieu NOE
Mathieu NOE 2024-12-10，11:00
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

类别

在 Help Center 和 File Exchange 中查找有关 Multirate Signal Processing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by