%Pulse Spectrum display
% An m file for the display of a spectrum for EE3700 tutorial problem1
% A pulse of normalise width PulseW (=0.5 for squarewave) and amplitude 1
% us used. The m file caalculates the corresponding spectrum.
% Author C. J. kikkert March 2005
clear;
OS=4; %Oversampling ratio
NB=1024;   %Basic length of period
NB2=NB/2;
N= OS*NB; % FFT length can change tis for better accuracy if needed.
PulsW=0.334; %Normalised Pulse Width 
PW=2*round((PulsW*NB-1)/2);    %odd pulse width centred around time =0
HPW=PW/2;
PW=PW+1;
PWP=PW/NB
Ft=zeros(1,N);  %fill time array with zeros
for II=1:OS    
    for I=1:HPW
            Ft((II-1)*NB+1)=1;   %Centre of pulses
            Ft((II-1)*NB+I+1)=1;
            Ft((II)*NB-I+1)=1;
            Ft((II-1)*NB+NB2+1)=-1;   %Centre of pulses
            Ft((II-1)*NB+NB2+I+1)=-1;
            Ft((II)*NB-NB2-I+1)=-1;
    end;
end;
Ft(1)=1;    %DC component Matlab arrays from 1 to N
figure(1);
str=['Time Waveform'];
plot(real(Ft(1:N)));
axis([0 N -1.1 1.1])   
xlabel('Time (Samples)');ylabel('Amplitude');
grid on
title(str);
Ff=fft(Ft)/N;
Ff=fftshift(Ff);
figure(2);
str=['Spectrum Waveform'];
Span=round(4*N/PW);
if (Span>(N/2-1)) Span=N/2-1; end;
Fmax=N/2+Span;
Fmin=N/2-Span;
for I=1:2*Span+1
    Freq(I)=(I-Span-2)/OS;
end;
plot(Freq,real(Ff(Fmin:Fmax)));
xlabel('Frequency');ylabel('Amplitude');
grid on
title(str);