/src/vorbis/lib/mdct.c

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