Explore DFT Resolution, need help with answer interpretation

Hello there, I'm relatively new to Matlab and Digital Signal Processing. I'm trying to analyze the frequency resolution of the DFT as a function of the DFT size assuming a rectangular window. Using the a sum of sinusoids sampled at 5 kHz vary the frequency spacing

of the sinusoids to show the frequency resolution for at least 4 different DFT lengths, 256, 512, 1024, and 4096.

I got two versions of code from chatgpt. But I don't know how to interpret them.

1st version

% MATLAB code to analyze the frequency resolution of the DFT using dftmtx and a rectangular window

fs = 5000; % Sampling frequency (5 kHz)

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];

% Frequencies of the sum of sinusoids (vary frequency spacing)

f1 = 1000; % Frequency of the first sinusoid (1 kHz)

f_spacing = [1, 5, 10, 20]; % Frequency spacing between the two sinusoids

f_end = f1 + f_spacing; % Frequency of the second sinusoid

% Prepare figure

figure;

hold on;

for N = N_dft

figure; % Create a new figure for each DFT size

hold on;

for spacing = f_spacing

f2 = f1 + spacing; % Second sinusoid frequency

% Generate the sum of two sinusoids with frequencies f1 and f2

x = sin(2*pi*f1*t) + sin(2*pi*f2*t);

% Apply rectangular window (by taking the first N samples)

x_windowed = x(1:N); % Select the first N samples for DFT calculation

% Generate DFT matrix for size N using dftmtx

DFT_matrix = dftmtx(N);

% Manually compute the DFT using the DFT matrix

X = DFT_matrix * x_windowed(:); % Compute DFT of the windowed signal

% Compute the frequency axis for the current DFT

freq_axis = (0:N-1)*(fs/N);

% Plot the magnitude of the DFT

plot(freq_axis, abs(X), 'DisplayName', ['Spacing = ', num2str(spacing), ' Hz']);

end

% Add labels and legend

xlabel('Frequency (Hz)');

ylabel('Magnitude');

title(['DFT Magnitude Spectrum (N = ', num2str(N), ')']);

legend('show');

grid on;

hold off;

end

which gives me these resulted figures.

All the peaks are so similar that I can't distinguish them at all.

2nd version of code is here:

% Parameters

fs = 5000; % Sampling frequency (Hz)

T = 1; % Signal duration (seconds)

t = 0:1/fs:T-1/fs; % Time vector

% Sinusoids: Frequencies and amplitudes

f1 = 400; % Frequency of first sinusoid (Hz)

f2 = 450; % Frequency of second sinusoid (Hz)

f3 = 500; % Frequency of third sinusoid (Hz)

% Create the signal: sum of three sinusoids

x = sin(2*pi*f1*t) + sin(2*pi*f2*t) + sin(2*pi*f3*t);

% Define DFT lengths

DFT_lengths = [256, 512, 1024, 4096];

% Plot frequency resolution for each DFT length using a rectangular window

figure;

for i = 1:length(DFT_lengths)

N = DFT_lengths(i); % DFT size

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

% Compute the N-point DFT

X = (dftmtx(N)*x_windowed(:)).'; % Compute the N-point DFT

X_mag = abs(X); % Magnitude of DFT

% Frequency vector

f = (0:N-1) * (fs / N);

% Plot the DFT magnitude (first half)

subplot(length(DFT_lengths), 1, i);

plot(f(1:N/2), X_mag(1:N/2));

title(['DFT Magnitude Spectrum with Rectangular Window, N = ', num2str(N)]);

xlabel('Frequency (Hz)');

ylabel('|X(f)|');

grid on;

end

% Adjust plot layout

sgtitle('Frequency Resolution for Different DFT Sizes Using Rectangular Window');

%Part 5

% MATLAB code to analyze zero-padding in the DFT using a rectangular window

fs = 5000; % Sampling frequency (5 kHz)

t_duration = 1; % Signal duration in seconds

t = 0:1/fs:t_duration-1/fs; % Time vector

% Window sizes to analyze

window_sizes = [256, 512];

% Zero-padded DFT sizes to analyze

N_dft = [1024, 2048, 4096];

% Frequencies of the sum of sinusoids (vary frequency spacing)

f1 = 1000; % Frequency of the first sinusoid (1 kHz)

f_spacing = [5, 10]; % Frequency spacing between the two sinusoids

f_end = f1 + f_spacing; % Frequency of the second sinusoid

% Prepare figure

for window_size = window_sizes

figure; % Create a new figure for each window size

hold on;

for N = N_dft

for spacing = f_spacing

f2 = f1 + spacing; % Second sinusoid frequency

% Generate the sum of two sinusoids with frequencies f1 and f2

x = sin(2*pi*f1*t) + sin(2*pi*f2*t);

% Apply rectangular window (by taking the first window_size samples)

x_windowed = x(1:window_size); % Select the first window_size samples

% Zero-pad the signal if DFT size is larger than window size

x_padded = [x_windowed, zeros(1, N - window_size)];

% Generate DFT matrix for size N using dftmtx

DFT_matrix = dftmtx(N);

% Manually compute the DFT using the DFT matrix

X = DFT_matrix * x_padded(:); % Compute DFT of the windowed and zero-padded signal

% Compute the frequency axis for the current DFT

freq_axis = (0:N-1)*(fs/N);

% Plot the magnitude of the DFT

plot(freq_axis, abs(X), 'DisplayName', ['Spacing = ', num2str(spacing), ' Hz, N = ', num2str(N)]);

end

end

% Add labels and legend

xlabel('Frequency (Hz)');

ylabel('Magnitude');

title(['Zero-Padded DFT Magnitude Spectrum (Window Size = ', num2str(window_size), ')']);

legend('show');

grid on;

hold off;

end

with the second version of code I can see the differences among different sample sizes, as more sampling yields more sharp signals, but why is that and what differences it will make?

Answer 1

[Paul]

Hi Zhen,

The first code gives the expected result, you just need to zoom in to see it.

If we look at how freq_axis is defined, we see that the difference between successive frequencies is:

fs/N

fs = 5000; % Sampling frequency (5 kHz)
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];
% Frequencies of the sum of sinusoids (vary frequency spacing)
f1 = 1000; % Frequency of the first sinusoid (1 kHz)
f_spacing = [1, 5, 10, 20]; % Frequency spacing between the two sinusoids
f_end = f1 + f_spacing; % Frequency of the second sinusoid

The ratio fs/N is

fs./N_dft

ans = 1×4

19.5312 9.7656 4.8828 1.2207

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

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

which, I think, explains the plots below (along with the fact that as N increases we'll have points in freq_axis that get closer to f1 and f2).

% Prepare figure
for N = N_dft
    figure; % Create a new figure for each DFT size
    hold on;
    for spacing = f_spacing
        f2 = f1 + spacing; % Second sinusoid frequency
        
        % Generate the sum of two sinusoids with frequencies f1 and f2
        x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
        
        % Apply rectangular window (by taking the first N samples)
        x_windowed = x(1:N); % Select the first N samples for DFT calculation
        
        % Generate DFT matrix for size N using dftmtx
        DFT_matrix = dftmtx(N);
        
        % Manually compute the DFT using the DFT matrix
        X = DFT_matrix * x_windowed(:); % Compute DFT of the windowed signal
        
        % Compute the frequency axis for the current DFT
        freq_axis = (0:N-1)*(fs/N);
        
        % Plot the magnitude of the DFT
        plot(freq_axis, abs(X),'-o', 'DisplayName', ['Spacing = ', num2str(spacing), ' Hz']);
    end
    
    % Add labels and legend
    xlabel('Frequency (Hz)');
    ylabel('Magnitude');
    title(['DFT Magnitude Spectrum (N = ', num2str(N), ')']);
    legend('show');
    grid on;
    hold off;
    xlim([f1-10, f2+10])
end

[Zhen]

Hi Paul:
Thanks so much! You are a life saver! Endorsement!
I'll compare my codes with your codes and see if I can understand. I'll let you know if I still have questions.
Best!

[Zhen]

Hi Paul:

Thanks so much for your answer! Now I understand that when number of points in DFT increases, the frequency resolution is smaller, which means any frequency spacing under this resolution can be revealed. In this case I have two sin functions add up, which have two signals, and any frequency resolution fall below the spacing between two sin functions will clearly give us two peaks instead of one.

I do have one more questions about zero-padding:

I'm doing a project where I need to provide an analysis of zero padding impacts by using the sum of sinusoids sampled at 5 kHz vary the frequency spacing of the sinusoids and show DFT of lengths, 256, 512, 1024, and 4096, using the window sizes of 256, and 512, Assume a rectangular window.

Which means DFT size is larger than window size, and we are zero-padding the samples.

I got these two figures from my code, but I don't know how to interpret the impacts of zero-padding. It seems that at window size of 256, no matter how you increase DFT size, the results are always not distinguishable between two peaks of sinusoids. While my instructor said frequency accuracy depends on window size, and frequency resolution depends on DFT size. But here when window size is too small, we can't distinguish between peaks even the resolution is small.

Here is my code:

%Part 5

% MATLAB code to analyze zero-padding in the DFT using a rectangular window

fs = 5000; % Sampling frequency (5 kHz)

t_duration = 1; % Signal duration in seconds

t = 0:1/fs:t_duration-1/fs; % Time vector

% Window sizes to analyze

window_sizes = [256, 512];

% Zero-padded DFT sizes to analyze

N_dft = [1024, 2048, 4096];

% Frequencies of the sum of sinusoids (vary frequency spacing)

f1 = 1000; % Frequency of the first sinusoid (1 kHz)

f_spacing = [5, 10]; % Frequency spacing between the two sinusoids

f_end = f1 + f_spacing; % Frequency of the second sinusoid

% Prepare figure

for window_size = window_sizes

figure; % Create a new figure for each window size

hold on;

for N = N_dft

for spacing = f_spacing

f2 = f1 + spacing; % Second sinusoid frequency

% Generate the sum of two sinusoids with frequencies f1 and f2

x = sin(2*pi*f1*t) + sin(2*pi*f2*t);

% Apply rectangular window (by taking the first window_size samples)

x_windowed = x(1:window_size); % Select the first window_size samples

% Zero-pad the signal if DFT size is larger than window size

x_padded = [x_windowed, zeros(1, N - window_size)];

% Generate DFT matrix for size N using dftmtx

DFT_matrix = dftmtx(N);

% Manually compute the DFT using the DFT matrix

X = DFT_matrix * x_padded(:); % Compute DFT of the windowed and zero-padded signal

% Compute the frequency axis for the current DFT

freq_axis = (0:N-1)*(fs/N);

% Plot the magnitude of the DFT

plot(freq_axis, abs(X), 'DisplayName', ['Spacing = ', num2str(spacing), ' Hz, N = ', num2str(N)]);

end

end

% Add labels and legend

xlabel('Frequency (Hz)');

ylabel('Magnitude');

title(['Zero-Padded DFT Magnitude Spectrum (Window Size = ', num2str(window_size), ')']);

legend('show');

grid on;

hold off;

xlim([f1-10, f2+10])

end

I posed a new questions about it and you are very welcome to answer it there.

Again thanks so much for your answer, which is very clear and helps me understand frequency resolution very much!

Best

Zhen