Planet

navi

home

PPS

about

screenshots

download

development

forum

Context Navigation

source: downloads/libvorbis-1.2.0/lib/mdct.c @ 35

Last change on this file since 35 was 16, checked in by landauf, 18 years ago
added libvorbis
File size: 14.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: