% TestBin Test coding of binary sequence
% We have a sequence of binary symbols (0/1), where the
% probability of a 1 is p1 and the probability of a 0 is p0=1-p1
% Direct arithmetic coding (DAC) of this sequence (p1 is given) is compared to
% [adaptive arithmetic coding (AAC) and] count down arithmetic coding (CDAC).
% The value we are interested in is the bit rate, which will be a random
% variable depending on the coding method (DAC/AAC/CDAC),
% the length of the sequence (N), and the probability of a 1 (p1).
% The expectation (and sometimes the distribution) for this random variable
% is interesting.
%----------------------------------------------------------------------
% Copyright (c) 2000. Karl Skretting. All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail: karl.skretting@tn.his.no Homepage: http://www.ux.his.no/~karlsk/
%
% HISTORY:
% Ver. 1.0 14.05.2001 KS: function made
%----------------------------------------------------------------------
clear all;
p1=0.2;
p0=1-p1;
N=100;
% All sequences are possible. N1 is the number of ones, N0 the number of zeros,
% in the sequence of length N. We have N1+N0=N
% we make an example sequence
rand('state',301);
x=floor(rand(N,1)+p1);
I=find(x);
N1=length(I)
N0=N-N1;
return;
% some different code segments and more or less meaningful comments
% For DAC (p1 is given) the expected bit rate is given by the entropy
ent=-(p0*log2(p0)+p1*log2(p1));
if 1
% the entropy (or expected bit rate for DAC) is a function of p1
p=linspace(0.001,0.999,159);
entp=-((1-p).*log2((1-p))+p.*log2(p));
figure(1);
plot(p,entp);
title('The entropy function');
xlabel('Probability');
ylabel('Entropy');
grid on;
end
% but the bit rate for an actual sequence will depend on the number
% of zeros and ones in the sequence.
% The probability that N1=n depends on p1 and N, P{N1==n}=binopdf(n,N,p1)
% The bit rate given N1 (and N) and p1 is
bits=N1*log2(1/p1)+N0*log2(1/p0);
br=bits/N;
% the pdf for br may be calculated using the binomial probability function
[m,v]=binostat(N,p1);
% we call this punction brpdf, this function will be zero
% except at some (N+1) distict points, these points are
brpdfx=((0:N)*log2(1/p1)+(N:(-1):0)*log2(1/p0))/N;
% we note that these are evenly distributed between min and max
% the step size is: log2((1-p1)/p1)/N and range is log2((1-p1)/p1)
brpdfx=linspace(log2(1/p0),log2(1/p1),N+1);
% and the distribution is the binomial one
brpdf=binopdf(0:N,N,p1);
plot(brpdfx,brpdf);
% form the graph we see that for p1=0.2 the bit rate will probably be between 0.5 and 1
% the mean of this function is the same as the entropy
% the plot for the entropy function may be extended to also display the
% confidence interval, the range which the actual value is more
% than pc*100% likely to be whitin. Note that this depends on N
if 0
N=102;
pc=0.98;
pcd=(1-pc)/2; % equal on both sides
p=linspace(0.001,0.999,159);
entp=-((1-p).*log2((1-p))+p.*log2(p));
e1=p;e2=p;
if N>100
% approximate by normal distribution
t=norminv(pcd,0,1);
% for p(i), the number of ones will be between e1(i) and e2(i)
e1=(N*p+t*sqrt(N*(p.*(1-p)))); % could be floor ?
e2=(N*p-t*sqrt(N*(p.*(1-p)))); % could be ceil ?
I=find(e1<0); e1(I)=0;
I=find(e2>N); e2(I)=N;
else
for i=1:length(p)
t=binoinv(pcd,N,p(i));
t2=binocdf(t,N,p(i));
t=max([0,t-1]);
t1=binocdf(t,N,p(i));
if (t2>t1)
e1(i)=t+(pcd-t1)/(t2-t1);
else
e1(i)=t;
end
end
end
e1=-(e1.*log2(p)+(N-e1).*log2(1-p))/N;
e2=fliplr(e1);
figure(1);
plot(p,entp,p,e1,p,e2);
title(['The entropy function with ',num2str(pc*100),' percent intervall (N=',int2str(N),')']);
xlabel('Probability');
ylabel('Entropy');
grid on;
end
% The actual case will usually be that we have a sequence to be coded,
% and we do not know the probability for ones (p1) in the generating process,
% we do not even know if the elements are independent nor if they are
% equally distributed but this latter we assume now.
% A method to code the sequence is to first code the number of ones, as
% before we assume that the length of the sequence N is known, then
% to code N1, 0<= N1 <= N, will need log2(N+1) bits (we simplify and assume
% that all values of N1 have the same probability, which is not true but practical).
% For the first element to be coded we then have P{x(1)==1}=N1/N
% First make a random sequence
N=100; % N is a known value
p1=0.2; % p1 and p0 is (unknown) parameters just used to make the random
p0=1-p1; % sequence x
rand('state',301);
x=floor(rand(N,1)+p1); % x is a random sequence, x(n) is a random variable
I=find(x);
N1=length(I); % N1 is a random variable (but it will be given in the beginning)
N0=N-N1; % so is N0
% Let us define some random variables, M(n) n=1:N, where M(n) is number of
% ones in x(1:n), obviously we will have M(n) <= M(n+1) and M(N)=N1.
% We define an array PM of size Nx(N1+1) where PM(n,n1)=probability{M(n)==(n1-1)}
% We have PM(n,n1) = P{ M(n)==(n1-1) }
% PM(1,1) = P{ M(1)==0 } = N0/N = 1-N1/N
% PM(1,2) = P{ M(1)==1 } = N1/N
% PM(1,k) = P{ M(1)==(k-1) } = 0 for k>2
% We note that we should have sum(PM(n,:))==1 for n=1:N
% PM(n,n1+1) = P{ M(n)==n1 }
% = P{ M(n-1)==n1 } * P{ x(n)==0 | M(n-1)==n1) }
% + P{ M(n-1)==(n1-1) } * P{ x(n)==1 | M(n-1)==(n1-1) }
% = PM(n-1,n1+1) * (1-(N1-n1)/(N-n+1))
% + PM(n-1,n1) * (N1-n1+1)/(N-n+1)
PM=zeros(N,N1+1);
PM(1,1)=1-N1/N;
PM(1,2)=N1/N;
for n=2:N
PM(n,1)=PM(n-1,1)*((N-N1-n+1)/(N-n+1));
for k=2:(N1+1)
PM(n,k)=PM(n-1,k)*((N-N1-n+k)/(N-n+1))+PM(n-1,k-1)*(N1-k+2)/(N-n+1);
end
end
% how many bits will we need to code the sequence x
N=100; % N is a known value
p1=0.2; % p1 and p0 is (unknown) parameters just used to make the random
p0=1-p1; % sequence x
btot=0;
LC=200;
b1=zeros(LC,1);
b2=zeros(LC,1);
for count=1:LC
rand('state',301+count);
x=floor(rand(N,1)+p1); % x is a random sequence, x(n) is a random variable
I=find(x);
N1=length(I); % N1 is a random variable (but it will be given in the beginning)
N0=N-N1; % so is N0
bits=log2(N+1); % to store N1
onesleft=N1;
totleft=N;
for n=1:N
if x(n)
bits=bits-log2(onesleft/totleft);
onesleft=onesleft-1;
else
bits=bits-log2(1-onesleft/totleft);
end
totleft=totleft-1;
end
% disp(['bits needed = ',num2str(bits)]);
btot=btot+bits;
b1(count)=bits/N;
b2(count)=-(N1*log2(p1)+(N-N1)*log2(1-p1))/N;
end
disp(['average bits needed = ',num2str(btot/LC)]);
mean([b1,b2]) % b2 has lower mean
var([b1,b2]) % b1 has often lower variance, for short sequences
max([b1,b2]) % usually max is lowest for b1
N=100; % N is a known value
p1=0.1; % p1 and p0 is (unknown) parameters just used to make the random
p0=1-p1; % sequence x
btot=0;
LC=20;
b1=zeros(LC,1);
b2=zeros(LC,1);
rand('state',300);
x=floor(rand(N,1)+p1); % x is a random sequence, x(n) is a random variable
I=find(x);
N1=length(I); % N1 is a random variable (but it will be given in the beginning)
N0=N-N1; % so is N0
% x=[ones(1,N1),zeros(1,N0)];
for count=1:LC
temp=randn(N,1);
[temp,I]=sort(temp);
x=x(I); % shuffle the sequence
bits=log2(N+1); % to store N1
onesleft=N1;
totleft=N;
for n=1:N
if x(n)
bits=bits-log2(onesleft/totleft);
onesleft=onesleft-1;
else
bits=bits-log2(1-onesleft/totleft);
end
totleft=totleft-1;
end
disp(['bits needed = ',num2str(bits)]);
btot=btot+bits;
b1(count)=bits/N;
b2(count)=(log2(N+1)-N1*log2(N1/N)-(N-N1)*log2((N-N1)/N))/N;
end
disp(['average bits needed = ',num2str(btot/LC)]);
mean(b1)
var(b1)
max(b1)
% this shows that for a sequence where length N, and number of ones N1 is
% given the number of bits to code this is predictable
% We set up B1(N,N1) as the number of bits needed to code a binary sequence
% of length N where there is N1 ones. (we code as above, adapting the
% probabilities as we go along, which saves ~3 bits)
% Since the order of the elements in x do not matter we put the ones first
N=10;
N1=1;
N1=min(N1,N-N1); % 0 <= N1 <= N/2
n=0:(N1-1);
B1=log2(N+1)+sum(log2(N-n))-sum(log2(N1-n));
B1=B1/N; % bit rate
disp(['B1(',int2str(N),',',int2str(N1),')=',num2str(B1)]);
% while using a fixed probability N1/N will give
B0=(log2(N+1)-N1*log2(N1)-(N-N1)*log2(N-N1))/N+log2(N);
disp(['B0(',int2str(N),',',int2str(N1),')=',num2str(B0)]);
disp(['The difference is ',num2str(N*(B0-B1))]);
% plot of bits saved by using adaptive probability
N=1000
if 1
if (N<101)
NN=1:ceil(N/2);
else
NN=ceil(linspace(1,N/2,25));
end
B1=zeros(size(NN));
B2=zeros(size(NN));
for i=1:length(NN)
N1=NN(i);
n=0:(N1-1);
B1(i)=log2(N+1)+sum(log2(N-n))-sum(log2(N1-n)); % bits
B2(i)=log2(N+1)-N1*log2(N1)-(N-N1)*log2(N-N1)+N*log2(N); % bits
end
figure(1);
plot(log2(NN),B2-B1);
title(['Bits saved by adapting probability, N=',int2str(N)]);
xlabel('Number of ones');
ylabel('Bits');
grid on;
end
p1=0.5;
if 1
NN=[10:10:100,200:100:1000];
B1=zeros(size(NN));
B2=zeros(size(NN));
for i=1:length(NN)
N=NN(i);
N1=floor(NN(i)*p1);
n=0:(N1-1);
B1(i)=log2(N+1)+sum(log2(N-n))-sum(log2(N1-n)); % bits
B2(i)=log2(N+1)-N1*log2(N1)-(N-N1)*log2(N-N1)+N*log2(N); % bits
end
figure(1);
plot(log2(NN),B2-B1);
title(['Bits saved by adapting probability, p1=',num2str(p1)]);
xlabel('log2(N)');
ylabel('Bits');
grid on;
end
% test of approximation
for N=[50,100,500,1000,5000,10000];
p=1/2;
N1=floor(p*N);
% first the correct answer
B1=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1)))); % bits
% then the approximation
N0=N-N1;
% The teoretic approximation
% b1est=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)+s2+...
% (1/(N+1/2)+(N+1/2)/(N+1/4)-(N0+1/2)/(N0+1/4)-(N1+1/2)/(N1+1/4)-10)*log2(exp(1));
disp(['N=',int2str(N),', N1=',int2str(N1)]);
b1est=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-0.010631*log2(log2(N))-1.28817;
disp([' 1000*Difference 1 is ',num2str(1000*(B1-b1est))]);
b2est=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-1.3226;
disp([' 1000*Difference 2 is ',num2str(1000*(B1-b2est))]);
end
% both approximations are well useful
% try to find "optimal coefficients
A=zeros(0,1);b=zeros(0,1);
for N=[50,100,250,500,1000,5000,10000];
for p=[0.10,0.2,0.3,0.4,0.45,0.5];
N1=floor(p*N);
b1=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1)))); % bits
b2=N*log2(N)-(N-N1)*log2(N-N1)-N1*log2(N1);
b3=1.5*log2(N)-0.5*log2(N-N1)-0.5*log2(N1);
% A=[A;[log2(N),log2(N-N1),log2(N1),1]];
A=[A;1];
b=[b;b1-b2-b3];
end
end
x=A\b;
% log2(N), log2(N-N1), log2(N1) 1
% x= 1.48892478527407 -0.49461589408971 -0.49616311290196 -1.29382354586074
% norm(b-A*x) = 0.02512
% log2(N), log2(N-N1), log2(N1) 1
% x= [ 1.5; -0.5; -0.5; -1.32113849612329 ];
% norm(b-A*x) = 0.04150
% log2(N), log2(N-N1), log2(N1), log2(log2(N)), 1
% x= 1.498973 -0.494827 -0.49623245 -0.0616557 -1.190184
% norm(b-A*x) = 0.01417
% we should not count the smallest N1
A=zeros(0,2);b=zeros(0,1);
for N=[50,70,100,250,500,750,1000];
% for p=[0.1,0.2,0.3,0.4,0.45,0.5];
for p=0.05:0.05:0.5;
N1=floor(p*N);
if N1>20
b1=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1)))); % bits
b2=N*log2(N)-(N-N1)*log2(N-N1)-N1*log2(N1);
b3=1.5*log2(N)-0.5*log2(N-N1)-0.5*log2(N1);
A=[A;[log2(log2(N)),1]];
b=[b;b1-b2-b3];
end
end
end
x=A\b;
% log2(N), log2(N-N1), log2(N1) 1
% x= [ 1.5; -0.5; -0.5; -1.32234151562252 ];
% norm(b-A*x) = 0.03216
% x= [ 1.5; -0.5; -0.5; -1.32298392950550 ];
% norm(b-A*x) = 0.03547 % using more p-values
% x= [ 1.50283174633830 -0.504191846468 -0.49998398295 -1.31180783192693 ];
% norm(b-A*x) = 0.02252 % using more p-values
% log2(N), log2(N-N1), log2(N1), log2(log2(N)), 1
% x= 1.5; -0.5; -0.5; -0.010631 -1.28817076478349
% norm(b-A*x) = 0.02134 % using more p-values
% x= 1.509169 -0.50062 -0.49959 -0.069397 -1.18453873231815
% norm(b-A*x) = 0.00757 % using more p-values
% we should not count the smallest N1
A=zeros(0,5);b=zeros(0,1);
for N=[50,100,250,500,1000,5000,10000];
% for p=[0.1,0.2,0.3,0.4,0.45,0.5];
for p=0.05:0.05:0.5;
N1=floor(p*N);
if N1>20
b1=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1)))); % bits
b2=N*log2(N)-(N-N1)*log2(N-N1)-N1*log2(N1);
b3=1.5*log2(N)-0.5*log2(N-N1)-0.5*log2(N1);
A=[A;[log2(N),log2(N-N1),log2(N1),log2(log2(N)),1]];
% A=[A;[log2(log2(N)),1]];
b=[b;b1-b2];
end
end
end
x=A\b;
% log2(N), log2(N-N1), log2(N1) 1
% x= [ 1.5; -0.5; -0.5; -1.32234151562252 ];
% norm(b-A*x) = 0.03216
% x= [ 1.5; -0.5; -0.5; -1.32298392950550 ];
% norm(b-A*x) = 0.03547 % using more p-values
% x= [ 1.50283174633830 -0.504191846468 -0.49998398295 -1.31180783192693 ];
% norm(b-A*x) = 0.02252 % using more p-values
% log2(N), log2(N-N1), log2(N1), log2(log2(N)), 1
% x= 1.5; -0.5; -0.5; -0.010631 -1.28817076478349
% norm(b-A*x) = 0.02134 % using more p-values
% x= 1.509169 -0.50062 -0.49959 -0.069397 -1.18453873231815
% norm(b-A*x) = 0.00757 % using more p-values
% the very good approximation function could then be
function b=BitEst(N,N1);
if (N1>(N/2)); N1=N-N1; end;
N0=N-N1;
if (N>1000)
b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-1.3256;
elseif (N1>20)
b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-0.020984*log2(log2(N))-1.25708;
else
b=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1))));
end
return;
% test this
for N=[50,100,500,1000,5000,10000];
p=1/2;
N1=floor(p*N);
% first the correct answer
B1=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1)))); % bits
% then the approximation
if (N1>(N/2)); N1=N-N1; end;
N0=N-N1;
if (N>1000)
b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-1.3256;
elseif (N1>20)
b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-0.020984*log2(log2(N))-1.25708;
else
b=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1))));
end
% display results
disp(['N=',int2str(N),', N1=',int2str(N1),' bits=',num2str(B1),...
' approx.bits=',num2str(b)]);
disp([' 1000*Difference is ',num2str(1000*(B1-b))]);
end
% --------------------------------------------
% Test how a sequence could be split to make more compact (and slower) compression
clear all;
p1=0.42;
p0=1-p1;
ent=-(p0*log2(p0)+p1*log2(p1));
N=300;
% All sequences are possible. N1 is the number of ones, N0 the number of zeros,
% in the sequence of length N. We have N1+N0=N
% we make an example sequence
if 0
rand('state',301);
x=floor(rand(N,1)+p1);
% x=[x';x';x']; % make this not so much random sequence
x=x(:);
else
x1=abs(floor(filter(1,[1,-0.97],randn(1,N))+0.5)); % an AR-1 signal
x=[bitget(x1,4);bitget(x1,3);bitget(x1,2);bitget(x1,1)]; % split has no effect
x=[bitget(x1,4),bitget(x1,3),bitget(x1,2),bitget(x1,1)]; % split has effect
x=x(:);
end
% a sequence, x, of length N should be coded
N=length(x);
I=find(x);
N1=length(I);
if (N1>(N/2)); N1=N-N1; end; % then N1 is number of zeros
N0=N-N1;
b=N1*log2(N/N1)+N0*log2(N/N0);
disp(['N=',int2str(N),', N1=',int2str(N1),', (e)bits=',num2str(b)]);
% alternative 1, direct
if (N>1000)
b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-1.3256;
elseif (N1>20)
b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-0.020984*log2(log2(N))-1.25708;
else
b=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1))));
end
disp(['N=',int2str(N),', N1=',int2str(N1),' N1/N=',num2str(N1/N),' bits=',num2str(b)]);
bits1=b;
% alternativ 2, splitt in 2 sequences
I=find(x(1:(N-1)));
J=find(x(1:(N-1))==0);
N11=length(I);
% store N11, 0 <= N11 <= (N-1)
bits2=log2(N); % bits needed to store N11
if N11<(N/2)
x1=x(I+1);
x2=[x(1);x(J+1)];
else
x1=[x(1);x(I+1)];
x2=x(J+1);
end
Ntemp=N;N1temp=N1;N=length(x1);N1=length(find(x1));N1=min([N1,N-N1]);N0=N-N1;
if (N>1000)
b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-1.3256;
elseif (N1>20)
b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-0.020984*log2(log2(N))-1.25708;
else
b=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1))));
end
disp(['N=',int2str(N),', N1=',int2str(N1),' N1/N=',num2str(N1/N),' bits=',num2str(b)]);
N=Ntemp;N1=N1temp;
bits21=b;
Ntemp=N;N1temp=N1;N=length(x2);N1=length(find(x2));N1=min([N1,N-N1]);N0=N-N1;
if (N>1000)
b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-1.3256;
elseif (N1>20)
b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-0.020984*log2(log2(N))-1.25708;
else
b=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1))));
end
disp(['N=',int2str(N),', N1=',int2str(N1),' N1/N=',num2str(N1/N),' bits=',num2str(b)]);
N=Ntemp;N1=N1temp;
bits22=b;
bits2=bits2+bits21+bits22;
% we must also indicate which method used
bits1=bits1+0.415; % probability 3/4
bits2=bits2+2; % probability 1/4
disp(['N=',int2str(N),', N1=',int2str(N1),', bits1=',num2str(bits1),...
', bits2=',num2str(bits2)]);
% This demonstrate some of the features of this coding technique
% it is only effective if the dependencies are to the previous symbol
% and not to the symbol 2 or more positions before.
% Also all possible dependicies where position has information will
% not be resolved.
% It is an effective way of coding completly random binary sequences (where p is
% the same for all symbols and no depenencies).
% Also if the dependency is to the preceeding bit this is effective coding.
% Note that if we force a split in several levels this is not the same as
% splitting for previous symbols '00', '01', '10' and '11'.
% For effective coding, it will always be important to resolve as much as possible
% of the dependencies early in the coding process, i.e. when
% we have most (general) knowledge of the symbols (and any dependencies).
