/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	18.3k	void mdct_init(mdct_lookup *lookup,int n){
52	18.3k	int bitrev=_ogg_malloc(sizeof(bitrev)*(n/4));
53	18.3k	DATA_TYPE T=_ogg_malloc(sizeof(T)*(n+n/4));
54
55	18.3k	int i;
56	18.3k	int n2=n>>1;
57	18.3k	int log2n=lookup->log2n=rint(log((float)n)/log(2.f));
58	18.3k	lookup->n=n;
59	18.3k	lookup->trig=T;
60	18.3k	lookup->bitrev=bitrev;
61
62		/* trig lookups... */
63
64	5.69M	for(i=0;i<n/4;i++){
65	5.67M	T[i2]=FLOAT_CONV(cos((M_PI/n)(4*i)));
66	5.67M	T[i2+1]=FLOAT_CONV(-sin((M_PI/n)(4*i)));
67	5.67M	T[n2+i2]=FLOAT_CONV(cos((M_PI/(2n))(2i+1)));
68	5.67M	T[n2+i2+1]=FLOAT_CONV(sin((M_PI/(2n))(2i+1)));
69	5.67M	}
70	2.85M	for(i=0;i<n/8;i++){
71	2.83M	T[n+i2]=FLOAT_CONV(cos((M_PI/n)(4i+2)).5);
72	2.83M	T[n+i2+1]=FLOAT_CONV(-sin((M_PI/n)(4i+2)).5);
73	2.83M	}
74
75		/* bitreverse lookup... */
76
77	18.3k	{
78	18.3k	int mask=(1<<(log2n-1))-1,i,j;
79	18.3k	int msb=1<<(log2n-2);
80	2.85M	for(i=0;i<n/8;i++){
81	2.83M	int acc=0;
82	32.4M	for(j=0;msb>>j;j++)
83	29.5M	if((msb>>j)&i)acc\|=1<<j;
84	2.83M	bitrev[i*2]=((~acc)&mask)-1;
85	2.83M	bitrev[i*2+1]=acc;
86
87	2.83M	}
88	18.3k	}
89	18.3k	lookup->scale=FLOAT_CONV(4.f/n);
90	18.3k	}
91
92		/* 8 point butterfly (in place, 4 register) */
93	261M	STIN void mdct_butterfly_8(DATA_TYPE *x){
94	261M	REG_TYPE r0 = x[6] + x[2];
95	261M	REG_TYPE r1 = x[6] - x[2];
96	261M	REG_TYPE r2 = x[4] + x[0];
97	261M	REG_TYPE r3 = x[4] - x[0];
98
99	261M	x[6] = r0 + r2;
100	261M	x[4] = r0 - r2;
101
102	261M	r0 = x[5] - x[1];
103	261M	r2 = x[7] - x[3];
104	261M	x[0] = r1 + r0;
105	261M	x[2] = r1 - r0;
106
107	261M	r0 = x[5] + x[1];
108	261M	r1 = x[7] + x[3];
109	261M	x[3] = r2 + r3;
110	261M	x[1] = r2 - r3;
111	261M	x[7] = r1 + r0;
112	261M	x[5] = r1 - r0;
113
114	261M	}
115
116		/* 16 point butterfly (in place, 4 register) */
117	130M	STIN void mdct_butterfly_16(DATA_TYPE *x){
118	130M	REG_TYPE r0 = x[1] - x[9];
119	130M	REG_TYPE r1 = x[0] - x[8];
120
121	130M	x[8] += x[0];
122	130M	x[9] += x[1];
123	130M	x[0] = MULT_NORM((r0 + r1) * cPI2_8);
124	130M	x[1] = MULT_NORM((r0 - r1) * cPI2_8);
125
126	130M	r0 = x[3] - x[11];
127	130M	r1 = x[10] - x[2];
128	130M	x[10] += x[2];
129	130M	x[11] += x[3];
130	130M	x[2] = r0;
131	130M	x[3] = r1;
132
133	130M	r0 = x[12] - x[4];
134	130M	r1 = x[13] - x[5];
135	130M	x[12] += x[4];
136	130M	x[13] += x[5];
137	130M	x[4] = MULT_NORM((r0 - r1) * cPI2_8);
138	130M	x[5] = MULT_NORM((r0 + r1) * cPI2_8);
139
140	130M	r0 = x[14] - x[6];
141	130M	r1 = x[15] - x[7];
142	130M	x[14] += x[6];
143	130M	x[15] += x[7];
144	130M	x[6] = r0;
145	130M	x[7] = r1;
146
147	130M	mdct_butterfly_8(x);
148	130M	mdct_butterfly_8(x+8);
149	130M	}
150
151		/* 32 point butterfly (in place, 4 register) */
152	65.2M	STIN void mdct_butterfly_32(DATA_TYPE *x){
153	65.2M	REG_TYPE r0 = x[30] - x[14];
154	65.2M	REG_TYPE r1 = x[31] - x[15];
155
156	65.2M	x[30] += x[14];
157	65.2M	x[31] += x[15];
158	65.2M	x[14] = r0;
159	65.2M	x[15] = r1;
160
161	65.2M	r0 = x[28] - x[12];
162	65.2M	r1 = x[29] - x[13];
163	65.2M	x[28] += x[12];
164	65.2M	x[29] += x[13];
165	65.2M	x[12] = MULT_NORM( r0 * cPI1_8 - r1 * cPI3_8 );
166	65.2M	x[13] = MULT_NORM( r0 * cPI3_8 + r1 * cPI1_8 );
167
168	65.2M	r0 = x[26] - x[10];
169	65.2M	r1 = x[27] - x[11];
170	65.2M	x[26] += x[10];
171	65.2M	x[27] += x[11];
172	65.2M	x[10] = MULT_NORM(( r0 - r1 ) * cPI2_8);
173	65.2M	x[11] = MULT_NORM(( r0 + r1 ) * cPI2_8);
174
175	65.2M	r0 = x[24] - x[8];
176	65.2M	r1 = x[25] - x[9];
177	65.2M	x[24] += x[8];
178	65.2M	x[25] += x[9];
179	65.2M	x[8] = MULT_NORM( r0 * cPI3_8 - r1 * cPI1_8 );
180	65.2M	x[9] = MULT_NORM( r1 * cPI3_8 + r0 * cPI1_8 );
181
182	65.2M	r0 = x[22] - x[6];
183	65.2M	r1 = x[7] - x[23];
184	65.2M	x[22] += x[6];
185	65.2M	x[23] += x[7];
186	65.2M	x[6] = r1;
187	65.2M	x[7] = r0;
188
189	65.2M	r0 = x[4] - x[20];
190	65.2M	r1 = x[5] - x[21];
191	65.2M	x[20] += x[4];
192	65.2M	x[21] += x[5];
193	65.2M	x[4] = MULT_NORM( r1 * cPI1_8 + r0 * cPI3_8 );
194	65.2M	x[5] = MULT_NORM( r1 * cPI3_8 - r0 * cPI1_8 );
195
196	65.2M	r0 = x[2] - x[18];
197	65.2M	r1 = x[3] - x[19];
198	65.2M	x[18] += x[2];
199	65.2M	x[19] += x[3];
200	65.2M	x[2] = MULT_NORM(( r1 + r0 ) * cPI2_8);
201	65.2M	x[3] = MULT_NORM(( r1 - r0 ) * cPI2_8);
202
203	65.2M	r0 = x[0] - x[16];
204	65.2M	r1 = x[1] - x[17];
205	65.2M	x[16] += x[0];
206	65.2M	x[17] += x[1];
207	65.2M	x[0] = MULT_NORM( r1 * cPI3_8 + r0 * cPI1_8 );
208	65.2M	x[1] = MULT_NORM( r1 * cPI1_8 - r0 * cPI3_8 );
209
210	65.2M	mdct_butterfly_16(x);
211	65.2M	mdct_butterfly_16(x+16);
212
213	65.2M	}
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	2.52M	int points){
219
220	2.52M	DATA_TYPE *x1 = x + points - 8;
221	2.52M	DATA_TYPE *x2 = x + (points>>1) - 8;
222	2.52M	REG_TYPE r0;
223	2.52M	REG_TYPE r1;
224
225	130M	do{
226
227	130M	r0 = x1[6] - x2[6];
228	130M	r1 = x1[7] - x2[7];
229	130M	x1[6] += x2[6];
230	130M	x1[7] += x2[7];
231	130M	x2[6] = MULT_NORM(r1 * T[1] + r0 * T[0]);
232	130M	x2[7] = MULT_NORM(r1 * T[0] - r0 * T[1]);
233
234	130M	r0 = x1[4] - x2[4];
235	130M	r1 = x1[5] - x2[5];
236	130M	x1[4] += x2[4];
237	130M	x1[5] += x2[5];
238	130M	x2[4] = MULT_NORM(r1 * T[5] + r0 * T[4]);
239	130M	x2[5] = MULT_NORM(r1 * T[4] - r0 * T[5]);
240
241	130M	r0 = x1[2] - x2[2];
242	130M	r1 = x1[3] - x2[3];
243	130M	x1[2] += x2[2];
244	130M	x1[3] += x2[3];
245	130M	x2[2] = MULT_NORM(r1 * T[9] + r0 * T[8]);
246	130M	x2[3] = MULT_NORM(r1 * T[8] - r0 * T[9]);
247
248	130M	r0 = x1[0] - x2[0];
249	130M	r1 = x1[1] - x2[1];
250	130M	x1[0] += x2[0];
251	130M	x1[1] += x2[1];
252	130M	x2[0] = MULT_NORM(r1 * T[13] + r0 * T[12]);
253	130M	x2[1] = MULT_NORM(r1 * T[12] - r0 * T[13]);
254
255	130M	x1-=8;
256	130M	x2-=8;
257	130M	T+=16;
258
259	130M	}while(x2>=x);
260	2.52M	}
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	60.0M	int trigint){
267
268	60.0M	DATA_TYPE *x1 = x + points - 8;
269	60.0M	DATA_TYPE *x2 = x + (points>>1) - 8;
270	60.0M	REG_TYPE r0;
271	60.0M	REG_TYPE r1;
272
273	576M	do{
274
275	576M	r0 = x1[6] - x2[6];
276	576M	r1 = x1[7] - x2[7];
277	576M	x1[6] += x2[6];
278	576M	x1[7] += x2[7];
279	576M	x2[6] = MULT_NORM(r1 * T[1] + r0 * T[0]);
280	576M	x2[7] = MULT_NORM(r1 * T[0] - r0 * T[1]);
281
282	576M	T+=trigint;
283
284	576M	r0 = x1[4] - x2[4];
285	576M	r1 = x1[5] - x2[5];
286	576M	x1[4] += x2[4];
287	576M	x1[5] += x2[5];
288	576M	x2[4] = MULT_NORM(r1 * T[1] + r0 * T[0]);
289	576M	x2[5] = MULT_NORM(r1 * T[0] - r0 * T[1]);
290
291	576M	T+=trigint;
292
293	576M	r0 = x1[2] - x2[2];
294	576M	r1 = x1[3] - x2[3];
295	576M	x1[2] += x2[2];
296	576M	x1[3] += x2[3];
297	576M	x2[2] = MULT_NORM(r1 * T[1] + r0 * T[0]);
298	576M	x2[3] = MULT_NORM(r1 * T[0] - r0 * T[1]);
299
300	576M	T+=trigint;
301
302	576M	r0 = x1[0] - x2[0];
303	576M	r1 = x1[1] - x2[1];
304	576M	x1[0] += x2[0];
305	576M	x1[1] += x2[1];
306	576M	x2[0] = MULT_NORM(r1 * T[1] + r0 * T[0]);
307	576M	x2[1] = MULT_NORM(r1 * T[0] - r0 * T[1]);
308
309	576M	T+=trigint;
310	576M	x1-=8;
311	576M	x2-=8;
312
313	576M	}while(x2>=x);
314	60.0M	}
315
316		STIN void mdct_butterflies(mdct_lookup *init,
317		DATA_TYPE *x,
318	2.69M	int points){
319
320	2.69M	DATA_TYPE *T=init->trig;
321	2.69M	int stages=init->log2n-5;
322	2.69M	int i,j;
323
324	2.69M	if(--stages>0){
325	2.52M	mdct_butterfly_first(T,x,points);
326	2.52M	}
327
328	9.99M	for(i=1;--stages>0;i++){
329	67.3M	for(j=0;j<(1<<i);j++)
330	60.0M	mdct_butterfly_generic(T,x+(points>>i)*j,points>>i,4<<i);
331	7.30M	}
332
333	67.9M	for(j=0;j<points;j+=32)
334	65.2M	mdct_butterfly_32(x+j);
335
336	2.69M	}
337
338	18.3k	void mdct_clear(mdct_lookup *l){
339	18.3k	if(l){
340	18.3k	if(l->trig)_ogg_free(l->trig);
341	18.3k	if(l->bitrev)_ogg_free(l->bitrev);
342	18.3k	memset(l,0,sizeof(*l));
343	18.3k	}
344	18.3k	}
345
346		STIN void mdct_bitreverse(mdct_lookup *init,
347	2.69M	DATA_TYPE *x){
348	2.69M	int n = init->n;
349	2.69M	int *bit = init->bitrev;
350	2.69M	DATA_TYPE *w0 = x;
351	2.69M	DATA_TYPE *w1 = x = w0+(n>>1);
352	2.69M	DATA_TYPE *T = init->trig+n;
353
354	261M	do{
355	261M	DATA_TYPE *x0 = x+bit[0];
356	261M	DATA_TYPE *x1 = x+bit[1];
357
358	261M	REG_TYPE r0 = x0[1] - x1[1];
359	261M	REG_TYPE r1 = x0[0] + x1[0];
360	261M	REG_TYPE r2 = MULT_NORM(r1 * T[0] + r0 * T[1]);
361	261M	REG_TYPE r3 = MULT_NORM(r1 * T[1] - r0 * T[0]);
362
363	261M	w1 -= 4;
364
365	261M	r0 = HALVE(x0[1] + x1[1]);
366	261M	r1 = HALVE(x0[0] - x1[0]);
367
368	261M	w0[0] = r0 + r2;
369	261M	w1[2] = r0 - r2;
370	261M	w0[1] = r1 + r3;
371	261M	w1[3] = r3 - r1;
372
373	261M	x0 = x+bit[2];
374	261M	x1 = x+bit[3];
375
376	261M	r0 = x0[1] - x1[1];
377	261M	r1 = x0[0] + x1[0];
378	261M	r2 = MULT_NORM(r1 * T[2] + r0 * T[3]);
379	261M	r3 = MULT_NORM(r1 * T[3] - r0 * T[2]);
380
381	261M	r0 = HALVE(x0[1] + x1[1]);
382	261M	r1 = HALVE(x0[0] - x1[0]);
383
384	261M	w0[2] = r0 + r2;
385	261M	w1[0] = r0 - r2;
386	261M	w0[3] = r1 + r3;
387	261M	w1[1] = r3 - r1;
388
389	261M	T += 4;
390	261M	bit += 4;
391	261M	w0 += 4;
392
393	261M	}while(w0<w1);
394	2.69M	}
395
396	2.69M	void mdct_backward(mdct_lookup init, DATA_TYPE in, DATA_TYPE *out){
397	2.69M	int n=init->n;
398	2.69M	int n2=n>>1;
399	2.69M	int n4=n>>2;
400
401		/* rotate */
402
403	2.69M	DATA_TYPE *iX = in+n2-7;
404	2.69M	DATA_TYPE *oX = out+n2+n4;
405	2.69M	DATA_TYPE *T = init->trig+n4;
406
407	261M	do{
408	261M	oX -= 4;
409	261M	oX[0] = MULT_NORM(-iX[2] * T[3] - iX[0] * T[2]);
410	261M	oX[1] = MULT_NORM (iX[0] * T[3] - iX[2] * T[2]);
411	261M	oX[2] = MULT_NORM(-iX[6] * T[1] - iX[4] * T[0]);
412	261M	oX[3] = MULT_NORM (iX[4] * T[1] - iX[6] * T[0]);
413	261M	iX -= 8;
414	261M	T += 4;
415	261M	}while(iX>=in);
416
417	2.69M	iX = in+n2-8;
418	2.69M	oX = out+n2+n4;
419	2.69M	T = init->trig+n4;
420
421	261M	do{
422	261M	T -= 4;
423	261M	oX[0] = MULT_NORM (iX[4] * T[3] + iX[6] * T[2]);
424	261M	oX[1] = MULT_NORM (iX[4] * T[2] - iX[6] * T[3]);
425	261M	oX[2] = MULT_NORM (iX[0] * T[1] + iX[2] * T[0]);
426	261M	oX[3] = MULT_NORM (iX[0] * T[0] - iX[2] * T[1]);
427	261M	iX -= 8;
428	261M	oX += 4;
429	261M	}while(iX>=in);
430
431	2.69M	mdct_butterflies(init,out+n2,n2);
432	2.69M	mdct_bitreverse(init,out);
433
434		/* roatate + window */
435
436	2.69M	{
437	2.69M	DATA_TYPE *oX1=out+n2+n4;
438	2.69M	DATA_TYPE *oX2=out+n2+n4;
439	2.69M	DATA_TYPE *iX =out;
440	2.69M	T =init->trig+n2;
441
442	261M	do{
443	261M	oX1-=4;
444
445	261M	oX1[3] = MULT_NORM (iX[0] * T[1] - iX[1] * T[0]);
446	261M	oX2[0] = -MULT_NORM (iX[0] * T[0] + iX[1] * T[1]);
447
448	261M	oX1[2] = MULT_NORM (iX[2] * T[3] - iX[3] * T[2]);
449	261M	oX2[1] = -MULT_NORM (iX[2] * T[2] + iX[3] * T[3]);
450
451	261M	oX1[1] = MULT_NORM (iX[4] * T[5] - iX[5] * T[4]);
452	261M	oX2[2] = -MULT_NORM (iX[4] * T[4] + iX[5] * T[5]);
453
454	261M	oX1[0] = MULT_NORM (iX[6] * T[7] - iX[7] * T[6]);
455	261M	oX2[3] = -MULT_NORM (iX[6] * T[6] + iX[7] * T[7]);
456
457	261M	oX2+=4;
458	261M	iX += 8;
459	261M	T += 8;
460	261M	}while(iX<oX1);
461
462	2.69M	iX=out+n2+n4;
463	2.69M	oX1=out+n4;
464	2.69M	oX2=oX1;
465
466	261M	do{
467	261M	oX1-=4;
468	261M	iX-=4;
469
470	261M	oX2[0] = -(oX1[3] = iX[3]);
471	261M	oX2[1] = -(oX1[2] = iX[2]);
472	261M	oX2[2] = -(oX1[1] = iX[1]);
473	261M	oX2[3] = -(oX1[0] = iX[0]);
474
475	261M	oX2+=4;
476	261M	}while(oX2<iX);
477
478	2.69M	iX=out+n2+n4;
479	2.69M	oX1=out+n2+n4;
480	2.69M	oX2=out+n2;
481	261M	do{
482	261M	oX1-=4;
483	261M	oX1[0]= iX[3];
484	261M	oX1[1]= iX[2];
485	261M	oX1[2]= iX[1];
486	261M	oX1[3]= iX[0];
487	261M	iX+=4;
488	261M	}while(oX1>oX2);
489	2.69M	}
490	2.69M	}
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: 2025-07-12 06:11