clear all
close all
fclose all

len=2^12;
fs=44100;

wind=hann(len); # make the hann window
slen=len/2; # the shift length, i.e. 1/2 overlap.

x=wavread('01jt.wav');  # or you can specify any other wave file, of course.
whos   # shows the form of the file, n samples by 2 channels

flen=length(x); # gets the length of the file in samples

nblk=floor(flen*2/len)-1 # might clip the end of the file
nonzeroblocks=0;

lenspec=len/2+1; # length of the spectrum
sumspec(1:lenspec,1:2)=0; # to keep overall power spectrum

for ii = 1:1:nblk
	work=x( ((ii-1)*slen +1):((ii+1)*slen) , 1:2); # pick one block of data
	work(1:len,1)=work(1:len,1) .* wind;  # window left channel
	work(1:len,2)=work(1:len,2) .* wind;  # ditto right channel
	ss=sum(sum(abs(work)))   # is there any energy here?
	if ss > 1/32768  # if it's smaller than that there is no signal.
		nonzeroblocks=nonzeroblocks+1;
		workt=fft(work); # take the transform
		ps=workt(1:lenspec,1:2) .* conj(workt(1:lenspec,1:2));
		sumspec=sumspec + ps; # keep overall power spec sum.
		psmax=max(max(ps)); # find overall maximum
		ps=ps/psmax; #normalize to peak for presentation
		ps=max(ps,.00000000001); # put in minimum energy to avoid stuff
		ps=log10(ps)*10; # dB spectrum.
		subplot(2,1,1);
		plot(work);
		subplot(2,1,2)
		plot(ps);

	end

	fflush(stdout);  # this slows things down but makes sure
			 # we see the output from the program
	junk=input('hit cr');

end