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 ?

[oybozkurt]

If we have both acceleration (in g) and force (Newtons) acquired from accelerometer and impact hammer from a cantilever beam, how we can obtain FRF to evaluate natural frequency and damping ratio?

[Mathieu NOE]

hello @oybozkurt

i tried to reproduce the way FFT analysers works for hammer testing. The code below assume a simple mass / spring / dashpot system excited by hammer pulses .

As I don't have always experimental data to use for demo, I prefered here to generate some dummy data , but you can easily adpat the code to your own needs .

I used modalfit (from the Signal Processing Tbx) to give the damping ratio and eigenfrequency from the estimated FRF

hope it helps !

%% Modal Analysis ( hammer test ) - FRF computation 

%% 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)); % add some noise as in real experimental life 

% 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 (or use boxcar)
windowy = ones(nfft,1); % no window on response sensor signal (or use a exp window as suggested below if the signal decay rate is too slow).
% 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');

%% natural frequencies & damping ratioes
N = 1 ; % number of (dominant) modes to identify
% only use valide data (coherence above 0.9)
ind = Cxy>0.9;
[fn,dr] = modalfit(Txy(ind),f(ind),Fs,N,'FitMethod','pp');
T = table((1:N)',fn,dr,'VariableNames',{'Mode','Frequency','Damping'})



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

    

[oybozkurt]

Thank you for your quick reply. For example, the acceleration and force values obtained as a result of a hit with a hammer are available in the excel file. Can you check this sample data?

[Mathieu NOE]

hello

so I adapted my code for your experimental data - so there is only one hit in your record, I wish you could redo the test with multiple hits in the same record (make 5/6 hits with approx 5 seconds interval).

the benefit of making multiple hits is to reduce the noise effects and average your response curve.

nevertheless , this is the result obtained so far with your single hit record :

hope it helps !

Results

Mode Frequency Damping

____ _________ _________

1 66.147 0.0081933

2 307.36 0.031123

Code :

%% Modal Analysis ( hammer test ) - FRF computation 
%% load data
filename = "Acceleration-Force.xlsx";
headerlines = 10;
% force data
data = readmatrix(filename,"Sheet","Force","NumHeaderLines",headerlines,"ExpectedNumVariables",2);
timef = data(:,1);
force_signal = data(:,2);
force_unit = 'N';
% accel data
data = readmatrix(filename,"Sheet","Acceleration","NumHeaderLines",headerlines,"ExpectedNumVariables",2);
timea = data(:,1);
accel_signal = data(:,2);
accel_unit = 'g';

if sum(abs(timef - timea))>eps
    error('Time vectors don''t match - double check your data')
else
    time = timef;
    clear timef timea
end

dt = mean(diff(time));
Fs = 1/dt;
samples = numel(force_signal);
channels = size(accel_signal,2);

%% FFT processing parameters
nfft = 1024*4;

% 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; % (seconds) - 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
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
% 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)
title('force signal');
xlabel('Time (s)');
ylabel(force_unit);

% 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)
title('accel signal');
xlabel('Time (s)');
ylabel(accel_unit);

%% FFT processing 
% windowing
windowx = ones(nfft,1); % no window on force signal (or use boxcar)
windowy = ones(nfft,1); % no window on response sensor signal (or use a exp window as suggested below if the signal decay rate is too slow).
% 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);

            % 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

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

figure(4)
subplot(3,1,1),loglog(f,abs(Txy),'r');
title(['FRF estimation based on valid data chuncks # : ' num2str(data_valid')])
xlim([flow fhigh]);
ylabel(['Magnitude (' accel_unit ' / ' force_unit ')']);
subplot(3,1,2),semilogx(f,180/pi*(angle(Txy)),'r');
xlim([flow fhigh]);
ylabel('Phase (°)');
subplot(3,1,3),semilogx(f,Cxy,'r');
xlim([flow fhigh]);
xlabel('Frequency (Hz)');
ylabel('Coh');
ylim([0 1.1]);

%% natural frequencies & damping ratioes
N = 2 ; % number of (dominant) modes to identify
% only use valide data (coherence above 0.9)
ind = Cxy>0.9;

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

T = 2x3 table
    Mode    Frequency     Damping 
    ____    _________    _________

     1       66.147      0.0081933
     2       307.36       0.031123

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

    

    

Answer 2

[oybozkurt]

Hello Mathieu,

I have been busy trying your programme for a long time, it is very good but I need your help with a few things. Firstly, in some cases it takes points very close to each other as 1st and 2nd mode. Secondly, is there a document you can recommend to understand the technical content of your programme?

[Mathieu NOE]

hello @oybozkurt and happy new year !

would you like to share your data and the code ? I can have a look and make suggestions. I have myself found that using modalfit options can be a bit tricky to get the correct results. Maybe I'll come with my own function some day

regarding publications, you can find a lot of publications on the internet regarding modal analysis, most of them will describe how time domain data are processed - here some suggestions (see also in attachment)

Modal Impact Postprocessing: Getting the Best FRF

https://www.google.fr/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwiPmZW84OCKAxVZKvsDHa6YLDoQFnoECBYQAQ&url=https%3A%2F%2Fpublications.jrc.ec.europa.eu%2Frepository%2Fbitstream%2FJRC66055%2Freqno_jrc66055_jrcseriesrap8.pdf%5B1%5D.pdf&usg=AOvVaw2i032J_6HZbLwbFukfeOdF&opi=89978449

[oybozkurt]

Thanks for your reply. I attach the 3 different data files and the code used (given by you). Some of the frequencies selected by the programme in these files are not compatible with the graph presented by the programme.

Also, as far as I understand from the programme, the damping coefficient is calculated with the half power method, is it possible for the programme to calculate it with log decrement at the same time? For the specimens used in the form of cantilever beam for the experiments, can the mode shapes be drawn if the impact hammer strike point and accelerometer position are entered?

[Mathieu NOE]

hello again @oybozkurt

attached is a slightly modified version of my code. I simply used the "peak piking" method for modalfit here as the other options seems not to give the correct results - as you experienced.

the code also shows now the N (here N=2) identified modes as blue diamonds overlaid on the magnitude plot

Maybe modalfit is not the best code for the job, I was about to try with this instead : Rational Fraction Polynomial Method - File Exchange - MATLAB Central

the half power method works in frequency domain data , for multiple modes

the log decrement method works directly on time data but works only for one single dominant mode - you cannot separate multiple modes frequencies and damping values with this method

this code could be extended to multiple sensors responses ,and from there you could draw the mode shapes , but I am not there yet. The initial purpose was to simply demonstrate how we can generate the FRF data from multiple hits data, not to re-create an entire post processing tool (like many already available , like ABRAVIBE)

all the best

[oybozkurt]

Hi Mathieu,

I just find a signalexpress software and I checked the graph of it with your code. There are some differences between them. Firstly, the magnitudes for vertical axis of frequency graph are different. Secondly, there is a difference on natural frequency and also it seems that differences are found between damping ratios. Could you comment on these? I want to be understand these.

Thank you

[Mathieu NOE]

hello again

your excel file does not contain any hammer test data - can you provide the right file please ?

[oybozkurt]

The sheet "Force" is from hammer.

[Mathieu NOE]

oh yes indeed - I didn't see the excel file was different with more sheets :)

so I changed the code the load the correct sheets and also to display the magnitude plot now in dB as in signalexpress , I also did some changes for figure 3 display

the commented lines are the previous code version (y log scale) - the new code is magnitude converted in dB ref to the declared units.

then I opted for new modalfit parameters :

N = 1 ; % number of (dominant) modes to identify

FR = [5 300]; %

I see the same values of the dominant mode as in signalexpress (fn = 30 Hz, T = -34 dB ref 1 g/N)

T = 1×4 table

Mode Frequency Damping Gain (dB)

____ _________ _______ _________

1 30.028 0.12593 -33.616

Code updated in attachment

Answer 3

[oybozkurt]

Hello Mathieu,

First of all, thank you for your efforts. How accurately does the modalfit command work? When the graphs for ‘FitMethod’,‘lsrf’ and ‘FitMethod’,‘pp’ are carefully examined, the natural frequencies do not exactly coincide with the peaks, and there are significant differences between damping evaluated using ‘lsrf’ and ‘pp’ for the data I attached. Again, is the damping calculated by the method in the picture I attached or by another method.

[Mathieu NOE]

hehe, I have made the same observations as you. And either I am not using modalfit the correct way but I found a bit annoying to have such difficulties to make that function deliver what she's supposed to do

this publication Modal Analysis Using the Signal Processing Toolbox of Matlab 2017 shows also that that function gives different damping values depending on the method you selected , more or less in accordance with a "professionnal" grade software , so each method comes with its pro and cons and you need to try / find out what works best for you.

although the demo code provided by TMW seem to show that modalfit is capable of delivering good results (at least on the demo data), we have here a case where it seems to give unreliable results.

I have not too much used this function up to now , as I mostly focus on finding frequencies rather than the exact damping level (and therefore I did not use modalfit in the past) , but I agree with you, I was not really impressed with the results of modalfit, unless (again) there is something I did no see how to use it in a better way.

[oybozkurt]

Hımm. What do you think about the frequency values (previously available in ms excel) presented by Signal Express software during the experiment? How can you calculate natural frequency and damping using these values?

[oybozkurt]

I meant without using modalfit.

[Mathieu NOE]

So I tried to get some info's about how modalfit works but the editable code does not really help me

fortunately other people have done more investigations than me and here comes back some tips :

https://www.bing.com/ck/a?!&&p=8d9f6d5451788f18df885bdd4f6dc839062aa2e08a9ec9f5fa6a1ca414f9416aJmltdHM9MTczNjcyNjQwMA&ptn=3&ver=2&hsh=4&fclid=11718eea-1a77-63c7-15dc-9a801bb76261&psq=modalfit+gives+different+damping+values&u=a1aHR0cHM6Ly9vcHVzNC5rb2J2LmRlL29wdXM0LWhhdy9maWxlcy8zNTk1L0kwMDEzMTU4ODZUaGVzaXMucGRm&ntb=1

PP (peak picking or peak amplitude as described in the picture) method is fine in our case as it's used for SISO systems - it's simple , fast and usually gives fairly good results. It's the same as the "-3 dB" method you showed in your picture above

other methods like LSCE and LSRF are intended to be used for SIMO or MIMO systems , so not really our job here , but I have to read again carefully those publications to see if there is anything that speaks again for a SISO case. I prefer not to say something wrong now.

this is also an interesting publication in attachment

maybe it's time I need to create my own modalfit function :)

[Mathieu NOE]

forgot to say that there are some info's on the help page of modal fit

modalfit

but here I see some differences between how mathworks implemented the PP method vs how it's described in the litterature - seems there is more than one way to define PP method

[Mathieu NOE]

hello again @oybozkurt

finally I tried to replace modalfit with my own version the peak picking / - 3 dB bandwith algorithm

and I performed some comparisons between the code with modalfit vs my own code

attached find all the code versions + one new function peakfinder (using it as alternative to findpeaks for faster execution - can be downloaded from Fex : peakfinder(x0, sel, thresh, extrema, includeEndpoints, interpolate) - File Exchange - MATLAB Central or simply copy the attached file in your path)

summary of codes with modalfit :

• Modal_Analysis_FRF_computationP0.m : to be used with data file Pr1.xlsx,Pr2.xlsx,etc...
• Modal_Analysis_FRF_computationT0.m : to be used with data file Test.xlsx

summary of codes with my own modal extraction :

• Modal_Analysis_FRF_computationP1.m : to be used with data file Pr1.xlsx,Pr2.xlsx,etc..
• Modal_Analysis_FRF_computationT1.m : to be used with data file Test.xlsx

And my results so far : (with my own code I got better frequency accuracy and the damping values are matching modalfit results when using the pp option) :

>> Modal_Analysis_FRF_computationP0 (with modal fit)

filename = "Pr1.xlsx";

Mode Frequency Damping

____ _________ ________

1 35.121 0.013558

2 286.1 0.011218

>> Modal_Analysis_FRF_computationP1 (my own peak picking / -3dB algo)

filename = "Pr1.xlsx";

Mode Frequency Damping Gain (dB)

____ _________ ________ _________

1 35.149 0.015418 -22.365

2 286.01 0.011156 -12.514

3 320.5 0.012275 -23.44

Modal_Analysis_FRF_computationT0 (with modal fit)

filename = "Test.xlsx";

Mode Frequency Damping

____ _________ _______

1 30.028 0.12593

2 200.61 0.19108

Modal_Analysis_FRF_computationT1 (my own peak picking / -3dB algo)

filename = "Test.xlsx";

Mode Frequency Damping Gain (dB)

____ _________ _______ _________

1 29.49 0.12159 -33.526

[oybozkurt]

Thank you. I prepared an MS Excel file to evaluate natural frequency and damping ratio using FRF data taken from Signal Express. I applied an equivalent command to peakfinder(max) to perform -3dB bandwidth algorithm. I will compare your codes results with my Excel file results.

[Mathieu NOE]

Just fyi

as we have now a code that works , it's time to make the script shorter and better looking. So the "heavy lifting" is done by functions and the main script is now much lighter and readable.

functions are :

dotheFFT.m : basically to go from the time records to the FRF

find_zc.m : zero crossing detection in signals (used indotheFFT.m)

peak_picking_results.m : extract modal parameters from FRFs (peak picking algorithm)

[Jo98]

Hello Mr Mathieu, I would like to know if you know about Stability Lobe Diagram? cause I have some questions about a code of SLD.

[Mathieu NOE]

hello @Jo98

let's try - what is the question ?

[Jo98]

Hello, it's only to know what do you think about this two approach for obtaining the SLD?

First method delete the negative values from blim variable (chip width), and the second method is applying the Fourier equations (Altintas method).

All of these codes are based on "Mechanical Vibrations" by Tony Schmitz, I'm gonna sending you the base code (book) and my own code.

The book codes are: p_4_8_1 and p_4_9_1.

The else is mine code, and I am using my own FRF instead of creating a signal from parameters (damping, stiffness and mass).

[Mathieu NOE]

hello again

I know SLD because I read quite a lot of publications about machining / chatter issues , but I don't realize SLD computations or predictions myself because that is more the problem of the consultant or tool suppliers or the machine designer.

I define my role as a supplier of damping systems that will improve the milling machine performance - in attachment you see the active damping system we have designed that allows usually +200% performance in rough milling.

I do not have the book you are mentionning , but maybe I can get it somehow.

Is your question related about if your code does exactly follow the book, I will be unable to tell right now as I don't have the book.

[Jo98]

Hello,

I understand it thanks for your help, I was trying to create a code on matlab for creating an SLD from the FRF signal.

As always thanks for your help.

Best Regards.

[Mathieu NOE]

ok - all the best for the future !