%  plot periodic square wave and its DFT

N = 20 ;    %  the period,  let's make it even
m = 7 ;     % 2m+1 1's in a period,  2m+1 <= N


np =2 ;	%  display np periods to left of origin, np to right, and one at origin 
nz = 3 ;    %  number of zeros to display to left and right of the np periods

t = [-nz - m - np * N : np * N + m + nz ] ;  % the times to display
x = [ zeros(1,nz) ones(1,2*m+1) ];
for i=1:2*np
	x = [x zeros(1,N-(2*m+1)) ones(1,2*m+1) ] ;
end
x = [x zeros(1,nz)] ;
subplot(2,1,1)
stem(t,x)
v = axis ;
axis([v(1) v(2) 0 1.5]) ;
ht = title('x[n]')
set(ht,'Position',[-np*N 1.54 17])
xlabel('n')

xp = [ones(1,m+1) zeros(1, N-(2*m+1)) ones(1,m) ] ;  % one period

XX = fft( xp ) / N ;
subplot(2,1,2)
stem([0:N-1], XX) ;
ht = title('X[k]')
vv = get(ht,'Position')
set(ht,'Position',[1 vv(2) vv(3)])
xlabel('k') ;
line([0 N-1], [0 0]) ;

ff = [0 : .01 : N-1] ;
envel = (1/N) * sin( (2*m+1) * ff *pi/N ) ./ sin ( ff*pi/N) ;
hold on
plot(ff,envel)
hold off
