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
