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)')

[Matt J]

In what way does it seem "off"?

[Devin Hunter]

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:

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.

[Jan]

Note: This is not twitter. No # before the tags. Thanks.

Answer 1

[Matt J]

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)')

[Devin Hunter]

That makes so much more sense. Thanks man!