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