- 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
application of fft analyzer
15 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Mathieu NOE
on 30 May 2023
hello
there are multiple aspects in vibration / dynamic system analysis :
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
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
More Answers (0)
See Also
Categories
Find more on Spectral Estimation in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!