lbg2.c

full search vector quantizers with codebooks of sizes

/* This program designs full search vector quantizers with codebooks of 
   sizes 1,2,4,.....,m , block-length of k, and based on a training
   sequence of length l samples. It writes the generated codebooks in an
   output file whose name has to be supplied. To run the program use:
   fsvq <k> <m> <filename for training sequence> <l>
	<perturbation factor> <distortion threshold> <output filename>
   For a description of the algorithm see the following two papers:
   1. Y.Linde, A.Buzo, and R.M.Gray,"AN Algorithm for Vector Quantizer
      Design,"IEEE-COM,Vol.28,1/80,pp.84-95.
   2. R.M.Gray, and Y.Linde,"Vector Quantizers and Predictive Quantizers
      for Gauss-Markov Sources,"IEEE-COM,Vol.30,2/82,pp.381-389.
   Nader Moayeri ;  1/8/84                                                */
#include <math.h>
#include <stdio.h>
#include <sys/types.h>
#include <sys/times.h>
#define PMODE 0644
#define MINCON 8.5e+37
float x[4000000], code[4096], centroid[4096] ;
int count[2048] ;
main(argc,argv)
int argc ;
char *argv[] ;
{  register int i ;
   register int j ;
   register int kk ;
   register float *xp ;
   register float *codep ;
   int k, l, m, n, ii, jj, ll, mm, iterate, nv, nzero, fd, fd2, max ;
   int zero[10] ;
   float a, delta, pert, dis_th, d, dd, dis, min ;
   struct tms time_buf ;

   k = atoi(argv[1]);
   m = atoi(argv[2]);
   l = atoi(argv[4]);
   if ( ( fd2 = creat(argv[7],PMODE) )  ==  -1 )
      {  printf("Error in creating the output file.\n") ;
	 exit(1) ;
      }
   if ( ( i = write(fd2,&k,4) )  !=  4 )  exit(1) ;
   if ( ( i = write(fd2,&m,4) )  !=  4 )  exit(1) ;
   if ( ( i = write(fd2,&l,4) )  !=  4 )  exit(1) ;
   n = 1 ;
   delta = atof(argv[5]) ;
   dis_th = atof(argv[6]) ;
   pert = 1. + delta ;  
   printf("k=%d   l=%d   delta=%f   dis_th=%f\n\n",k,l,delta,dis_th) ;
   if ( ( fd = open(argv[3],0) )  ==  -1 )
      {  printf("Error in opening the file of training sequence.\n") ;
	 exit(1) ;
      }
   lseek (fd, 400L, 0);
   if ( ( i = read(fd,x,4*l) )  !=  (4*l) )
      {  printf("Error in reading the training sequence.\n") ;
	 exit(1) ;
      }
   nv = l / k ;
   codep = &code[0];
   for (i = 0; i < k; i++);
        (*(codep++)) = 0;
   xp = &x[0];
   for (i = 0; i < nv; i++){
        codep = &code[0];
        for(j = 0; j < k; j++)
                (*(codep++)) += ((*(xp++))/nv);
   }
  

	 while ( iterate == 1 )
	    {  ++mm ;
               ii = 0 ;
               for ( i = 0; i < n; ++i )
                  {  count[i] = 0 ;
	             codep = &centroid[ii] ;
	             for ( j = 0; j < k; ++codep,++j )
	                (*codep) = 0. ;
	             ii += k ;
                  }
               ii = 0 ;
               dd = 0. ; /* per letter distortion of a transitory code */
               for ( i = 0; i < nv; ++i )
                  {  min = MINCON ;
                     jj = 0 ;
	             for ( j = 0; j < n; ++j )
	                {  dis = 0. ; /* per letter distortion between
	        	                 training vector i and codeword j */
	                   xp = &x[ii] ;
	                   codep = &code[jj] ;
	                   for ( kk = 0; kk < k; ++xp,++codep,++kk )
		              {  a = (*xp) - (*codep) ;
		                 dis += a*a ;
		              }
	                   if ( dis < min )
		              {  min = dis ;
		                 ll = j ; /* index for the closest codeword
		 	                     to training vector i          */
		              }
	                   jj += k ;
	                }
	             ++count[ll] ;
	             jj = k*ll ;
	             xp = &x[ii] ;
	             codep = &centroid[jj] ;
	             for ( j = 0; j < k; ++xp,++codep,++j )
	                (*codep) += (*xp) ;
	             ii += k ;
	             dd += min ;
                  }
	       dd /= (k*nv) ;
	       printf("# of iterations=%d   distortion=%f\n",mm,dd) ;
	       fflush(stdout) ;
	       if (  ((d-dd)/dd)  <  dis_th  )
		  {  iterate = 0 ;
                     printf("n=%d   distortion=%f\n",n,dd) ;
		     ii = 0 ;
                     for ( i = 0; i < n; ++i )
			{  codep = &code[ii] ;
			   for ( j = 0; j < k; ++codep,++j )
			      {  a = (*codep) ;
                                 printf("%f  ",a) ;
				 if ( ( jj = write(fd2,&a,4) )  !=  4 )
				    {  printf("Error in writing to the ") ;
				       printf("output file.\n") ;
				       exit(1) ;
                                    }
			      }
			   ii += k ;
			   printf("\n") ;
			}
                     printf("\n\n\n") ;
		  }
	       else
		  {  d = dd ;
		     nzero = 0 ; /* number of decision regions without
				    any training vector fallen into them */
		     for ( i = 0; i < n; ++i )
			if ( count[i] == 0 )
			   {  zero[nzero] = i ;
			      ++nzero ;
			   }
			else
			   {  ii = i * k ;
			      codep = &centroid[ii] ;
			      for ( j = 0; j < k; ++codep,++j )
				 (*codep) /= count[i] ;
			   }
		     if ( nzero > m/10+1 ) /*(dfl)*/
			{  printf("n=%d   nzero=%d\n",n,nzero) ;
	                   fflush(stdout) ;
			   for ( i = 0; i < nzero; ++i )
			      {  max = 0 ;
				 for ( j = 0; j < n; ++j )
				    if ( count[j] > max )
				       {  max = count[j] ;
					  ii = j ;
				       }
				 count[ii] = 0 ;
				 ii *= k ;
				 jj = k*zero[i] ;
				 xp = &centroid[ii] ;
				 codep = &centroid[jj] ;
				 for ( j = 0; j < k; ++xp,++codep,++j )
				    {  (*codep) = pert*(*xp) ;
				       (*xp) /= pert ;
				    }
			      }
			}
		     ii = 0 ;
		     for ( i = 0; i < n; ++i )
			{  xp = &centroid[ii] ;
			   codep = &code[ii] ;
			   for ( j = 0; j < k; ++xp,++codep,++j )
			      (*codep) = (*xp) ;
			   ii += k ;
			}
		  }
	    }
      }
   times(&time_buf) ;
   printf("%s   %s     user-time=%f   system-time=%f\n",
	  argv[3],argv[7],time_buf.tms_utime/60.,time_buf.tms_stime/60.) ;
   printf("******************************************************\n\n\n") ;   
}