application of fft analyzer

15 views (last 30 days)
Haithem
Haithem on 30 May 2023
Commented: Mathieu NOE on 30 May 2023
Hello
I'm currently studying a vibratory system, each time I pick up the acceleration as a function of time, to have easy-to-read curves, I've applied the fourier transform and then the magnitude and phase calculation to plot them, but I still had curves with a lot of noise and that didn't translate my values even though I applied a smoothing. Is there another solution to make a good study of vibratory systems, also you will find my matlab code here, thanks for your help.
% Lecture des données d'accélération à partir d'un fichier Excel
clc;clear;
a = xlsread('test(30,50,-100,0,0,0).xlsx', 'C:C');% colonne A contient les données d'accélération
Time = xlsread('test(30,50,-100,0,0,0).xlsx', 'A:A');
%Tracé de la courbe d'accélération en fonction du temps
%figure(1);
%plot(Time, a);
%xlabel('fréquence');
%ylabel('Accélération (m/s^2)');
%FFT
L = length(a);
Fs = find(Time == interp1(Time, Time , 1.0,'nearest'))-1; %compute sample freq, Fs
%Fs = 1 / (Time(2) - Time(1));
n=2^nextpow2(L);
%n = 2^24;
X=fft(a, n);
P2 = abs(X/L);
P1 = P2(1:n/2+1);
P1(:,2:end-1) = 2*P1(:,2:end-1);
P1 = P1/max(P1);
f = Fs*(0:(n/2))/n;
%figure(2);
%plot(f,P1, 'b')
%grid on
% Calcul de la magnitude et de la phase
magnitude = abs(X(1:L/2+1));
phase = angle(X(1:L/2+1));
% Tracé du diagramme de Bode
figure;
subplot(2, 1, 1);
semilogx(f(1:L/2+1), 20*log10(magnitude)); % Tracé de la magnitude en échelle logarithmique
xlabel('Fréquence (Hz)');
ylabel('Magnitude (dB)');
title('Diagramme de Bode - Magnitude');
subplot(2, 1, 2);
semilogx(f(1:L/2+1), rad2deg(phase)); % Tracé de la phase en échelle logarithmique
xlabel('Fréquence (Hz)');
ylabel('Phase (degrés)');
title('Diagramme de Bode - Phase');
%Lissage magnetude
%freq = f(1:L/2+1);
%mag = 20*log10(magnitude);
%windowSize = 100; % Taille de la fenêtre de moyenne mobile
%b = (1/windowSize)*ones(1,windowSize);
%a = 1;
%amplitude_lisse = filter(b, a, mag);
%figure(4);
%plot(freq,amplitude_lisse)
%grid on
%lissage de phase
%freq = f(1:L/2+1);
%PH = rad2deg(phase);
%windowSize = 100; % Taille de la fenêtre de moyenne mobile
%b = (1/windowSize)*ones(1,windowSize);
%a = 1;
%amplitude_lisse = filter(b, a, PH);
%figure(4);
%plot(freq,amplitude_lisse)
%grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lissage de l'amplitude par moyenne mobile
%windowSize = 100; % Taille de la fenêtre de moyenne mobile
%b = (1/windowSize)*ones(1,windowSize);
%a = 1;
%amplitude_lisse = filter(b, a, P1);
%figure(3);
%plot(f,amplitude_lisse)
%grid on

Accepted Answer

Mathieu NOE
Mathieu NOE on 30 May 2023
hello
there are multiple aspects in vibration / dynamic system analysis :
  • response signal analysis either in time or frequency domain (fft spectrum, spectrogram, etc...)
  • transfer functions (=> Bode diagrams) based on the analysis of input and output data
I don't see in your code how you can do a Bode plot of you system based only on one signal; you need both input and output to compute a transfer function
see example below (the data is attached)
clc
clearvars
data = load('beam_experiment.mat');
x = transpose(data.x); %input
y = transpose(data.y); %output
fs = data.fs; % sampling frequency
NFFT = 2048;
NOVERLAP = round(0.75*NFFT); % 75 percent overlap
%% solution 1 with tfestimate (requires Signal Processing Tbx)
% [Txy,F] = tfestimate(x,y,hanning(NFFT),NOVERLAP,NFFT,fs);
%% alternative with supplied sub function
[Txy,Cxy,F] = mytfe_and_coh(x,y,NFFT,fs,hanning(NFFT),NOVERLAP);
% Txy = transfer function (complex), Cxy = coherence, F = freq vector
% Bode plots
figure(1),
subplot(3,1,1),plot(F,20*log10(abs(Txy)));
ylabel('Mag (dB)');
subplot(3,1,2),plot(F,180/pi*(angle(Txy)));
ylabel('Phase (°)');
subplot(3,1,3),plot(F,Cxy);
xlabel('Frequency (Hz)');
ylabel('Coh');
%%% damping ratioes for modes
N = 2 ; % number of (dominant) modes to identify
[fn,dr] = modalfit(Txy,F,fs,N,'FitMethod','pp');
T = table((1:N)',fn,dr,'VariableNames',{'Mode','Frequency','Damping'})
%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = mytfe_and_coh(x,y,nfft,Fs,window,noverlap)
% Transfer Function and Coherence Estimate
% compute PSD and CSD
window = window(:);
n = length(x); % Number of data points
nwind = length(window); % length of window
if n < nwind % zero-pad x , y if length is less than the window length
x(nwind)=0;
y(nwind)=0;
n=nwind;
end
x = x(:); % Make sure x is a column vector
y = y(:); % Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;
Pxx = zeros(nfft,1);
Pyy = zeros(nfft,1);
Pxy = zeros(nfft,1);
for i=1:k
xw = window.*x(index);
yw = window.*y(index);
index = index + (nwind - noverlap);
Xx = fft(xw,nfft);
Yy = fft(yw,nfft);
Xx2 = abs(Xx).^2;
Yy2 = abs(Yy).^2;
Xy2 = Yy.*conj(Xx);
Pxx = Pxx + Xx2;
Pyy = Pyy + Yy2;
Pxy = Pxy + Xy2;
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);
Pyy = Pyy(select);
Pxy = Pxy(select);
else
select = 1:nfft;
end
Txy = Pxy ./ Pxx; % transfer function estimate
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;
end
  2 Comments
Haithem
Haithem on 30 May 2023
Thank you for your reply. I have used a swept frequency between 0 and 150 Hz to stimulate my system. Can I use this function directly [Txy, F] = tfestimate(x, y, hanning(NFFT), NOVERLAP, NFFT, fs); or must I declare that my frequency is variable during the simulation?
here's the code I've tried:
clc;
clearvars;
x = xlsread('testpos(16,-23,113,-26,0,0)sans_delais.xlsx', 'F:F'); % Input
y = xlsread('testpos(16,-23,113,-26,0,0)sans_delais.xlsx', 'B:B'); % output
fs = 300; % Sampling frequency
% Check the length of the signals
minLength = min(length(x), length(y));
NFFT = 500; % Reduce the segment length (choose a power of 2)
NOVERLAP = round(0.75 * NFFT); % 50 percent overlap (adjust if needed)
% Calculate transfer function using tfestimate
[Txy, F] = tfestimate(x, y, hanning(NFFT), NOVERLAP, NFFT, fs);
% Calculate coherence using mscohere
Cxy = mscohere(x, y, hanning(NFFT), NOVERLAP, NFFT, fs);
% Bode plots
figure(1);
subplot(3, 1, 1);
plot(F, 20 * log10(abs(Txy)));
ylabel('Mag (dB)');
subplot(3, 1, 2);
plot(F, 180/pi * angle(Txy));
ylabel('Phase (°)');
subplot(3, 1, 3);
plot(F, Cxy);
xlabel('Frequency (Hz)');
ylabel('Coh');
Mathieu NOE
Mathieu NOE on 30 May 2023
hello again
seems to me your code is correct but If you want me to test it, I need your data file as well
there is no need to declare what kind of signal or varaible frequency info to tfestimate
simply make sure that your excitation signal covers the right frequnecy range, has sufficient energy (amplitude) and for low damping systems use slow sweeps (or use random, random bursts signals)
tx

Sign in to comment.

More Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!