I'm trying to compare the estimation of the power spectral density (in this case of a fake cosine signal) by two different methods: 1st, MATLAB's fft function, and 2nd, a direct calculation from the definition, $\int f(t)\cdot e^{-i\omega t}\mathrm{d}t$. The problem I find is that, for whatever reason, the FFT result is slightly shifted towards the right in frequency, and I'm not able to see why.

I'm adding the full code so that if anyone wants to try it, it should directly work and plot the results in MATLAB, something like

>> compare_FT(10);

As is clear from the code, the peak should be exactly at $\omega = 6$: that is indeed the case for the calculation from definition, but not exactly so for the FFT.

enter image description here

function[] = compare_FT(fs)

%% signal parameters (now fake, then from lab)
f = 6/(2*pi);
w_d = 2*pi*f;
T = 100;
num_samples = T*fs;

%% vector definitions
t = linspace(0, T, num_samples);
xx = cos(w_d*t); % the "fake" signal
dt = t(2) - t(1); % or, equivalently, dt = 1/fs
% omega vector for DTFT in a specific range:
ww_size = num_samples; % the bigger this size, the more points in the estimated FT
w_start = 5; % start calculating the FT at this omega
w_end = 7;  % stop at this omega
ww = linspace(w_start, w_end, ww_size);
Sx = zeros(1, ww_size); % prealocate space for integration

% DTFT method
for i = 1:ww_size
    int_value = trapz(xx.*exp(-1i*ww(i)*t))*dt;
    Sx(i) = int_value.*conj(int_value)/T; % estimate PSD

%% FFT method
n_fft = 2^14; % number of points in FFT (more points -> better FT discretization)
signal_fft = fft(xx, n_fft); % the fft itself
X2_fft = abs(signal_fft).^2; %signal_fft.*conj(signal_fft);
L = length(X2_fft);
Sx_fft = X2_fft(1:L/2+1)/(fs*num_samples); % proper scaling FFT to estimate PSD
w_fft = 2*pi*fs*(0:L/2)/L; 

%% Plot results to compare both methods
hold on;
% Plot theoretical value at top
plot(ww, ones(1,length(ww))*T/4, '--', 'Color', 'g', 'LineWidth', 3.5);
% Plot "direct" calculation
plot(ww, Sx, '.', 'Color', 'b', 'LineWidth', 3.5);
% Plot "FFT" calculation
plot(w_fft, Sx_fft, '-', 'Color', 'r', 'LineWidth', 3.5);
set(gca,'FontSize',18,'FontName', 'CMU Sans Serif');
xlabel('\omega', 'fontsize', 24, 'FontName', 'CMU Sans Serif');
ylabel('$S_x$', 'fontsize', 24, 'Interpreter', 'latex', 'FontName', 'CMU Sans Serif');
xlim([w_start, w_end]);
grid on;

Your vector

t = linspace(0, T, num_samples);

is the problem. It should be

t = (0:num_samples-1)/fs;

Note that, even though you claim it in a comment in the code, in the current version, dt = t(2) - t(1); is not equal to 1/f_s.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.