%% example_DFT_FFT_square.m
% EE-3221 Digital Signal Processing
% 5 May 2020
% Dr. Prust <prust@msoe.edu>

%% x[n] = 1 for n = 0,1,2,...,M-1
%% x[n] = u[n] - u[n-M]
%% X(e^{j Omega}) = sin(Omega M /2) / (sin(Omega/2)) * exp(-j Omega (M-1)/2)

close all
clear all

%% Plot DTFT using closed form solution
Omega = linspace(-2*pi, 2*pi, 1000);  % see periodicity
M = 8;
X_DTFT = sin(Omega*M/2)./ sin(Omega/2) .* exp(-j*Omega*(M-1)/2);
plot(Omega, abs(X_DTFT),'LineWidth',2) % use "plot" since we are approximating a continuous function
title('DTFT X(e^{j\Omega})')
xlabel('\Omega')
set(gca,'fontsize',14)

%% Calculate DFT (using DFT definition)
x = ones(1,M);
n = 0:(M-1);
N_0 = M;    % Will do a N_0 = M point DFT
Omega_0 = 2*pi/N_0;
k = 0:(N_0-1);  
X_DFT = zeros(1,N_0);
for ii=1:N_0
    X_DFT(ii) = sum(x.*exp(-j*k(ii)*(2*pi/N_0)*n));
end

%% Generate plot showing DTFT and DFT
figure;
plot(Omega, abs(X_DTFT), 'LineWidth',2)
hold on
stem(k*Omega_0, abs(X_DFT), 'LineWidth',2)
title('DTFT X(e^{j\Omega}) and DFT X[k]')
xlabel('\Omega')
set(gca,'fontsize',14)

%% Can get a finer (more dense) sampling of DTFT by "zero padding"
x = [ones(1,M) zeros(1,8)];  % Add zeros by "zero-padding"
N_0 = length(x);
n = 0:(N_0-1);
Omega_0 = 2*pi/N_0;
k = 0:(N_0-1);  
X_DFT = zeros(1,N_0);
for ii=1:N_0
    X_DFT(ii) = sum(x.*exp(-j*k(ii)*(2*pi/N_0)*n));
end

%% Generate plot showing DTFT and zero-padded DFT
figure;
plot(Omega, abs(X_DTFT), 'LineWidth',2)
hold on
stem(k*Omega_0, abs(X_DFT), 'LineWidth',2)
title('DTFT X(e^{j\Omega}) and zero-padded DFT X[k]')
xlabel('\Omega')
set(gca,'fontsize',14)

%% Compute 1024 point using FFT
N_0 = 1024;
Omega_0 = 2*pi/N_0;
k = 0:(N_0-1);
x = ones(1,M);
X = fft(x,N_0);
figure
plot(k*Omega_0, abs(X))
title('1024 point FFT')
xlabel('\Omega')
set(gca,'fontsize',14)

%% Same calculation, but plotted from -pi to pi
k = -N_0/2 : 1 : (N_0/2 - 1);
X = fftshift(fft(x,N_0));
figure
plot(k*Omega_0, abs(X))
title('1024 point FFT (with fftshift)')
xlabel('\Omega')
set(gca,'fontsize',14)
