ffmpeg / libavcodec / rdft.c @ 2881c831
History  View  Annotate  Download (4.15 KB)
1 
/*


2 
* (I)RDFT transforms

3 
* Copyright (c) 2009 Alex Converse <alex dot converse at gmail dot com>

4 
*

5 
* This file is part of FFmpeg.

6 
*

7 
* FFmpeg is free software; you can redistribute it and/or

8 
* modify it under the terms of the GNU Lesser General Public

9 
* License as published by the Free Software Foundation; either

10 
* version 2.1 of the License, or (at your option) any later version.

11 
*

12 
* FFmpeg is distributed in the hope that it will be useful,

13 
* but WITHOUT ANY WARRANTY; without even the implied warranty of

14 
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

15 
* Lesser General Public License for more details.

16 
*

17 
* You should have received a copy of the GNU Lesser General Public

18 
* License along with FFmpeg; if not, write to the Free Software

19 
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA

20 
*/

21 
#include <stdlib.h> 
22 
#include <math.h> 
23 
#include "libavutil/mathematics.h" 
24 
#include "fft.h" 
25  
26 
/**

27 
* @file libavcodec/rdft.c

28 
* (Inverse) Real Discrete Fourier Transforms.

29 
*/

30  
31 
/* sin(2*pi*x/n) for 0<=x<n/4, followed by n/2<=x<3n/4 */

32 
#if !CONFIG_HARDCODED_TABLES

33 
SINTABLE(16);

34 
SINTABLE(32);

35 
SINTABLE(64);

36 
SINTABLE(128);

37 
SINTABLE(256);

38 
SINTABLE(512);

39 
SINTABLE(1024);

40 
SINTABLE(2048);

41 
SINTABLE(4096);

42 
SINTABLE(8192);

43 
SINTABLE(16384);

44 
SINTABLE(32768);

45 
SINTABLE(65536);

46 
#endif

47 
SINTABLE_CONST FFTSample * const ff_sin_tabs[] = {

48 
NULL, NULL, NULL, NULL, 
49 
ff_sin_16, ff_sin_32, ff_sin_64, ff_sin_128, ff_sin_256, ff_sin_512, ff_sin_1024, 
50 
ff_sin_2048, ff_sin_4096, ff_sin_8192, ff_sin_16384, ff_sin_32768, ff_sin_65536, 
51 
}; 
52  
53 
static void ff_rdft_calc_c(RDFTContext* s, FFTSample* data); 
54  
55 
av_cold int ff_rdft_init(RDFTContext *s, int nbits, enum RDFTransformType trans) 
56 
{ 
57 
int n = 1 << nbits; 
58 
int i;

59 
const double theta = (trans == DFT_R2C  trans == DFT_C2R ? 1 : 1)*2*M_PI/n; 
60  
61 
s>nbits = nbits; 
62 
s>inverse = trans == IDFT_C2R  trans == DFT_C2R; 
63 
s>sign_convention = trans == IDFT_R2C  trans == DFT_C2R ? 1 : 1; 
64  
65 
if (nbits < 4  nbits > 16) 
66 
return 1; 
67  
68 
if (ff_fft_init(&s>fft, nbits1, trans == IDFT_C2R  trans == IDFT_R2C) < 0) 
69 
return 1; 
70  
71 
ff_init_ff_cos_tabs(nbits); 
72 
s>tcos = ff_cos_tabs[nbits]; 
73 
s>tsin = ff_sin_tabs[nbits]+(trans == DFT_R2C  trans == DFT_C2R)*(n>>2);

74 
#if !CONFIG_HARDCODED_TABLES

75 
for (i = 0; i < (n>>2); i++) { 
76 
s>tsin[i] = sin(i*theta); 
77 
} 
78 
#endif

79 
s>rdft_calc = ff_rdft_calc_c; 
80 
return 0; 
81 
} 
82  
83 
/** Map one real FFT into two parallel real even and odd FFTs. Then interleave

84 
* the two real FFTs into one complex FFT. Unmangle the results.

85 
* ref: http://www.engineeringproductivitytools.com/stuff/T0001/PT10.HTM

86 
*/

87 
static void ff_rdft_calc_c(RDFTContext* s, FFTSample* data) 
88 
{ 
89 
int i, i1, i2;

90 
FFTComplex ev, od; 
91 
const int n = 1 << s>nbits; 
92 
const float k1 = 0.5; 
93 
const float k2 = 0.5  s>inverse; 
94 
const FFTSample *tcos = s>tcos;

95 
const FFTSample *tsin = s>tsin;

96  
97 
if (!s>inverse) {

98 
ff_fft_permute(&s>fft, (FFTComplex*)data); 
99 
ff_fft_calc(&s>fft, (FFTComplex*)data); 
100 
} 
101 
/* i=0 is a special case because of packing, the DC term is real, so we

102 
are going to throw the N/2 term (also real) in with it. */

103 
ev.re = data[0];

104 
data[0] = ev.re+data[1]; 
105 
data[1] = ev.redata[1]; 
106 
for (i = 1; i < (n>>2); i++) { 
107 
i1 = 2*i;

108 
i2 = ni1; 
109 
/* Separate even and odd FFTs */

110 
ev.re = k1*(data[i1 ]+data[i2 ]); 
111 
od.im = k2*(data[i1 ]data[i2 ]); 
112 
ev.im = k1*(data[i1+1]data[i2+1]); 
113 
od.re = k2*(data[i1+1]+data[i2+1]); 
114 
/* Apply twiddle factors to the odd FFT and add to the even FFT */

115 
data[i1 ] = ev.re + od.re*tcos[i]  od.im*tsin[i]; 
116 
data[i1+1] = ev.im + od.im*tcos[i] + od.re*tsin[i];

117 
data[i2 ] = ev.re  od.re*tcos[i] + od.im*tsin[i]; 
118 
data[i2+1] = ev.im + od.im*tcos[i] + od.re*tsin[i];

119 
} 
120 
data[2*i+1]=s>sign_convention*data[2*i+1]; 
121 
if (s>inverse) {

122 
data[0] *= k1;

123 
data[1] *= k1;

124 
ff_fft_permute(&s>fft, (FFTComplex*)data); 
125 
ff_fft_calc(&s>fft, (FFTComplex*)data); 
126 
} 
127 
} 
128  
129 
av_cold void ff_rdft_end(RDFTContext *s)

130 
{ 
131 
ff_fft_end(&s>fft); 
132 
} 