Hi Zhen,
Consider the following, which is one case from the overall analysis.
%Part 1
% MATLAB code to analyze the frequency accuracy of the DFT manually
fs = 5000; % Sampling frequency 5 kHz
f_start = 1002; % Start frequency of sinusoid (1 kHz)
f_end = 1002; % End frequency of sinusoid (1.05 kHz)
f_step = 1; % Frequency step size (small increments)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% DFT sizes to analyze
N_dft = [256, 512, 1024, 4096];
% Prepare figure
%figure;
%hold on;
for N = N_dft(1)
freq_accuracy = []; % Store accuracy results for each DFT size
for f = f_start:f_step:f_end
% Generate sinusoidal waveform
x = sin(2*pi*f*t);
rectangular_window = rectwin(N)';
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;
% Manually compute the DFT
X = x_windowed*dftmtx(N);
% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);
% Find the index of the peak in the magnitude of the DFT
[~, peak_idx] = max(abs(X));
% Estimated frequency from the DFT
f_estimated = freq_axis(peak_idx);
% Store the frequency error (absolute difference between true and estimated frequency)
freq_error = abs(f_estimated - f);
freq_accuracy = [freq_accuracy, freq_error];
end
% Plot the frequency accuracy for the current DFT size
%plot(f_start:f_step:f_end, freq_accuracy, 'DisplayName', ['N = ', num2str(N)]);
end
% Add labels and legend
%xlabel('True Frequency (Hz)');
%ylabel('Frequency Error (Hz)');
%title('Frequency Accuracy of the DFT for Different DFT Sizes (Manual DFT)');
%legend('show');
%grid on;
%hold off;
Plotting the spectrum shows there are two peaks, as expected for a sinusoid.
It turns out that sometimes the second peak is very, very slightly higher than the first, presumbably just due to numerical effects.
figure
plot(freq_axis,abs(X));
grid
xline(f_estimated)
We can see that the first peak is just slightly smaller than the second.
[temp,idx] = sort(abs(X),'descend');
format long e
temp(1:2)
temp(1) - temp(2)
freq_axis(idx(1:2))
Maybe the code should be changed to only search for the max peak at frequencies below the Nyquist frequency.
