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 {
  30   float cos;
  31   float msin;
  32 } dct_table_type;
  33
  34 static float dct_core_320[100];
  35 static float dct_core_640[100];
  36 static dct_table_type dct_table_5[5];
  37 static dct_table_type dct_table_10[10];
  38 static dct_table_type dct_table_20[20];
  39 static dct_table_type dct_table_40[40];
  40 static dct_table_type dct_table_80[80];
  41 static dct_table_type dct_table_160[160];
  42 static dct_table_type dct_table_320[320];
  43 static dct_table_type dct_table_640[640];
  44 static dct_table_type *dct_tables[8] = { dct_table_5,
  45   dct_table_10,
  46   dct_table_20,
  47   dct_table_40,
  48   dct_table_80,
  49   dct_table_160,
  50   dct_table_320,
  51   dct_table_640
  52 };
  53
  54 static int dct4_initialized = 0;
  55
  56 void
  57 siren_dct4_init ()
  58 {
  59   int i, j = 0;
  60   double scale_320 = (float) sqrt (2.0 / 320);
  61   double scale_640 = (float) sqrt (2.0 / 640);
  62   double angle;
  63   double scale;
  64
  65   /* set up dct4 tables */
  66   for (i = 0; i < 10; i++) {
  67     angle = (float) ((i + 0.5) * PI);
  68     for (j = 0; j < 10; j++) {
  69       dct_core_320[(i * 10) + j] =
  70           (float) (scale_320 * cos ((j + 0.5) * angle / 10));
  71       dct_core_640[(i * 10) + j] =
  72           (float) (scale_640 * cos ((j + 0.5) * angle / 10));
  73     }
  74   }
  75
  76   for (i = 0; i < 8; i++) {
  77     scale = (float) (PI / ((5 << i) * 4));
  78     for (j = 0; j < (5 << i); j++) {
  79       angle = (float) (j + 0.5) * scale;
  80       dct_tables[i][j].cos = (float) cos (angle);
  81       dct_tables[i][j].msin = (float) -sin (angle);
  82     }
  83   }
  84
  85   dct4_initialized = 1;
  86 }
  87
  88
  89 void
  90 siren_dct4 (float *Source, float *Destination, int dct_length)
  91 {
  92   int log_length = 0;
  93   float *dct_core = NULL;
  94   dct_table_type **dct_table_ptr_ptr = NULL;
  95   dct_table_type *dct_table_ptr = NULL;
  96   float OutBuffer1[640];
  97   float OutBuffer2[640];
  98   float *Out_ptr;
  99   float *NextOut_ptr;
 100   float *In_Ptr = NULL;
 101   float *In_Ptr_low = NULL;
 102   float *In_Ptr_high = NULL;
 103   float In_val_low;
 104   float In_val_high;
 105   float *Out_ptr_low = NULL;
 106   float *Out_ptr_high = NULL;
 107   float mult1, mult2, mult3, mult4, mult5, mult6, mult7, mult8, mult9, mult10;
 108   int i, j;
 109
 110   if (dct4_initialized == 0)
 111     siren_dct4_init ();
 112
 113   if (dct_length == 640) {
 114     log_length = 5;
 115     dct_core = dct_core_640;
 116   } else {
 117     log_length = 4;
 118     dct_core = dct_core_320;
 119   }
 120
 121   Out_ptr = OutBuffer1;
 122   NextOut_ptr = OutBuffer2;
 123   In_Ptr = Source;
 124   for (i = 0; i <= log_length; i++) {
 125     for (j = 0; j < (1 << i); j++) {
 126       Out_ptr_low = Out_ptr + (j * (dct_length >> i));
 127       Out_ptr_high = Out_ptr + ((j + 1) * (dct_length >> i));
 128       do {
 129         In_val_low = *In_Ptr++;
 130         In_val_high = *In_Ptr++;
 131         *Out_ptr_low++ = In_val_low + In_val_high;
 132         *--Out_ptr_high = In_val_low - In_val_high;
 133       } while (Out_ptr_low < Out_ptr_high);
 134     }
 135
 136     In_Ptr = Out_ptr;
 137     Out_ptr = NextOut_ptr;
 138     NextOut_ptr = In_Ptr;
 139   }
 140
 141   for (i = 0; i < (2 << log_length); i++) {
 142     for (j = 0; j < 10; j++) {
 143       mult1 = In_Ptr[(i * 10)] * dct_core[j * 10];
 144       mult2 = In_Ptr[(i * 10) + 1] * dct_core[(j * 10) + 1];
 145       mult3 = In_Ptr[(i * 10) + 2] * dct_core[(j * 10) + 2];
 146       mult4 = In_Ptr[(i * 10) + 3] * dct_core[(j * 10) + 3];
 147       mult5 = In_Ptr[(i * 10) + 4] * dct_core[(j * 10) + 4];
 148       mult6 = In_Ptr[(i * 10) + 5] * dct_core[(j * 10) + 5];
 149       mult7 = In_Ptr[(i * 10) + 6] * dct_core[(j * 10) + 6];
 150       mult8 = In_Ptr[(i * 10) + 7] * dct_core[(j * 10) + 7];
 151       mult9 = In_Ptr[(i * 10) + 8] * dct_core[(j * 10) + 8];
 152       mult10 = In_Ptr[(i * 10) + 9] * dct_core[(j * 10) + 9];
 153       Out_ptr[(i * 10) + j] = mult1 + mult2 + mult3 + mult4 +
 154           mult5 + mult6 + mult7 + mult8 + mult9 + mult10;
 155     }
 156   }
 157
 158
 159   In_Ptr = Out_ptr;
 160   Out_ptr = NextOut_ptr;
 161   NextOut_ptr = In_Ptr;
 162   dct_table_ptr_ptr = dct_tables;
 163   for (i = log_length; i >= 0; i--) {
 164     dct_table_ptr_ptr++;
 165     for (j = 0; j < (1 << i); j++) {
 166       dct_table_ptr = *dct_table_ptr_ptr;
 167       if (i == 0)
 168         Out_ptr_low = Destination + (j * (dct_length >> i));
 169       else
 170         Out_ptr_low = Out_ptr + (j * (dct_length >> i));
 171
 172       Out_ptr_high = Out_ptr_low + (dct_length >> i);
 173
 174       In_Ptr_low = In_Ptr + (j * (dct_length >> i));
 175       In_Ptr_high = In_Ptr_low + (dct_length >> (i + 1));
 176       do {
 177         *Out_ptr_low++ =
 178             (*In_Ptr_low * (*dct_table_ptr).cos) -
 179             (*In_Ptr_high * (*dct_table_ptr).msin);
 180         *--Out_ptr_high =
 181             (*In_Ptr_high++ * (*dct_table_ptr).cos) +
 182             (*In_Ptr_low++ * (*dct_table_ptr).msin);
 183         dct_table_ptr++;
 184         *Out_ptr_low++ =
 185             (*In_Ptr_low * (*dct_table_ptr).cos) +
 186             (*In_Ptr_high * (*dct_table_ptr).msin);
 187         *--Out_ptr_high =
 188             (*In_Ptr_low++ * (*dct_table_ptr).msin) -
 189             (*In_Ptr_high++ * (*dct_table_ptr).cos);
 190         dct_table_ptr++;
 191       } while (Out_ptr_low < Out_ptr_high);
 192     }
 193
 194     In_Ptr = Out_ptr;
 195     Out_ptr = NextOut_ptr;
 196     NextOut_ptr = In_Ptr;
 197   }
 198
 199 }