gst/siren/dct4.c

   1 /*
   2  * Siren Encoder/Decoder library
   3  *
   4  *   @author: Youness Alaoui <kakaroto@kakaroto.homelinux.net>
   5  *
   6  * This library is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Library General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2 of the License, or (at your option) any later version.
  10  *
  11  * This library is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Library General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Library General Public
  17  * License along with this library; if not, write to the
  18  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19  * Boston, MA 02111-1307, USA.
  20  */
  21
  22
  23 #include "siren7.h"
  24
  25
  26 #define PI 3.1415926
  27
  28 typedef struct {
  29   float cos;
  30   float msin;
  31 } dct_table_type;
  32
  33 static float dct_core_320[100];
  34 static float dct_core_640[100];
  35 static dct_table_type dct_table_5[5];
  36 static dct_table_type dct_table_10[10];
  37 static dct_table_type dct_table_20[20];
  38 static dct_table_type dct_table_40[40];
  39 static dct_table_type dct_table_80[80];
  40 static dct_table_type dct_table_160[160];
  41 static dct_table_type dct_table_320[320];
  42 static dct_table_type dct_table_640[640];
  43 static dct_table_type *dct_tables[8] = {dct_table_5,
  44                                         dct_table_10,
  45                                         dct_table_20,
  46                                         dct_table_40,
  47                                         dct_table_80,
  48                                         dct_table_160,
  49                                         dct_table_320,
  50                                         dct_table_640};
  51
  52 static int dct4_initialized = 0;
  53
  54 void siren_dct4_init() {
  55   int i, j = 0;
  56   double scale_320 = (float) sqrt(2.0/320);
  57   double scale_640 = (float) sqrt(2.0/640);
  58   double angle;
  59   double scale;
  60
  61   /* set up dct4 tables */
  62   for(i = 0; i < 10; i++) {
  63     angle = (float) ((i + 0.5) *  PI);
  64     for (j = 0 ; j < 10; j++) {
  65       dct_core_320[(i*10)+j] = (float) (scale_320 * cos((j + 0.5) * angle / 10));
  66       dct_core_640[(i*10)+j] = (float) (scale_640 * cos((j + 0.5) * angle / 10));
  67     }
  68   }
  69
  70   for(i = 0; i < 8; i++) {
  71     scale = (float) (PI / ((5 << i) * 4));
  72     for (j = 0 ; j < (5 << i); j++) {
  73       angle = (float) (j + 0.5) * scale;
  74       dct_tables[i][j].cos = (float) cos(angle);
  75       dct_tables[i][j].msin = (float)  -sin(angle);
  76     }
  77   }
  78
  79   dct4_initialized = 1;
  80 }
  81
  82
  83 void siren_dct4(float *Source, float *Destination, int dct_length) {
  84   int log_length = 0;
  85   float * dct_core = NULL;
  86   dct_table_type ** dct_table_ptr_ptr = NULL;
  87   dct_table_type * dct_table_ptr = NULL;
  88   float OutBuffer1[640];
  89   float OutBuffer2[640];
  90   float *Out_ptr;
  91   float *NextOut_ptr;
  92   float *In_Ptr = NULL;
  93   float *In_Ptr_low = NULL;
  94   float *In_Ptr_high = NULL;
  95   float In_val_low;
  96   float In_val_high;
  97   float *Out_ptr_low = NULL;
  98   float *Out_ptr_high = NULL;
  99   float mult1, mult2, mult3, mult4, mult5, mult6, mult7, mult8, mult9, mult10;
 100   int i,j;
 101
 102   if (dct4_initialized == 0)
 103     siren_dct4_init();
 104
 105   if (dct_length == 640) {
 106     log_length = 5;
 107     dct_core = dct_core_640;
 108   } else {
 109     log_length = 4;
 110     dct_core = dct_core_320;
 111   }
 112
 113   Out_ptr = OutBuffer1;
 114   NextOut_ptr = OutBuffer2;
 115   In_Ptr = Source;
 116   for (i = 0; i <= log_length; i++) {
 117     for (j = 0; j < (1 << i); j++) {
 118       Out_ptr_low = Out_ptr + (j * (dct_length >> i));
 119       Out_ptr_high = Out_ptr + ( (j+1) * (dct_length >> i));
 120       do {
 121         In_val_low = *In_Ptr++;
 122         In_val_high = *In_Ptr++;
 123         *Out_ptr_low++ = In_val_low + In_val_high;
 124         *--Out_ptr_high = In_val_low - In_val_high;
 125       } while (Out_ptr_low < Out_ptr_high);
 126     }
 127
 128     In_Ptr = Out_ptr;
 129     Out_ptr = NextOut_ptr;
 130     NextOut_ptr = In_Ptr;
 131   }
 132
 133   for (i = 0; i < (2 << log_length); i++) {
 134     for (j = 0 ; j < 10 ;  j ++) {
 135       mult1 = In_Ptr[(i*10)] * dct_core[j*10];
 136       mult2 = In_Ptr[(i*10) + 1] * dct_core[(j*10) + 1];
 137       mult3 = In_Ptr[(i*10) + 2] * dct_core[(j*10) + 2];
 138       mult4 = In_Ptr[(i*10) + 3] * dct_core[(j*10) + 3];
 139       mult5 = In_Ptr[(i*10) + 4] * dct_core[(j*10) + 4];
 140       mult6 = In_Ptr[(i*10) + 5] * dct_core[(j*10) + 5];
 141       mult7 = In_Ptr[(i*10) + 6] * dct_core[(j*10) + 6];
 142       mult8 = In_Ptr[(i*10) + 7] * dct_core[(j*10) + 7];
 143       mult9 = In_Ptr[(i*10) + 8] * dct_core[(j*10) + 8];
 144       mult10 = In_Ptr[(i*10) + 9] * dct_core[(j*10) + 9];
 145       Out_ptr[(i*10)+j] = mult1 + mult2 + mult3 + mult4 +
 146           mult5 + mult6 + mult7 + mult8 +
 147           mult9 + mult10;
 148     }
 149   }
 150
 151
 152   In_Ptr = Out_ptr;
 153   Out_ptr = NextOut_ptr;
 154   NextOut_ptr = In_Ptr;
 155   dct_table_ptr_ptr = dct_tables;
 156   for (i = log_length; i >= 0; i--) {
 157     dct_table_ptr_ptr++;
 158     for (j = 0; j < (1 << i); j++) {
 159       dct_table_ptr = *dct_table_ptr_ptr;
 160       if ( i == 0 )
 161         Out_ptr_low = Destination + (j * (dct_length >> i));
 162       else
 163         Out_ptr_low = Out_ptr + (j * (dct_length >> i));
 164
 165       Out_ptr_high = Out_ptr_low + (dct_length >> i);
 166
 167       In_Ptr_low = In_Ptr + (j * (dct_length >> i));
 168       In_Ptr_high = In_Ptr_low + (dct_length >> (i+1));
 169       do {
 170         *Out_ptr_low++ = (*In_Ptr_low * (*dct_table_ptr).cos) - (*In_Ptr_high * (*dct_table_ptr).msin);
 171         *--Out_ptr_high = (*In_Ptr_high++ * (*dct_table_ptr).cos) + (*In_Ptr_low++ * (*dct_table_ptr).msin);
 172         dct_table_ptr++;
 173         *Out_ptr_low++ = (*In_Ptr_low * (*dct_table_ptr).cos) + (*In_Ptr_high * (*dct_table_ptr).msin);
 174         *--Out_ptr_high = (*In_Ptr_low++ * (*dct_table_ptr).msin) - (*In_Ptr_high++ * (*dct_table_ptr).cos);
 175         dct_table_ptr++;
 176       } while (Out_ptr_low < Out_ptr_high);
 177     }
 178
 179     In_Ptr = Out_ptr;
 180     Out_ptr = NextOut_ptr;
 181     NextOut_ptr = In_Ptr;
 182   }
 183
 184 }