Help with DFT Algorithm (No FFT)

9 views (last 30 days)
Devin Hunter
Devin Hunter on 3 Aug 2022
Commented: Devin Hunter on 3 Aug 2022
(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 Comments
Jan
Jan on 3 Aug 2022
Note: This is not twitter. No # before the tags. Thanks.

Sign in to comment.

Accepted Answer

Matt J
Matt J on 3 Aug 2022
Edited: Matt J on 3 Aug 2022
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)')
  1 Comment
Devin Hunter
Devin Hunter on 3 Aug 2022
That makes so much more sense. Thanks man!

Sign in to comment.

More Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!