% IMPLEMENTATION OF DFT USING TRIGONOMETRIC FORM
% N = Number of samples
% T = Record length in time
% y = Sampled signal
% 
clear
N=128;
T=5;
dt=T/N;
t=(1:N)*dt;
f0=30;
M=N/2;
y(1,1:M)=ones(1,M); y(1,M+1:N)=zeros(1,N-M);
y=sin(2*pi*f0*t).*y;


% Compute arguments of trigonometric functions
for n=1:N
    for p=0:N/2
        Phi(p+1,n)=2*pi*p*n/N;
    end
end  

% Compute Fourier Coefficients
for p=1:N/2+1
    A(p)=2/N*y*cos(Phi(p,:))';
    B(p)=2/N*y*sin(Phi(p,:))';
end
A(1)=2/N*sum(y);
A(N/2+1)=A(N/2+1)/2;
B(N/2+1)=0;

% Reconstruct Signal - pmax is number of frequencies used in increasing order
pmax=N/2;
for n=1:N
    ynew(n)=A(1)/2;
    for p=2:pmax
        ynew(n)=ynew(n)+A(p)*cos(Phi(p,n))+B(p)*sin(Phi(p,n));
    end
end    

% Plot Data
figure(1)
plot(t,y,'o')

% Plot reconstruction over data
hold on
plot(t,ynew,'r')
xlabel('Time')
ylabel('Signal Height')
title('Reconstructed Signal')
hold off

% Compute Frequencies
n2=N/2;
f=(0:n2)/(n2*2*dt);

% Plot Fourier Coefficients
figure(2)
subplot(2,1,1)
stem(f,A)
xlabel('Frequency')
ylabel('Amplitude')
title('A')
subplot(2,1,2)
stem(f,B)
xlabel('Frequency')
ylabel('Amplitude')
title('B')

% Plot Fourier Spectrum
figure(3)
Power=sqrt(A.^2+B.^2);
stem(f,Power(1:n2+1))
xlabel('Frequency')
ylabel('Power')
title('Periodogram')

% Show Figure 1
figure(1)