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


lenspec=len/2+1; # length of the spectrum

nblk=500
sumspec(1:lenspec,1:nblk)=0; # to keep each

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
	workt=fft(work); # take the transform
	ps=workt(1:lenspec,1:2) .* conj(workt(1:lenspec,1:2));
	t=ps(1:lenspec,1)+ps(1:lenspec,2); # add the two channels
	sumspec(1:lenspec,nblk)=t;
end

t=max(max(sumspec))
fflush(stdout);

sumspec=sumspec/t;
sumspec=max(sumspec,.0000000001);
sumspec=10*log10(sumspec);

surface(sumspec);