How to make an FRF with my accelerometer Data and FFT?

Dear, Matlab Users.

I have some inquiries about my accelerometer Data, that was processed by a Dewesoft Software.

So, I am gonna share my code and data, and I would like to know if there any suggestion to this resolution.

Sensitivy: 102 mV/g, or 10.40 mV/m/s^2

%Load data
clear all
load("C:\Users\joant\Desktop\Field_Research\Matlab_Code\data_test.mat")
whos
close all

%Plot the array according to the data
%Data of Channel 1
g1= transpose(Data1_352C33OLD)
t1= transpose(Data1_time_352C33OLD); % time
dt=1/(Sample_rate) % period
L=length(g1)
%Iterate
% Figure of Channel 1
figure(1)
plot(t1,g1,'b.-')
set(gca,'Fontsize',14)
xlabel('time(s)')
ylabel('g')
xlim([min(t1) max(t1)]) 
%ylim([-6 6])
grid on

% Calculate Fourier transform of force
points=round((t1(L))/dt)  %Calculate the points
A_F= g1; % This is the Force for setting in the fft.
freq= (0:1/(dt*points):(1-1/(2*points))/dt)'; %Hz
freq=freq(1:round(points/2+1),:); %Hz
length(freq)
length(fft(A_F))

%In this part I need help
A_fft=fft(g1); %without the conversion by 9.8 m/s^2
abs_A=abs(A_fft)
figure(4)
plot(freq,abs_A(1:length(freq)))
xlim([0 6]) % Because 5 Hz is 300 RPM
%ylim([-400 400])
xlabel('Frequency (Hz)')
ylabel('g')

%According to the practice in the laboratory the result is the following
%image.

This was my process to obtain the FFT, I would like to know if this process is correct?. Because when I plot this data is not the same values in y-axis that I had in the Dewesoft software.

Dewesoft fft(g) (attached) vs matlab fft(g) (attached)

Question #2

According to a Mechanical vibration book, I have the next equation to find the FRF. I applied this but my FRF y axis value is too high so, I don't know If I didn't applied correctly. That information is in my code and I use the following equation.

w= natural frequency = 31.42 rad/s^2.

I assumpted that A/F is g because according to the book, it says the next information: "Because the acceloremeter produces a voltage that is proportional to acceleration, we obtain the accelerance, A/F(w), from the dynamic signal analyzer".

My graph result:

This graph it doesn't have sense, because I need values aroun of 0.15-0.17 in the first vibration mode. So, I would like to know if there any suggestion how to resolve this.

Best Regards,

Jo98.

Answer 1

[Mathieu NOE]

hello

I got the correct result if I assume that your data is in g and not volts (so your analyser has already used the sensor sensivity when storing the data)

there can be some minor difference in terms of frequency and amplitude as I am not aware of what type of window / nfft is used on your side

load("data_test.mat")
%Plot the array according to the data
%Data of Channel 1
signal = double(Data1_352C33OLD);
time= Data1_time_352C33OLD; % time

[samples,channels] = size(signal);
dt = mean(diff(time));
Fs = 1/dt;

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 10;
if decim>1
    Fs = Fs/decim;
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    end
   signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000;    % for maximal resolution take NFFT = signal length (but then no averaging possible)
OVERLAP = 0.75; % overlap will be used only if NFFT is shorter than signal length


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
plot(time,signal,'b');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);

figure(2),plot(freq,spectrum_raw,'b');grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude');


function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

[Mathieu NOE]

regarding your second question, accelerance is the ratio of acceleration divided by force , but I don't see any force signal in your data

also, when we want a good FRF estimate from time signals , they must contain enough energy accross the entire spectrum , otherwise the result (coherence) will be poor and inaccurate.

so either we use impulse or stationnary signals (random, chirp,...) but not a single frequency as it seems to be the case in your data. Unles you want the FRF computation valid only at this particular frequency

[Mathieu NOE]

Dynamic Stiffness, Compliance, Mobility, and more... (siemens.com)

[Jo98]

Hello, thanks so much for your answer,it was really helpful. Also I'm gonna try if with this result I can plot the FRF.

So, Could I ask some questions about the code by DM or email?

Best Regards.

[Mathieu NOE]

yes I am available to further help you either on this forum or per email - as you prefer

but also , if my answer here has helped you , are you keen to accept it ? tx

[Jo98]

Yes, I'm gonna contact you by DM, asking you the email.

Regard with your FRF, I replace this vector with the formula that I post in the original post, and it works. So A_F is (spectrum_raw*9.8 m/s^2). This because the conversion between g and m/s^2.

I assure that A_F is g because the theory of the book that I learnt about FRF says:" Because the acceloremeter produces a voltage that is proportional to acceleration, we obtain the accelerance, A/F(w), from the dynamic signal analyzer".

Best Regards.

[Mathieu NOE]

tx for accepting my answer

NB that what we display above is only the autospectrum of the acceleration signal

for the accelerance FRF we need a force signal , as already mentionned above, it's not in your data file.

Also why do you want an accelerance FRF at only one frequency ?