function [A,grk,gtk] = statcom2(A,bsize);

%-------------------------------------------------------------
% This function performs a statistical quantization of
% the normalized transformed matrix. The function
% takes the normalized transformed matrix and a vector
% which defines the number of bits to be considered
% at each level of the quantization
%-------------------------------------------------------------
u=sym('u');
countr=1;
countt=1;
clear grk
clear gtk
tnum=length(bsize)-1; % number of transformations
n=length(A); % dimension of input matrix

sizesq(1)=n/(2^tnum); % sizesq defines the sizes of the
sizesq(2)=sizesq(1); % various submatrices formed during
% the transformation
for i=3:tnum+1
sizesq(i)=2*sizesq(i-1);
end

dimhold=0; % dimhold will be used to locate the
% submatrix inside of the matrix
for i=1:tnum+1
dim=sizesq(i); % dimension of submatrix
bnum=bsize(i); % number of bits used to represent submatrix

if i==1 % at the highest level there is only one
k=1; % submatrix (upper lefthand corner)
else
k=3; % at all other levels there are 3 submatrices
end

for j=1:k
clear square

if j==1 % defines which submatrix will be quantized
square(1:dim,1:dim)=A(dimhold+1:dimhold+dim,dimhold+1:dimhold+dim);
elseif j==2
square(1:dim,1:dim)=A(dimhold+1:dimhold+dim,1:dim);
else
square(1:dim,1:dim)=A(1:dim,dimhold+1:dimhold+dim);
end

if bnum==0 % handles zero bit representation
square=zeros(dim,dim);
grk(countr)=0;
gtk(countt)=0;
countr=countr+1;
countt=countt+1;

else
square
sqmean=sum(sum(square))/(dim^2); % finds mean of submatrix
sqvar=std(square(:))^2; % finds variance of submatrix
tk=zeros(2^bnum+1,1);
rk=zeros(2^bnum,1);
tk(1)=min(min(square));
range=max(max(square))-min(min(square));
gtk(countt)=tk(1);
countt=countt+1;
for p=1:2^bnum
tk(p+1)=tk(1)+range*int(exp(((-1/3)*-(u-sqmean)^2)/(2*sqvar)),u,tk(1),tk(1)+ ...
range*(p)/(2^bnum))/int(exp(((-1/3)*-(u-sqmean)^2)/(2*sqvar)),u,tk(1),max(max(square)));
rk(p)=(tk(p+1)+tk(p))/2;
grk(countr)=rk(p);
gtk(countt)=tk(p+1);
countt=countt+1;
countr=countr+1;
end
end

value=0;
for l=1:dim % this loop determines at which decision level
for m=1:dim % the normalized value lies and assigns a value
if square(l,m)<=tk(2)
value=0;
elseif square(l,m)>=tk(length(tk))
value=length(tk)-2;
else
for o=2:length(tk)-1
if square(l,m)>= tk(o)
value=o-1;
end
end
end
if j==1 % this loop places the values for the decision level into the full matrix
A(dimhold+l,dimhold+m)=value;
elseif j==2
A(dimhold+l,m)=value;
elseif j==3
A(l,dimhold+m)=value;
end
end
end
end
dimhold=dimhold+dim;
end