% program to plot a pulse train and its spectrum using the FFT         
% (revision 2.0 ==> 4/99)
clear
bk=2*pi;
% typical values: 
Np=3;          % number of pulse widths in a pulse period
tau=1e-6;      % pulse width
fghz=0.01;     % carrier frequency in GHz
N=3;           % number of pulses in the train
Nt=256;        % number of time samples per pulse width 2^(integer)
ans=input('use default values (y/n)? ','s');
if ans=='n'
  fghz=input('carrier frequency in GHz: ');
  N=input('number of pulses in train: ');
  tau=input('input pulse width, tau: ');
  Np=input('pulse period in terms of pulse width: ');
  Nt=input('number of time steps: ');
end
T=Np*tau;       % period
fr=1/T;         % PRF
fc=fghz*1e9;
wc=bk*fc;       % carrier angular frequency
wo=bk*fr;       % PRF angular frequency
%-----------------------------------------------------
% plot a pulse train 
  dt=tau/Nt;      
  it=T/dt;
  k=0;
  for n=1:N
    for i=1:it
      k=k+1;
      to=(i-1)*dt;
      t(k)=to+(n-1)*T;
      S(k)=0;
       if to <= tau
        S(k)=sin(wc*to);
       end
    end
end
% there are now k=Nt*Np*N time samples; increase it by a factor of 
% NN before taking the fft (by zero packing)
NN=2;
S(k+1:NN*k)=zeros(1,(NN-1)*k);
t(1:NN*k)=dt*[0:NN*k-1];
A=fft(S,length(S)); Amax=max(abs(A));

figure(1)
subplot(211)
Smax=1;
disp(['number of time samples = ',num2str(k)])
plot(t/1e-6,S/Smax,'k'),xlabel('time (microseconds)'),ylabel('amplitude')
title(['pulse train packed with zeros (total of ',num2str(length(t)),' time samples)'])
grid,axis([0,max(t)/1e-6,-1,1])
subplot(212)
plot([1:length(A)],abs(A)/Amax,'k'),grid,xlabel('index'),ylabel('normalized magnitude')
title(['returned from FFT --  ',num2str(length(A)),'  frequency samples'])
F=fftshift(A);
delf=1/max(t); Fmax=max(abs(F));
f=[0:length(F)/2-1]*delf;
Fp=F(length(F)/2+1:length(F));

figure(2)
subplot(211)
plot([1:length(F)],abs(F)/Fmax,'k'),grid,xlabel('index'),ylabel('normalized magnitude')
title(['returned from FFTSHIFT --  ',num2str(length(A)),'  frequency samples'])
subplot(212)
plot(f/(1e9),abs(Fp)/Fmax,'k')
xlabel('frequency, GHz'),ylabel('normalized magnitude')
title('positive frequencies -- fft sample number converted to frequency')
grid 
 BW=10e6;  % bandwidth of frequencies to plot
 nn=0; Lim1=fc-BW/2; Lim2=fc+BW/2;
 for i=1:length(f)
	if f(i) >=Lim1 & f(i) <= Lim2
		nn=nn+1;
		fp(nn)=f(i); FFp(nn)=abs(Fp(i))/Fmax;
	end
 end
 subplot(212)
 plot(fp/(1e9),FFp,'k')
 xlabel('frequency, GHz'),ylabel('normalized magnitude')
 title('positive frequencies -- fft sample number converted to frequency')
 grid 
 axis([Lim1/1e9,Lim2/1e9,0,1])

figure(3)
subplot(211)
Smax=1;
disp(['number of time samples = ',num2str(k)])
plot(t/1e-6,S/Smax,'k'),xlabel('time (microseconds)'),ylabel('amplitude')
title(['pulse train packed with zeros (total of ',num2str(length(t)),' time samples)'])
grid,axis([0,max(t)/1e-6,-1,1])
 subplot(212)
 plot(fp/(1e9),FFp,'k')
 xlabel('frequency, GHz'),ylabel('normalized magnitude')
 title('positive frequencies -- fft sample number converted to frequency')
 grid 
 axis([Lim1/1e9,Lim2/1e9,0,1])