

% for GNU Octave 

sinc_phs=0:1:1023;

osr=124.382;   % square wave frequency is 
               % sampling frequency divided by 
	       % this number 

square_phs=sinc_phs/osr;


sgn=((mod(square_phs,1)<0.5)*2)-1; % 
dte=mod(square_phs*2-0.0,1)-0.5;
sincterms=sinc(dte*osr*0.45); % cutoff frequency 
                    %  0.25-0.50 is a good range
		    %  0.5 is Nyquist frequency
		    %  but you may want something lower
blit=cumsum(sincterms.*sgn)-0.5;

plot( dte )
hold on
plot( sgn )
plot( sincterms )
plot( blit )

hold off

if 1    % dB spectrum plot BLIT
	wblit=blit.*hamming(length(blit))';
	fftblit=10*log10(abs(fft(wblit)));
	plot(fftblit(1:end/2));
end
if 1    % dB spectrum plot square wave
	wsgn=sgn.*hamming(length(sgn))';
	fftsgn=10*log10(abs(fft(wsgn)));
	hold on
	plot(fftsgn(1:end/2));
	hold off
	legend('BLIT','Naive Square')
end
