Explore DFT frequency accuracy, got stuck

I'm trying to analyze the frequency accuracy of the DFT as a function of DFT size assuming a rectangular

window. Using single sinusoidal waveform sampled at 5 kHz vary its frequency in small

increments to show the how the frequency accuracy changes for at least 4 different DFT length,

256, 512, 1024, and 4096.

Here is my code:

%Part 1

% MATLAB code to analyze the frequency accuracy of the DFT manually

fs = 5000; % Sampling frequency 5 kHz

f_start = 1000; % Start frequency of sinusoid (1 kHz)

f_end = 1050; % 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

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;

However, the result I got is really weird.

At some frequency the estimate error between true frequency can be 3000, while most of the time it's under 1.

However, it I substitue this line: X = x_windowed*dftmtx(N);

with this line: X = fft(x_windowed, N);

Then I got this image

It seems that this is the correct figure.

I read from this website that fft and dftmtx function can be interchangeble

https://www.mathworks.com/help/signal/ref/dftmtx.html

But my results doesn't approve it. I want to know why.

Besides, am I applying windowing function in a correct way?

The code is as following:

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;

Answer 1

[Paul]

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)

ans = 1×2

1.092540172992740e+02 1.092540172992740e+02

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

temp(1) - temp(2)

ans = 
     1.421085471520200e-14

freq_axis(idx(1:2)) 

ans = 1×2

1.0e+00 * 4.003906250000000e+03 9.960937500000000e+02

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Maybe the code should be changed to only search for the max peak at frequencies below the Nyquist frequency.

[Zhen]

Hi Paul:

Thanks so much for your quick response! You saved my life! Actually I saw your deleted answer and modified my code so that now ths FFT function and dftmtx yield the same plots.

My question is: how is your answer now related to your previous deleted answer?

Is this slightly difference between first and second peak of sinusoid resulting in the frequency error figure?

[Paul]

I'm glad I could save your life.

My deleted answer was correct in that the DFT can be computed as

X1 = dftmtx(N)*x(:)

but it was incorrect in that the DFT can also be computed as (assuming x is a row vector)

X2 = x*dftmtx(N);

In theory, the elements of X1 and X2 will be equal to each other due to the symmetry of the DFT matrix. In practice, they might not be depending on how mtimes, * works for those two cases.

I don't know why the first approach (X1) always yields the same results as using fft, nor do I know why fft always returns (in your examples) the first spectral peak as being just slightly higher than the second.

The best approach for this particular problem is to limit the search for the max peak to frequencies at or less than the Nyquist frequency so you don't have to worry about issues related to numerical accuracy.