%This file calculates the Fourier series coefficients numerically
Omegao=1.0e3;
No=20; %20 samples per cycle
To=2*pi*1.0e-3; %periodic time
T=To/No;  %this is the sampling interval
Data=[ones(1,5)  0.5  zeros(1,9)  0.5  ones(1,5)]'; %This is the shape of the signal
OmegaD=Omegao*T;  %this is the discrete time frequency (rad/sample)
D=zeros(101,1); %matrix for storing the values of the Fourier Series coefficients
%start of loop for caluculating the coefficients
for n=-50:1:50
    Sum=0; %initialize summation
    for k=1:20  %repeat for all data points
        Sum=Sum+Data(k)*exp(-i*n*OmegaD*(k-1));
    end
    D(n+51,1)=Sum/No; %store the coefficients
end
plot(real(D));