2 * Siren Encoder/Decoder library
4 * @author: Youness Alaoui <kakaroto@kakaroto.homelinux.net>
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.
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.
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.
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,
54 static int dct4_initialized = 0;
60 double scale_320 = (float) sqrt (2.0 / 320);
61 double scale_640 = (float) sqrt (2.0 / 640);
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));
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);
90 siren_dct4 (float *Source, float *Destination, int dct_length)
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];
100 float *In_Ptr = NULL;
101 float *In_Ptr_low = NULL;
102 float *In_Ptr_high = NULL;
105 float *Out_ptr_low = NULL;
106 float *Out_ptr_high = NULL;
107 float mult1, mult2, mult3, mult4, mult5, mult6, mult7, mult8, mult9, mult10;
110 if (dct4_initialized == 0)
113 if (dct_length == 640) {
115 dct_core = dct_core_640;
118 dct_core = dct_core_320;
121 Out_ptr = OutBuffer1;
122 NextOut_ptr = OutBuffer2;
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));
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);
137 Out_ptr = NextOut_ptr;
138 NextOut_ptr = In_Ptr;
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;
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--) {
165 for (j = 0; j < (1 << i); j++) {
166 dct_table_ptr = *dct_table_ptr_ptr;
168 Out_ptr_low = Destination + (j * (dct_length >> i));
170 Out_ptr_low = Out_ptr + (j * (dct_length >> i));
172 Out_ptr_high = Out_ptr_low + (dct_length >> i);
174 In_Ptr_low = In_Ptr + (j * (dct_length >> i));
175 In_Ptr_high = In_Ptr_low + (dct_length >> (i + 1));
178 (*In_Ptr_low * (*dct_table_ptr).cos) -
179 (*In_Ptr_high * (*dct_table_ptr).msin);
181 (*In_Ptr_high++ * (*dct_table_ptr).cos) +
182 (*In_Ptr_low++ * (*dct_table_ptr).msin);
185 (*In_Ptr_low * (*dct_table_ptr).cos) +
186 (*In_Ptr_high * (*dct_table_ptr).msin);
188 (*In_Ptr_low++ * (*dct_table_ptr).msin) -
189 (*In_Ptr_high++ * (*dct_table_ptr).cos);
191 } while (Out_ptr_low < Out_ptr_high);
195 Out_ptr = NextOut_ptr;
196 NextOut_ptr = In_Ptr;