FFmpeg
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
rdft.c
Go to the documentation of this file.
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 02110-1301 USA
20  */
21 #include <stdlib.h>
22 #include <math.h>
23 #include "libavutil/mathematics.h"
24 #include "rdft.h"
25 
26 /**
27  * @file
28  * (Inverse) Real Discrete Fourier Transforms.
29  */
30 
31 /** Map one real FFT into two parallel real even and odd FFTs. Then interleave
32  * the two real FFTs into one complex FFT. Unmangle the results.
33  * ref: http://www.engineeringproductivitytools.com/stuff/T0001/PT10.HTM
34  */
36 {
37  int i, i1, i2;
38  FFTComplex ev, od;
39  const int n = 1 << s->nbits;
40  const float k1 = 0.5;
41  const float k2 = 0.5 - s->inverse;
42  const FFTSample *tcos = s->tcos;
43  const FFTSample *tsin = s->tsin;
44 
45  if (!s->inverse) {
46  s->fft.fft_permute(&s->fft, (FFTComplex*)data);
47  s->fft.fft_calc(&s->fft, (FFTComplex*)data);
48  }
49  /* i=0 is a special case because of packing, the DC term is real, so we
50  are going to throw the N/2 term (also real) in with it. */
51  ev.re = data[0];
52  data[0] = ev.re+data[1];
53  data[1] = ev.re-data[1];
54 
55 #define RDFT_UNMANGLE(sign0, sign1) \
56  for (i = 1; i < (n>>2); i++) { \
57  i1 = 2*i; \
58  i2 = n-i1; \
59  /* Separate even and odd FFTs */ \
60  ev.re = k1*(data[i1 ]+data[i2 ]); \
61  od.im = -k2*(data[i1 ]-data[i2 ]); \
62  ev.im = k1*(data[i1+1]-data[i2+1]); \
63  od.re = k2*(data[i1+1]+data[i2+1]); \
64  /* Apply twiddle factors to the odd FFT and add to the even FFT */ \
65  data[i1 ] = ev.re + od.re*tcos[i] sign0 od.im*tsin[i]; \
66  data[i1+1] = ev.im + od.im*tcos[i] sign1 od.re*tsin[i]; \
67  data[i2 ] = ev.re - od.re*tcos[i] sign1 od.im*tsin[i]; \
68  data[i2+1] = -ev.im + od.im*tcos[i] sign1 od.re*tsin[i]; \
69  }
70 
71  if (s->negative_sin) {
72  RDFT_UNMANGLE(+,-)
73  } else {
74  RDFT_UNMANGLE(-,+)
75  }
76 
77  data[2*i+1]=s->sign_convention*data[2*i+1];
78  if (s->inverse) {
79  data[0] *= k1;
80  data[1] *= k1;
81  s->fft.fft_permute(&s->fft, (FFTComplex*)data);
82  s->fft.fft_calc(&s->fft, (FFTComplex*)data);
83  }
84 }
85 
86 av_cold int ff_rdft_init(RDFTContext *s, int nbits, enum RDFTransformType trans)
87 {
88  int n = 1 << nbits;
89  int ret;
90 
91  s->nbits = nbits;
92  s->inverse = trans == IDFT_C2R || trans == DFT_C2R;
93  s->sign_convention = trans == IDFT_R2C || trans == DFT_C2R ? 1 : -1;
94  s->negative_sin = trans == DFT_C2R || trans == DFT_R2C;
95 
96  if (nbits < 4 || nbits > 16)
97  return AVERROR(EINVAL);
98 
99  if ((ret = ff_fft_init(&s->fft, nbits-1, trans == IDFT_C2R || trans == IDFT_R2C)) < 0)
100  return ret;
101 
102  ff_init_ff_cos_tabs(nbits);
103  s->tcos = ff_cos_tabs[nbits];
104  s->tsin = ff_cos_tabs[nbits] + (n >> 2);
105  s->rdft_calc = rdft_calc_c;
106 
107  if (ARCH_ARM) ff_rdft_init_arm(s);
108 
109  return 0;
110 }
111 
113 {
114  ff_fft_end(&s->fft);
115 }
av_cold void ff_rdft_end(RDFTContext *s)
Definition: rdft.c:112
const char * s
Definition: avisynth_c.h:768
int nbits
Definition: rdft.h:29
Definition: avfft.h:75
ptrdiff_t const GLvoid * data
Definition: opengl_enc.c:101
void(* fft_permute)(struct FFTContext *s, FFTComplex *z)
Do the permutation needed BEFORE calling fft_calc().
Definition: fft.h:101
FFTSample re
Definition: avfft.h:38
av_cold void ff_rdft_init_arm(RDFTContext *s)
Definition: rdft_init_arm.c:27
#define av_cold
Definition: attributes.h:82
#define RDFT_UNMANGLE(sign0, sign1)
RDFTransformType
Definition: avfft.h:71
#define AVERROR(e)
Definition: error.h:43
Definition: avfft.h:73
float FFTSample
Definition: avfft.h:35
void(* rdft_calc)(struct RDFTContext *s, FFTSample *z)
Definition: rdft.h:38
#define ff_fft_init
Definition: fft.h:149
Definition: avfft.h:72
int n
Definition: avisynth_c.h:684
Definition: avfft.h:74
const FFTSample * tsin
Definition: rdft.h:35
FFTContext fft
Definition: rdft.h:37
int negative_sin
Definition: rdft.h:36
#define ff_init_ff_cos_tabs
Definition: fft.h:141
#define ff_fft_end
Definition: fft.h:150
void(* fft_calc)(struct FFTContext *s, FFTComplex *z)
Do a complex FFT with the parameters defined in ff_fft_init().
Definition: fft.h:106
int sign_convention
Definition: rdft.h:31
static void rdft_calc_c(RDFTContext *s, FFTSample *data)
Map one real FFT into two parallel real even and odd FFTs.
Definition: rdft.c:35
int inverse
Definition: rdft.h:30
av_cold int ff_rdft_init(RDFTContext *s, int nbits, enum RDFTransformType trans)
Set up a real FFT.
Definition: rdft.c:86
const FFTSample * tcos
Definition: rdft.h:34