clc;
clear all;
close all;

fname='idiot.wav'


len=1024;

wind=hann(len);

fs=44100;

fs2=fs/2;

step=2^.2;
lpow=1/3.5;


low(1:100)=0;
high(1:100)=0;
nbins=0;

fl=20;
fhigh=70;
nbins=nbins+1;
low(nbins)=fl;
high(nbins)=fhigh;

f=fhigh;
while (f<fs2)
	f=fhigh+50;
	if ( fhigh*step > f)
		f=fhigh*step;
	end
	nbins=nbins+1;
	low(nbins)=fhigh;
	high(nbins)=f;
	fhigh=f;
end

nbins=nbins-1; % we can ignore the last bin.
low(nbins);
high(nbins);

for ii=1:nbins
	low(ii)= round(low(ii)/fs*len)+1;
	high(ii)=round(high(ii)/fs*len)+1;
end
high(nbins)=len/2+1;   % lump everything into that last ERB, thank you
pl(1:nbins)=0;

sig=wavread(fname);
slen=length(sig);

shift=len/2;

nblks=round(slen/shift-1.999)


hloud(1:nblks+2)=0;

powblocks=0;
lmin=3;

for iblk=1:nblks
	l=sig( ((iblk-1)*shift+1):(iblk+1)*shift,1);
	r=sig( ((iblk-1)*shift+1):(iblk+1)*shift,2);
	l=l .* wind;
	r=r .* wind;
	lt=fft(l);
	rt=fft(r);
	pl(1:nbins)=0;
	for ii=1:nbins
		pl(ii)=0;
		pl(ii)=pl(ii)+sum(lt(low(ii):high(ii)).*conj(lt(low(ii):high(ii))));
		pl(ii)=pl(ii)+sum(rt(low(ii):high(ii)).*conj(rt(low(ii):high(ii))));
		pl(ii)=(pl(ii)+.0001).^lpow;
	end
	temp=sum(pl);
	if (temp > lmin)
		powblocks=powblocks+1;
		hloud(powblocks)=temp;
		[iblk,temp];
		fflush(1);
	end
end

hloud(powblocks+1)=0;
hloud(powblocks+2)=400;
 

loudmean=sum(hloud(1:powblocks))/powblocks

dl(1:powblocks-1)=0;

for ii=1:powblocks-1;
	dl(ii)= 2*(hloud(ii+1)-hloud(ii))/(hloud(ii+1)+hloud(ii));
end


%dl= ( hloud(1:powblocks-1) - hloud(2:powblocks))/loudmean;

dh=sum(abs(dl))/(powblocks-1)
dlmax=max(dl);
dlmin=min(dl);

[b1,b2]=hist(hloud(1:powblocks+2),401,1);


tt=0;
ii=1;
while tt<.05
	tt=tt+b1(ii);
	ii=ii+1;
end
tt

ht=0;
ij=401;
while (ht<.05)
	ht=ht+b1(ij);
	ij=ij-1;
end
ht

subplot(4,1,1)
bar(b2,b1);
hold
ttm(1:401)=0;
t=max(b1);
ttm(ii)=-t/4;
ttm(ij)=-t/4;
ttm(round(loudmean))=-t/4;
bar(b2,ttm);
hold

tmp=sum(sum(sig .* sig));
lll=length(sig);
rms=sqrt(tmp/2/lll)
pk=max(max(abs(sig)))
pkr=pk/rms


subplot(4,1,2)

plot(1:powblocks,hloud(1:powblocks));
axis([ 1 powblocks 0 400 ]);

subplot(4,1,3)
plot(1:powblocks-1,dl);
axis([ 1 powblocks -1 1]);

subplot(4,1,4)
sig=round(sig*32768);

his=histc(sig(:,1),-32768:32767); % channel 1
his=his+histc(sig(:,2),-32768:32767); %channel 2

tot=sum(his);
his=his/tot;
his=max(his, .0000000001);
xax=-32768:32767;

semilogy(xax,his);
axis([-40000 40000 1e-8 1e-1]);

infostr=sprintf('filel %s mean loudness = %6.2f mean block to block change = %6.2f rms = %8.5f peak = %8.5f peak to rms = %6.2f number of blocks = %d of length %d samples',fname, loudmean, dh, rms, pk, pkr, powblocks, len);
xlabel(infostr);
