ffmpeg / libavcodec / rdft.c @ 2912e87a
History  View  Annotate  Download (4.12 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 Libav.

6 
*

7 
* Libav 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 
* Libav 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 Libav; 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

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 
static 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 
/** Map one real FFT into two parallel real even and odd FFTs. Then interleave

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

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

56 
*/

57 
static void ff_rdft_calc_c(RDFTContext* s, FFTSample* data) 
58 
{ 
59 
int i, i1, i2;

60 
FFTComplex ev, od; 
61 
const int n = 1 << s>nbits; 
62 
const float k1 = 0.5; 
63 
const float k2 = 0.5  s>inverse; 
64 
const FFTSample *tcos = s>tcos;

65 
const FFTSample *tsin = s>tsin;

66  
67 
if (!s>inverse) {

68 
ff_fft_permute(&s>fft, (FFTComplex*)data); 
69 
ff_fft_calc(&s>fft, (FFTComplex*)data); 
70 
} 
71 
/* i=0 is a special case because of packing, the DC term is real, so we

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

73 
ev.re = data[0];

74 
data[0] = ev.re+data[1]; 
75 
data[1] = ev.redata[1]; 
76 
for (i = 1; i < (n>>2); i++) { 
77 
i1 = 2*i;

78 
i2 = ni1; 
79 
/* Separate even and odd FFTs */

80 
ev.re = k1*(data[i1 ]+data[i2 ]); 
81 
od.im = k2*(data[i1 ]data[i2 ]); 
82 
ev.im = k1*(data[i1+1]data[i2+1]); 
83 
od.re = k2*(data[i1+1]+data[i2+1]); 
84 
/* Apply twiddle factors to the odd FFT and add to the even FFT */

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

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

89 
} 
90 
data[2*i+1]=s>sign_convention*data[2*i+1]; 
91 
if (s>inverse) {

92 
data[0] *= k1;

93 
data[1] *= k1;

94 
ff_fft_permute(&s>fft, (FFTComplex*)data); 
95 
ff_fft_calc(&s>fft, (FFTComplex*)data); 
96 
} 
97 
} 
98  
99 
av_cold int ff_rdft_init(RDFTContext *s, int nbits, enum RDFTransformType trans) 
100 
{ 
101 
int n = 1 << nbits; 
102 
int i;

103 
const double theta = (trans == DFT_R2C  trans == DFT_C2R ? 1 : 1)*2*M_PI/n; 
104  
105 
s>nbits = nbits; 
106 
s>inverse = trans == IDFT_C2R  trans == DFT_C2R; 
107 
s>sign_convention = trans == IDFT_R2C  trans == DFT_C2R ? 1 : 1; 
108  
109 
if (nbits < 4  nbits > 16) 
110 
return 1; 
111  
112 
if (ff_fft_init(&s>fft, nbits1, trans == IDFT_C2R  trans == IDFT_R2C) < 0) 
113 
return 1; 
114  
115 
ff_init_ff_cos_tabs(nbits); 
116 
s>tcos = ff_cos_tabs[nbits]; 
117 
s>tsin = ff_sin_tabs[nbits]+(trans == DFT_R2C  trans == DFT_C2R)*(n>>2);

118 
#if !CONFIG_HARDCODED_TABLES

119 
for (i = 0; i < (n>>2); i++) { 
120 
s>tsin[i] = sin(i*theta); 
121 
} 
122 
#endif

123 
s>rdft_calc = ff_rdft_calc_c; 
124  
125 
if (ARCH_ARM) ff_rdft_init_arm(s);

126  
127 
return 0; 
128 
} 
129  
130 
av_cold void ff_rdft_end(RDFTContext *s)

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