issue during FFT of a discrete data

조회 수: 7 (최근 30일)
hamid k
hamid k 2023년 11월 29일
댓글: Mathieu NOE 2023년 12월 11일
Hello, everyone.
I’m using Ansys Maxwell, which is finite element software. I plotted Force (N) versus Angular position (Deg), and the software calculated FFT, but it seems the result is not correct. Therefore, I have exported the Force (N) versus Angular position (Deg) data of the wave to the Matlab workspace. I have attempted FFT of discrete data but have been unsuccessful. How can I create an FFT with a wide harmonic order and magnitude?
You can find the data in the attached Excel file.

답변 (1개)

Mathieu NOE
Mathieu NOE 2023년 11월 29일
hello Hamid and welcome back again
you can see below that fft and my home made DFT computations give the same spectrum. Also notice here that I limited the computation to orders from 0 (dc term) to 128; you can of course change that upper limit, but if you reduce it you will notice that the reconstructed signal (in red) will be more approximate (as you would get by applying a low pass filter)
once we are in the frequency domain, the x vector is the inverse of your angular (revolution) - we can either speak in "order" (orders extraction, synchronous fft it's all the same stuff) or in 1/rev frequencies which is identical.
raw and reconstructed signals :
fft vs home made dft spectrum computation :
notice they are identical !
code :
%Load Excel file
data = readmatrix('Force_deg.xlsx'); % Pos(deg) Force(N)
%Extract theta and y columns
theta = data(:,1); % theta (0 - 360 deg)
y = data(:,2); % y data
n = numel(theta);
%%%%%%%%%%%%% main code %%%%%%%%%%%%%%%%%
% orders extraction (DFT)
orders = 128;
% model fit : X = A*cos(order*theta) + B*sin(order*theta) + C
C = mean(y);
yfitsum = C; % init yfitsum with the dc term
yc = y-C;
for k = 1:orders
A(k) = 2/n*trapz(yc.*cosd(k*theta));
B(k) = 2/n*trapz(yc.*sind(k*theta));
tmp = A(k)*cosd(k*theta) + B(k)*sind(k*theta);
yfitsum = yfitsum + tmp;
end
figure(1),
plot(theta, y, 'b',theta, yfitsum, 'r')
legend('data','model fit - all orders');
% FFT plot
figure(2)
% dft (home made)
AB_mag = [C sqrt(A.^2+B.^2)]; % so from dc to "orders" max frequency
subplot(2,1,1),plot(0:orders, AB_mag, 'bo-')
title('Signal Spectrum - my DFT')
xlabel('orders');
ylabel('magnitude');
% with fft
[f1,fft_spectrum1] = do_fft(theta./max(theta),y); % nb time vector is here replaced by 1 revolution angular vector
fft_spectrum1(1) = fft_spectrum1(1)/2; % dc term amplitude must be divided by 2
ind = (f1<=orders);
subplot(2,1,2),plot(f1(ind),fft_spectrum1(ind),'-*')
title('Signal Spectrum - matlab fft')
ylabel('|X(f)|')
xlabel('Frequency[1/revolution]')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [freq_vector,fft_spectrum] = do_fft(time,data)
time = time(:);
data = data(:);
dt = mean(diff(time));
Fs = 1/dt;
nfft = length(data); % maximise freq resolution => nfft equals signal length
%% use windowing or not at your conveniance
% no window
fft_spectrum = abs(fft(data))*2/nfft;
% % hanning window
% window = hanning(nfft);
% window = window(:);
% fft_spectrum = abs(fft(data.*window))*4/nfft;
% 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
  댓글 수: 6
Mathieu NOE
Mathieu NOE 2023년 12월 4일
hello Hamid, sorry I don't understand what you mean by "virtual sample points" ?
Mathieu NOE
Mathieu NOE 2023년 12월 11일
hello Hamid
do you mind accepting may answer ? tx

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

카테고리

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

Community Treasure Hunt

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

Start Hunting!

Translated by