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