Help with DFT Algorithm (No FFT)
显示 更早的评论
(Please refer to my first comment for further details regarding the issue) I am currently programming a script that plots the magnitude spectrum of a sinusoid, x(t) = 16cos(2*pi*100t), at Fs = 2kHz using discrete fourier transform (DFT) without using the FFT function as a way to fully understand the concept. The program is mostly complete, but I suspect there is something wrong with the way that I am collecting the complex numbers in the sum as the plot for magnitude spectrum does not represent what I believe it should represent once additional zeros and cycles are incorporated (see my first comment). Any advice to put me on the right track would be appreciated. The code that executes the DFT computation will be displayed below:
%% Input Vector Initialization Process
Fs = 2000; % Hz
n = 0:20-1;
x = 16 * cos(2 * pi * (1/20) * n); % Input Vector
%% Adding Additional Cycles (if needed)
cycles = 20;
xx = x; % x (expanded) = x
if cycles > 1
for j=2:cycles
xx = [xx x];
end
end
% Number of Zeroes Padded for x(n)
Z = 200;
%% Zero-Padding Process
xx = [xx zeros(1,Z)];
%% DFT Algorithm
X = []; %X(omega)
N = length(xx);
sum = 0;
df = Fs / N;
fr = (0:N-1)*df;
for k=0:N-1
for n=0:N-1
sum = sum + xx(n+1) * exp(-1i*2*pi*k*n / N);
end
X = [X sum];
sum = 0;
end
%% Plots
% Plotting Within Nyquist Range
nyq = ceil(length(X) / 2);
X = X(1:nyq);
fr = fr(1:nyq);
stem(fr,abs(X))
xlabel('Frequency (Hz)')
ylabel('|X(\omega)|')
title('Amplitude Spectrum of X(\omega)')
3 个评论
Matt J
2022-8-3
In what way does it seem "off"?
Devin Hunter
2022-8-3
I apologize for my previous vagueness. When I referenced the plot of the amplitude spectrum as off. I meant that its peak value is different than what I expected when zeros and additional cycles are incorporated. For example, when the number of padded zeros, 
, and the number of cycles within the function x is 
. This is what my spectrum looks like:
, and the number of cycles within the function x is . This is what my spectrum looks like:
The number of samples, N, of function x that are used in the calculation of the DFT is 
where 
# of points in single cycle (20) × # of cycles (20) + # of padded zeros (200). To my knowledge of DFT, the peak value of the amplitude spectrum is equal to 
where 
amplitude of the sampled sinusoid. When I carry out this calculation I get the following: 
. I notice that this is considerably far off from what the peak value in the figure appears to be. I have an idea that there may be something wrong with way that I am calculating the DFT of the signal, but I am not sure how. Any advice on how I should go about this would be appreciated.
where # of points in single cycle (20) × # of cycles (20) + # of padded zeros (200). To my knowledge of DFT, the peak value of the amplitude spectrum is equal to where amplitude of the sampled sinusoid. When I carry out this calculation I get the following: . I notice that this is considerably far off from what the peak value in the figure appears to be. I have an idea that there may be something wrong with way that I am calculating the DFT of the signal, but I am not sure how. Any advice on how I should go about this would be appreciated.
Jan
2022-8-3
Note: This is not twitter. No # before the tags. Thanks.
采纳的回答
When I carry out this calculation I get the following: 
.
.No, because 1/3 of your xx are padded zeros. This will reduce your amplitude by a factor of 2/3, giving 3200.
Had you generated 30 cycles - giving you a duration of 600 samples with no zero padding - you would have seen an amplitude of 4800 as you predicted. This is confirmed below.
%% Input Vector Initialization Process
dT = 1; %<---Edited
n = 0:dT:20-1;%<---Edited
x = 16 * cos(2 * pi * (1/20) * n); % Input Vector
%% Adding Additional Cycles (if needed)
cycles = 30;
xx = x; % x (expanded) = x
if cycles > 1
for j=2:cycles
xx = [xx x];
end
end
% Number of Zeroes Padded for x(n)
Z = 0;
%% Zero-Padding Process
xx = [xx zeros(1,Z)];
N=numel(xx);
df=1/N/dT;%<---Edited
%% DFT Algorithm
F=exp(-1i*2*pi*(0:N-1)'.*(0:N-1) / N) ; %DFT matrix
X = F * xx(:);
%% Plots
% Plotting Within Nyquist Range
nyq = ceil(length(X) / 2);
X = X(1:nyq);
fr = (0:N-1)*df;
fr = fr(1:nyq);
stem(fr,abs(X))
xlabel('Frequency (Hz)')
ylabel('|X(\omega)|')
title('Amplitude Spectrum of X(\omega)')
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Spectral Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
