application of fft analyzer

조회 수: 15 (최근 30일)
Haithem
Haithem 2023년 5월 30일
댓글: Mathieu NOE 2023년 5월 30일
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

채택된 답변

Mathieu NOE
Mathieu NOE 2023년 5월 30일
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
Haithem
Haithem 2023년 5월 30일
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 2023년 5월 30일
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

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Spectral Estimation에 대해 자세히 알아보기

제품

Community Treasure Hunt

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

Start Hunting!

Translated by