/src/vorbis/lib/mdct.c

Line	Count	Source (jump to first uncovered line)
1		/********************************************************************
2		* *
3		* THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
4		* USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
5		* GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
6		* IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
7		* *
8		* THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 *
9		* by the Xiph.Org Foundation https://xiph.org/ *
10		* *
11		********************************************************************
12
13		function: normalized modified discrete cosine transform
14		power of two length transform only [64 <= n ]
15
16		Original algorithm adapted long ago from _The use of multirate filter
17		banks for coding of high quality digital audio_, by T. Sporer,
18		K. Brandenburg and B. Edler, collection of the European Signal
19		Processing Conference (EUSIPCO), Amsterdam, June 1992, Vol.1, pp
20		211-214
21
22		The below code implements an algorithm that no longer looks much like
23		that presented in the paper, but the basic structure remains if you
24		dig deep enough to see it.
25
26		This module DOES NOT INCLUDE code to generate/apply the window
27		function. Everybody has their own weird favorite including me... I
28		happen to like the properties of y=sin(.5PI*sin^2(x)), but others may
29		vehemently disagree.
30
31		********************************************************************/
32
33		/* this can also be run as an integer transform by uncommenting a
34		define in mdct.h; the integerization is a first pass and although
35		it's likely stable for Vorbis, the dynamic range is constrained and
36		roundoff isn't done (so it's noisy). Consider it functional, but
37		only a starting point. There's no point on a machine with an FPU */
38
39		#include <stdio.h>
40		#include <stdlib.h>
41		#include <string.h>
42		#include <math.h>
43		#include "vorbis/codec.h"
44		#include "mdct.h"
45		#include "os.h"
46		#include "misc.h"
47
48		/* build lookups for trig functions; also pre-figure scaling and
49		some window function algebra. */
50
51	0	void mdct_init(mdct_lookup *lookup,int n){
52	0	int bitrev=_ogg_malloc(sizeof(bitrev)*(n/4));
53	0	DATA_TYPE T=_ogg_malloc(sizeof(T)*(n+n/4));
54
55	0	int i;
56	0	int n2=n>>1;
57	0	int log2n=lookup->log2n=rint(log((float)n)/log(2.f));
58	0	lookup->n=n;
59	0	lookup->trig=T;
60	0	lookup->bitrev=bitrev;
61
62		/* trig lookups... */
63
64	0	for(i=0;i<n/4;i++){
65	0	T[i2]=FLOAT_CONV(cos((M_PI/n)(4*i)));
66	0	T[i2+1]=FLOAT_CONV(-sin((M_PI/n)(4*i)));
67	0	T[n2+i2]=FLOAT_CONV(cos((M_PI/(2n))(2i+1)));
68	0	T[n2+i2+1]=FLOAT_CONV(sin((M_PI/(2n))(2i+1)));
69	0	}
70	0	for(i=0;i<n/8;i++){
71	0	T[n+i2]=FLOAT_CONV(cos((M_PI/n)(4i+2)).5);
72	0	T[n+i2+1]=FLOAT_CONV(-sin((M_PI/n)(4i+2)).5);
73	0	}
74
75		/* bitreverse lookup... */
76
77	0	{
78	0	int mask=(1<<(log2n-1))-1,i,j;
79	0	int msb=1<<(log2n-2);
80	0	for(i=0;i<n/8;i++){
81	0	int acc=0;
82	0	for(j=0;msb>>j;j++)
83	0	if((msb>>j)&i)acc\|=1<<j;
84	0	bitrev[i*2]=((~acc)&mask)-1;
85	0	bitrev[i*2+1]=acc;
86
87	0	}
88	0	}
89	0	lookup->scale=FLOAT_CONV(4.f/n);
90	0	}
91
92		/* 8 point butterfly (in place, 4 register) */
93	0	STIN void mdct_butterfly_8(DATA_TYPE *x){
94	0	REG_TYPE r0 = x[6] + x[2];
95	0	REG_TYPE r1 = x[6] - x[2];
96	0	REG_TYPE r2 = x[4] + x[0];
97	0	REG_TYPE r3 = x[4] - x[0];
98
99	0	x[6] = r0 + r2;
100	0	x[4] = r0 - r2;
101
102	0	r0 = x[5] - x[1];
103	0	r2 = x[7] - x[3];
104	0	x[0] = r1 + r0;
105	0	x[2] = r1 - r0;
106
107	0	r0 = x[5] + x[1];
108	0	r1 = x[7] + x[3];
109	0	x[3] = r2 + r3;
110	0	x[1] = r2 - r3;
111	0	x[7] = r1 + r0;
112	0	x[5] = r1 - r0;
113
114	0	}
115
116		/* 16 point butterfly (in place, 4 register) */
117	0	STIN void mdct_butterfly_16(DATA_TYPE *x){
118	0	REG_TYPE r0 = x[1] - x[9];
119	0	REG_TYPE r1 = x[0] - x[8];
120
121	0	x[8] += x[0];
122	0	x[9] += x[1];
123	0	x[0] = MULT_NORM((r0 + r1) * cPI2_8);
124	0	x[1] = MULT_NORM((r0 - r1) * cPI2_8);
125
126	0	r0 = x[3] - x[11];
127	0	r1 = x[10] - x[2];
128	0	x[10] += x[2];
129	0	x[11] += x[3];
130	0	x[2] = r0;
131	0	x[3] = r1;
132
133	0	r0 = x[12] - x[4];
134	0	r1 = x[13] - x[5];
135	0	x[12] += x[4];
136	0	x[13] += x[5];
137	0	x[4] = MULT_NORM((r0 - r1) * cPI2_8);
138	0	x[5] = MULT_NORM((r0 + r1) * cPI2_8);
139
140	0	r0 = x[14] - x[6];
141	0	r1 = x[15] - x[7];
142	0	x[14] += x[6];
143	0	x[15] += x[7];
144	0	x[6] = r0;
145	0	x[7] = r1;
146
147	0	mdct_butterfly_8(x);
148	0	mdct_butterfly_8(x+8);
149	0	}
150
151		/* 32 point butterfly (in place, 4 register) */
152	0	STIN void mdct_butterfly_32(DATA_TYPE *x){
153	0	REG_TYPE r0 = x[30] - x[14];
154	0	REG_TYPE r1 = x[31] - x[15];
155
156	0	x[30] += x[14];
157	0	x[31] += x[15];
158	0	x[14] = r0;
159	0	x[15] = r1;
160
161	0	r0 = x[28] - x[12];
162	0	r1 = x[29] - x[13];
163	0	x[28] += x[12];
164	0	x[29] += x[13];
165	0	x[12] = MULT_NORM( r0 * cPI1_8 - r1 * cPI3_8 );
166	0	x[13] = MULT_NORM( r0 * cPI3_8 + r1 * cPI1_8 );
167
168	0	r0 = x[26] - x[10];
169	0	r1 = x[27] - x[11];
170	0	x[26] += x[10];
171	0	x[27] += x[11];
172	0	x[10] = MULT_NORM(( r0 - r1 ) * cPI2_8);
173	0	x[11] = MULT_NORM(( r0 + r1 ) * cPI2_8);
174
175	0	r0 = x[24] - x[8];
176	0	r1 = x[25] - x[9];
177	0	x[24] += x[8];
178	0	x[25] += x[9];
179	0	x[8] = MULT_NORM( r0 * cPI3_8 - r1 * cPI1_8 );
180	0	x[9] = MULT_NORM( r1 * cPI3_8 + r0 * cPI1_8 );
181
182	0	r0 = x[22] - x[6];
183	0	r1 = x[7] - x[23];
184	0	x[22] += x[6];
185	0	x[23] += x[7];
186	0	x[6] = r1;
187	0	x[7] = r0;
188
189	0	r0 = x[4] - x[20];
190	0	r1 = x[5] - x[21];
191	0	x[20] += x[4];
192	0	x[21] += x[5];
193	0	x[4] = MULT_NORM( r1 * cPI1_8 + r0 * cPI3_8 );
194	0	x[5] = MULT_NORM( r1 * cPI3_8 - r0 * cPI1_8 );
195
196	0	r0 = x[2] - x[18];
197	0	r1 = x[3] - x[19];
198	0	x[18] += x[2];
199	0	x[19] += x[3];
200	0	x[2] = MULT_NORM(( r1 + r0 ) * cPI2_8);
201	0	x[3] = MULT_NORM(( r1 - r0 ) * cPI2_8);
202
203	0	r0 = x[0] - x[16];
204	0	r1 = x[1] - x[17];
205	0	x[16] += x[0];
206	0	x[17] += x[1];
207	0	x[0] = MULT_NORM( r1 * cPI3_8 + r0 * cPI1_8 );
208	0	x[1] = MULT_NORM( r1 * cPI1_8 - r0 * cPI3_8 );
209
210	0	mdct_butterfly_16(x);
211	0	mdct_butterfly_16(x+16);
212
213	0	}
214
215		/* N point first stage butterfly (in place, 2 register) */
216		STIN void mdct_butterfly_first(DATA_TYPE *T,
217		DATA_TYPE *x,
218	0	int points){
219
220	0	DATA_TYPE *x1 = x + points - 8;
221	0	DATA_TYPE *x2 = x + (points>>1) - 8;
222	0	REG_TYPE r0;
223	0	REG_TYPE r1;
224
225	0	do{
226
227	0	r0 = x1[6] - x2[6];
228	0	r1 = x1[7] - x2[7];
229	0	x1[6] += x2[6];
230	0	x1[7] += x2[7];
231	0	x2[6] = MULT_NORM(r1 * T[1] + r0 * T[0]);
232	0	x2[7] = MULT_NORM(r1 * T[0] - r0 * T[1]);
233
234	0	r0 = x1[4] - x2[4];
235	0	r1 = x1[5] - x2[5];
236	0	x1[4] += x2[4];
237	0	x1[5] += x2[5];
238	0	x2[4] = MULT_NORM(r1 * T[5] + r0 * T[4]);
239	0	x2[5] = MULT_NORM(r1 * T[4] - r0 * T[5]);
240
241	0	r0 = x1[2] - x2[2];
242	0	r1 = x1[3] - x2[3];
243	0	x1[2] += x2[2];
244	0	x1[3] += x2[3];
245	0	x2[2] = MULT_NORM(r1 * T[9] + r0 * T[8]);
246	0	x2[3] = MULT_NORM(r1 * T[8] - r0 * T[9]);
247
248	0	r0 = x1[0] - x2[0];
249	0	r1 = x1[1] - x2[1];
250	0	x1[0] += x2[0];
251	0	x1[1] += x2[1];
252	0	x2[0] = MULT_NORM(r1 * T[13] + r0 * T[12]);
253	0	x2[1] = MULT_NORM(r1 * T[12] - r0 * T[13]);
254
255	0	x1-=8;
256	0	x2-=8;
257	0	T+=16;
258
259	0	}while(x2>=x);
260	0	}
261
262		/* N/stage point generic N stage butterfly (in place, 2 register) */
263		STIN void mdct_butterfly_generic(DATA_TYPE *T,
264		DATA_TYPE *x,
265		int points,
266	0	int trigint){
267
268	0	DATA_TYPE *x1 = x + points - 8;
269	0	DATA_TYPE *x2 = x + (points>>1) - 8;
270	0	REG_TYPE r0;
271	0	REG_TYPE r1;
272
273	0	do{
274
275	0	r0 = x1[6] - x2[6];
276	0	r1 = x1[7] - x2[7];
277	0	x1[6] += x2[6];
278	0	x1[7] += x2[7];
279	0	x2[6] = MULT_NORM(r1 * T[1] + r0 * T[0]);
280	0	x2[7] = MULT_NORM(r1 * T[0] - r0 * T[1]);
281
282	0	T+=trigint;
283
284	0	r0 = x1[4] - x2[4];
285	0	r1 = x1[5] - x2[5];
286	0	x1[4] += x2[4];
287	0	x1[5] += x2[5];
288	0	x2[4] = MULT_NORM(r1 * T[1] + r0 * T[0]);
289	0	x2[5] = MULT_NORM(r1 * T[0] - r0 * T[1]);
290
291	0	T+=trigint;
292
293	0	r0 = x1[2] - x2[2];
294	0	r1 = x1[3] - x2[3];
295	0	x1[2] += x2[2];
296	0	x1[3] += x2[3];
297	0	x2[2] = MULT_NORM(r1 * T[1] + r0 * T[0]);
298	0	x2[3] = MULT_NORM(r1 * T[0] - r0 * T[1]);
299
300	0	T+=trigint;
301
302	0	r0 = x1[0] - x2[0];
303	0	r1 = x1[1] - x2[1];
304	0	x1[0] += x2[0];
305	0	x1[1] += x2[1];
306	0	x2[0] = MULT_NORM(r1 * T[1] + r0 * T[0]);
307	0	x2[1] = MULT_NORM(r1 * T[0] - r0 * T[1]);
308
309	0	T+=trigint;
310	0	x1-=8;
311	0	x2-=8;
312
313	0	}while(x2>=x);
314	0	}
315
316		STIN void mdct_butterflies(mdct_lookup *init,
317		DATA_TYPE *x,
318	0	int points){
319
320	0	DATA_TYPE *T=init->trig;
321	0	int stages=init->log2n-5;
322	0	int i,j;
323
324	0	if(--stages>0){
325	0	mdct_butterfly_first(T,x,points);
326	0	}
327
328	0	for(i=1;--stages>0;i++){
329	0	for(j=0;j<(1<<i);j++)
330	0	mdct_butterfly_generic(T,x+(points>>i)*j,points>>i,4<<i);
331	0	}
332
333	0	for(j=0;j<points;j+=32)
334	0	mdct_butterfly_32(x+j);
335
336	0	}
337
338	0	void mdct_clear(mdct_lookup *l){
339	0	if(l){
340	0	if(l->trig)_ogg_free(l->trig);
341	0	if(l->bitrev)_ogg_free(l->bitrev);
342	0	memset(l,0,sizeof(*l));
343	0	}
344	0	}
345
346		STIN void mdct_bitreverse(mdct_lookup *init,
347	0	DATA_TYPE *x){
348	0	int n = init->n;
349	0	int *bit = init->bitrev;
350	0	DATA_TYPE *w0 = x;
351	0	DATA_TYPE *w1 = x = w0+(n>>1);
352	0	DATA_TYPE *T = init->trig+n;
353
354	0	do{
355	0	DATA_TYPE *x0 = x+bit[0];
356	0	DATA_TYPE *x1 = x+bit[1];
357
358	0	REG_TYPE r0 = x0[1] - x1[1];
359	0	REG_TYPE r1 = x0[0] + x1[0];
360	0	REG_TYPE r2 = MULT_NORM(r1 * T[0] + r0 * T[1]);
361	0	REG_TYPE r3 = MULT_NORM(r1 * T[1] - r0 * T[0]);
362
363	0	w1 -= 4;
364
365	0	r0 = HALVE(x0[1] + x1[1]);
366	0	r1 = HALVE(x0[0] - x1[0]);
367
368	0	w0[0] = r0 + r2;
369	0	w1[2] = r0 - r2;
370	0	w0[1] = r1 + r3;
371	0	w1[3] = r3 - r1;
372
373	0	x0 = x+bit[2];
374	0	x1 = x+bit[3];
375
376	0	r0 = x0[1] - x1[1];
377	0	r1 = x0[0] + x1[0];
378	0	r2 = MULT_NORM(r1 * T[2] + r0 * T[3]);
379	0	r3 = MULT_NORM(r1 * T[3] - r0 * T[2]);
380
381	0	r0 = HALVE(x0[1] + x1[1]);
382	0	r1 = HALVE(x0[0] - x1[0]);
383
384	0	w0[2] = r0 + r2;
385	0	w1[0] = r0 - r2;
386	0	w0[3] = r1 + r3;
387	0	w1[1] = r3 - r1;
388
389	0	T += 4;
390	0	bit += 4;
391	0	w0 += 4;
392
393	0	}while(w0<w1);
394	0	}
395
396	0	void mdct_backward(mdct_lookup init, DATA_TYPE in, DATA_TYPE *out){
397	0	int n=init->n;
398	0	int n2=n>>1;
399	0	int n4=n>>2;
400
401		/* rotate */
402
403	0	DATA_TYPE *iX = in+n2-7;
404	0	DATA_TYPE *oX = out+n2+n4;
405	0	DATA_TYPE *T = init->trig+n4;
406
407	0	do{
408	0	oX -= 4;
409	0	oX[0] = MULT_NORM(-iX[2] * T[3] - iX[0] * T[2]);
410	0	oX[1] = MULT_NORM (iX[0] * T[3] - iX[2] * T[2]);
411	0	oX[2] = MULT_NORM(-iX[6] * T[1] - iX[4] * T[0]);
412	0	oX[3] = MULT_NORM (iX[4] * T[1] - iX[6] * T[0]);
413	0	iX -= 8;
414	0	T += 4;
415	0	}while(iX>=in);
416
417	0	iX = in+n2-8;
418	0	oX = out+n2+n4;
419	0	T = init->trig+n4;
420
421	0	do{
422	0	T -= 4;
423	0	oX[0] = MULT_NORM (iX[4] * T[3] + iX[6] * T[2]);
424	0	oX[1] = MULT_NORM (iX[4] * T[2] - iX[6] * T[3]);
425	0	oX[2] = MULT_NORM (iX[0] * T[1] + iX[2] * T[0]);
426	0	oX[3] = MULT_NORM (iX[0] * T[0] - iX[2] * T[1]);
427	0	iX -= 8;
428	0	oX += 4;
429	0	}while(iX>=in);
430
431	0	mdct_butterflies(init,out+n2,n2);
432	0	mdct_bitreverse(init,out);
433
434		/* roatate + window */
435
436	0	{
437	0	DATA_TYPE *oX1=out+n2+n4;
438	0	DATA_TYPE *oX2=out+n2+n4;
439	0	DATA_TYPE *iX =out;
440	0	T =init->trig+n2;
441
442	0	do{
443	0	oX1-=4;
444
445	0	oX1[3] = MULT_NORM (iX[0] * T[1] - iX[1] * T[0]);
446	0	oX2[0] = -MULT_NORM (iX[0] * T[0] + iX[1] * T[1]);
447
448	0	oX1[2] = MULT_NORM (iX[2] * T[3] - iX[3] * T[2]);
449	0	oX2[1] = -MULT_NORM (iX[2] * T[2] + iX[3] * T[3]);
450
451	0	oX1[1] = MULT_NORM (iX[4] * T[5] - iX[5] * T[4]);
452	0	oX2[2] = -MULT_NORM (iX[4] * T[4] + iX[5] * T[5]);
453
454	0	oX1[0] = MULT_NORM (iX[6] * T[7] - iX[7] * T[6]);
455	0	oX2[3] = -MULT_NORM (iX[6] * T[6] + iX[7] * T[7]);
456
457	0	oX2+=4;
458	0	iX += 8;
459	0	T += 8;
460	0	}while(iX<oX1);
461
462	0	iX=out+n2+n4;
463	0	oX1=out+n4;
464	0	oX2=oX1;
465
466	0	do{
467	0	oX1-=4;
468	0	iX-=4;
469
470	0	oX2[0] = -(oX1[3] = iX[3]);
471	0	oX2[1] = -(oX1[2] = iX[2]);
472	0	oX2[2] = -(oX1[1] = iX[1]);
473	0	oX2[3] = -(oX1[0] = iX[0]);
474
475	0	oX2+=4;
476	0	}while(oX2<iX);
477
478	0	iX=out+n2+n4;
479	0	oX1=out+n2+n4;
480	0	oX2=out+n2;
481	0	do{
482	0	oX1-=4;
483	0	oX1[0]= iX[3];
484	0	oX1[1]= iX[2];
485	0	oX1[2]= iX[1];
486	0	oX1[3]= iX[0];
487	0	iX+=4;
488	0	}while(oX1>oX2);
489	0	}
490	0	}
491
492	0	void mdct_forward(mdct_lookup init, DATA_TYPE in, DATA_TYPE *out){
493	0	int n=init->n;
494	0	int n2=n>>1;
495	0	int n4=n>>2;
496	0	int n8=n>>3;
497	0	DATA_TYPE w=alloca(nsizeof(w)); / forward needs working space */
498	0	DATA_TYPE *w2=w+n2;
499
500		/* rotate */
501
502		/* window + rotate + step 1 */
503
504	0	REG_TYPE r0;
505	0	REG_TYPE r1;
506	0	DATA_TYPE *x0=in+n2+n4;
507	0	DATA_TYPE *x1=x0+1;
508	0	DATA_TYPE *T=init->trig+n2;
509
510	0	int i=0;
511
512	0	for(i=0;i<n8;i+=2){
513	0	x0 -=4;
514	0	T-=2;
515	0	r0= x0[2] + x1[0];
516	0	r1= x0[0] + x1[2];
517	0	w2[i]= MULT_NORM(r1T[1] + r0T[0]);
518	0	w2[i+1]= MULT_NORM(r1T[0] - r0T[1]);
519	0	x1 +=4;
520	0	}
521
522	0	x1=in+1;
523
524	0	for(;i<n2-n8;i+=2){
525	0	T-=2;
526	0	x0 -=4;
527	0	r0= x0[2] - x1[0];
528	0	r1= x0[0] - x1[2];
529	0	w2[i]= MULT_NORM(r1T[1] + r0T[0]);
530	0	w2[i+1]= MULT_NORM(r1T[0] - r0T[1]);
531	0	x1 +=4;
532	0	}
533
534	0	x0=in+n;
535
536	0	for(;i<n2;i+=2){
537	0	T-=2;
538	0	x0 -=4;
539	0	r0= -x0[2] - x1[0];
540	0	r1= -x0[0] - x1[2];
541	0	w2[i]= MULT_NORM(r1T[1] + r0T[0]);
542	0	w2[i+1]= MULT_NORM(r1T[0] - r0T[1]);
543	0	x1 +=4;
544	0	}
545
546
547	0	mdct_butterflies(init,w+n2,n2);
548	0	mdct_bitreverse(init,w);
549
550		/* roatate + window */
551
552	0	T=init->trig+n2;
553	0	x0=out+n2;
554
555	0	for(i=0;i<n4;i++){
556	0	x0--;
557	0	out[i] =MULT_NORM((w[0]T[0]+w[1]T[1])*init->scale);
558	0	x0[0] =MULT_NORM((w[0]T[1]-w[1]T[0])*init->scale);
559	0	w+=2;
560	0	T+=2;
561	0	}
562	0	}

Coverage Report

Created: 2024-09-06 07:53