testbin.m

算术编码

% 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:  [e-mail]   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