Explore DFT Resolution, need help with answer interpretation

6 次查看(过去 30 天)
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?

采纳的回答

Paul
Paul 2024-10-7
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
  2 个评论
Zhen
Zhen 2024-10-7

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
Zhen 2024-10-8
编辑:Zhen 2024-10-8
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

请先登录,再进行评论。

更多回答(0 个)

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by