1 #ifdef NEWSUBBAND
2 #include "common.h"
3 #include "enwindow.h"
4 #include "subband.h"
5
6 #ifdef REFERENCECODE
7 /************************************************************************
8 *
9 * window_subband()
10 *
11 * PURPOSE: Overlapping window on PCM samples
12 *
13 * SEMANTICS:
14 * 32 16-bit pcm samples are scaled to fractional 2's complement and
15 * concatenated to the end of the window buffer #x#. The updated window
16 * buffer #x# is then windowed by the analysis window #c# to produce the
17 * windowed sample #z#
18 *
19 ************************************************************************/
20 #define MAXCHANNELS 14 //6 "normal" channels, plus all the multi-linguals (7 ML channels?)
window_subband(double ** buffer,double z[64],int k)21 void window_subband (double **buffer, double z[64], int k)
22 {
23 typedef double XX[MAXCHANNELS][HAN_SIZE];
24 static XX *x;
25 double *xk;
26 int i;
27 static int off[MAXCHANNELS] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
28 static char init = 0;
29 double t;
30 double *ep0, *ep1, *ep2, *ep3, *ep4, *ep5, *ep6, *ep7;
31 if (!init) {
32 x = (XX *) mem_alloc (sizeof (XX), "x");
33 memset (x, 0, MAXCHANNELS * HAN_SIZE * sizeof (double));
34 init = 1;
35 }
36 xk = (*x)[k];
37
38 //fprintf(stdout,"k %i\n",k);
39
40 /* replace 32 oldest samples with 32 new samples */
41 for (i = 0; i < 32; i++)
42 xk[31 - i + off[k]] = *(*buffer)++ / SCALE;
43
44 ep0 = &enwindow[0];
45 ep1 = &enwindow[64];
46 ep2 = &enwindow[128];
47 ep3 = &enwindow[192];
48 ep4 = &enwindow[256];
49 ep5 = &enwindow[320];
50 ep6 = &enwindow[384];
51 ep7 = &enwindow[448];
52
53 /* shift samples into proper window positions */
54 for (i = 0; i < 64; i++) {
55 t = xk[(i + off[k]) & (512 - 1)] * *ep0++;
56 t += xk[(i + 64 + off[k]) & (512 - 1)] * *ep1++;
57 t += xk[(i + 128 + off[k]) & (512 - 1)] * *ep2++;
58 t += xk[(i + 192 + off[k]) & (512 - 1)] * *ep3++;
59 t += xk[(i + 256 + off[k]) & (512 - 1)] * *ep4++;
60 t += xk[(i + 320 + off[k]) & (512 - 1)] * *ep5++;
61 t += xk[(i + 384 + off[k]) & (512 - 1)] * *ep6++;
62 t += xk[(i + 448 + off[k]) & (512 - 1)] * *ep7++;
63 z[i] = t;
64 }
65
66 off[k] += 480; /*offset is modulo (HAN_SIZE-1) */
67 off[k] &= HAN_SIZE - 1;
68
69 }
70
71
72 /************************************************************************
73 *
74 * filter_subband()
75 *
76 * PURPOSE: Calculates the analysis filter bank coefficients
77 *
78 * SEMANTICS:
79 * The windowed samples #z# is filtered by the digital filter matrix #m#
80 * to produce the subband samples #s#. This done by first selectively
81 * picking out values from the windowed samples, and then multiplying
82 * them by the filter matrix, producing 32 subband samples.
83 *
84 ************************************************************************/
filter_subband(double z[HAN_SIZE],double s[SBLIMIT])85 void filter_subband (double z[HAN_SIZE], double s[SBLIMIT])
86 {
87 double yprime[32];
88 register int i, j;
89
90 static double m[16][32];
91 static int init = 0;
92
93 if (init == 0) {
94 init++;
95 create_dct_matrix (m);
96 }
97
98 yprime[0] = z[16];
99 for (i = 1; i <= 16; i++)
100 yprime[i] = z[i + 16] + z[16 - i];
101 for (i = 17; i <= 31; i++)
102 yprime[i] = z[i + 16] - z[80 - i];
103
104 for (i = 15; i >= 0; i--) {
105 register double s0 = 0.0, s1 = 0.0;
106 register double *mp = m[i];
107 register double *xinp = yprime;
108 for (j = 0; j < 8; j++) {
109 s0 += *mp++ * *xinp++;
110 s1 += *mp++ * *xinp++;
111 s0 += *mp++ * *xinp++;
112 s1 += *mp++ * *xinp++;
113 }
114 s[i] = s0 + s1;
115 s[31 - i] = s0 - s1;
116 }
117 }
118 #endif //REFERENCECODE
119
create_dct_matrix(double filter[16][32])120 void create_dct_matrix (double filter[16][32])
121 {
122 register int i, k;
123
124 for (i = 0; i < 16; i++)
125 for (k = 0; k < 32; k++) {
126 if ((filter[i][k] = 1e9 * cos ((double) ((2 * i + 1) * k * PI64))) >= 0)
127 modf (filter[i][k] + 0.5, &filter[i][k]);
128 else
129 modf (filter[i][k] - 0.5, &filter[i][k]);
130 filter[i][k] *= 1e-9;
131 }
132 }
133
134 //____________________________________________________________________________
135 //____ WindowFilterSubband() _________________________________________
136 //____ RS&A - Feb 2003 _______________________________________________________
WindowFilterSubband(short * pBuffer,int ch,double s[SBLIMIT])137 void WindowFilterSubband (short *pBuffer, int ch, double s[SBLIMIT])
138 {
139 register int i, j;
140 int pa, pb, pc, pd, pe, pf, pg, ph;
141 double t;
142 double *dp, *dp2;
143 double *pEnw;
144 double y[64];
145 double yprime[32];
146
147 static double x[2][512];
148 static double m[16][32];
149 static int init = 0;
150 static int off[2];
151 static int half[2];
152
153 if (init == 0) {
154 init++;
155 off[0] = 0;
156 off[1] = 0;
157 half[0] = 0;
158 half[1] = 0;
159 for (i = 0; i < 2; i++)
160 for (j = 0; j < 512; j++)
161 x[i][j] = 0;
162 create_dct_matrix (m);
163 }
164
165 dp = x[ch] + off[ch] + half[ch] * 256;
166
167 /* replace 32 oldest samples with 32 new samples */
168 for (i = 0; i < 32; i++)
169 dp[(31 - i) * 8] = (double) pBuffer[i] / SCALE;
170
171 // looks like "school example" but does faster ...
172 dp = (x[ch] + half[ch] * 256);
173 pa = off[ch];
174 pb = (pa + 1) % 8;
175 pc = (pa + 2) % 8;
176 pd = (pa + 3) % 8;
177 pe = (pa + 4) % 8;
178 pf = (pa + 5) % 8;
179 pg = (pa + 6) % 8;
180 ph = (pa + 7) % 8;
181
182 for (i = 0; i < 32; i++) {
183 dp2 = dp + i * 8;
184 pEnw = enwindow + i;
185 t = dp2[pa] * pEnw[0];
186 t += dp2[pb] * pEnw[64];
187 t += dp2[pc] * pEnw[128];
188 t += dp2[pd] * pEnw[192];
189 t += dp2[pe] * pEnw[256];
190 t += dp2[pf] * pEnw[320];
191 t += dp2[pg] * pEnw[384];
192 t += dp2[ph] * pEnw[448];
193 y[i] = t;
194 }
195
196 yprime[0] = y[16]; // Michael Chen�s dct filter
197
198 dp = half[ch] ? x[ch] : (x[ch] + 256);
199 pa = half[ch] ? (off[ch] + 1) & 7 : off[ch];
200 pb = (pa + 1) % 8;
201 pc = (pa + 2) % 8;
202 pd = (pa + 3) % 8;
203 pe = (pa + 4) % 8;
204 pf = (pa + 5) % 8;
205 pg = (pa + 6) % 8;
206 ph = (pa + 7) % 8;
207
208 for (i = 0; i < 32; i++) {
209 dp2 = dp + i * 8;
210 pEnw = enwindow + i + 32;
211 t = dp2[pa] * pEnw[0];
212 t += dp2[pb] * pEnw[64];
213 t += dp2[pc] * pEnw[128];
214 t += dp2[pd] * pEnw[192];
215 t += dp2[pe] * pEnw[256];
216 t += dp2[pf] * pEnw[320];
217 t += dp2[pg] * pEnw[384];
218 t += dp2[ph] * pEnw[448];
219 y[i + 32] = t;
220 // 1st pass on Michael Chen�s dct filter
221 if (i > 0 && i < 17)
222 yprime[i] = y[i + 16] + y[16 - i];
223 }
224
225 // 2nd pass on Michael Chen�s dct filter
226 for (i = 17; i < 32; i++)
227 yprime[i] = y[i + 16] - y[80 - i];
228
229 for (i = 15; i >= 0; i--) {
230 register double s0 = 0.0, s1 = 0.0;
231 register double *mp = m[i];
232 register double *xinp = yprime;
233 for (j = 0; j < 8; j++) {
234 s0 += *mp++ * *xinp++;
235 s1 += *mp++ * *xinp++;
236 s0 += *mp++ * *xinp++;
237 s1 += *mp++ * *xinp++;
238 }
239 s[i] = s0 + s1;
240 s[31 - i] = s0 - s1;
241 }
242
243 half[ch] = (half[ch] + 1) & 1;
244 if (half[ch] == 1)
245 off[ch] = (off[ch] + 7) & 7;
246 }
247
248
249 #endif //NEWSUBBAND
250