issue during FFT of a discrete data

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.

Answer 1

[Mathieu NOE]

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

[hamid k]

Hello Mathieu,

Thank you for your response and consistent support. I am curious as to why the outcome of your code indicates a zero for the fundamental harmonic when we should expect to see either the 1st or 2nd order as the fundamental harmonic.

[Mathieu NOE]

hello Hamid

both methods do the same , the raw signal gets decomposed in Fourier coefficients , so first a dc term (mean value) + higher harmonics (1,2,3,...)

now we see that the first dominant peak correspond to the 4th order (or harmonic) , so I am not sure what you mean when you say "we should expect to see either the 1st or 2nd order as the fundamental harmonic".

for me "fundamental" is always harmonic / order 1 . but that does not imply it's amplitude is the largest - so normally people speak about the dominant frequency / harmonic / order which can be any integer between 0 (dc) and Inf , not only 1st or 2nd

but also, in many situations we don't really care about the fft/dft spectrum at f = 0 , so we could remove that point from the plot if you wish (and start x at 1 and not 0)

simply we shoudn't forget the dc term when you compare the raw and reconstructed signals , otherwise you can have an y offset problem if your dc term is not neglectable.

[Mathieu NOE]

yet another way to see that the dominant frequency is the 4th harmonic, it's simply because in one revolution (360°) your signal is the repetition (4 times) of (almost) the same profile that appears in each window of 90 °

[hamid k]

Hello Mathieu, thank you for your response. Why do you use virtual sample points to modify the original data in your method?

[Mathieu NOE]

hello Hamid, sorry I don't understand what you mean by "virtual sample points" ?

[Mathieu NOE]

hello Hamid

do you mind accepting may answer ? tx