FFmpeg  4.0
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, odsum;
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[i2 ]-data[i1 ]); \
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  odsum.re = od.re*tcos[i] sign0 od.im*tsin[i]; \
66  odsum.im = od.im*tcos[i] sign1 od.re*tsin[i]; \
67  data[i1 ] = ev.re + odsum.re; \
68  data[i1+1] = ev.im + odsum.im; \
69  data[i2 ] = ev.re - odsum.re; \
70  data[i2+1] = odsum.im - ev.im; \
71  }
72 
73  if (s->negative_sin) {
74  RDFT_UNMANGLE(+,-)
75  } else {
76  RDFT_UNMANGLE(-,+)
77  }
78 
79  data[2*i+1]=s->sign_convention*data[2*i+1];
80  if (s->inverse) {
81  data[0] *= k1;
82  data[1] *= k1;
83  s->fft.fft_permute(&s->fft, (FFTComplex*)data);
84  s->fft.fft_calc(&s->fft, (FFTComplex*)data);
85  }
86 }
87 
88 av_cold int ff_rdft_init(RDFTContext *s, int nbits, enum RDFTransformType trans)
89 {
90  int n = 1 << nbits;
91  int ret;
92 
93  s->nbits = nbits;
94  s->inverse = trans == IDFT_C2R || trans == DFT_C2R;
95  s->sign_convention = trans == IDFT_R2C || trans == DFT_C2R ? 1 : -1;
96  s->negative_sin = trans == DFT_C2R || trans == DFT_R2C;
97 
98  if (nbits < 4 || nbits > 16)
99  return AVERROR(EINVAL);
100 
101  if ((ret = ff_fft_init(&s->fft, nbits-1, trans == IDFT_C2R || trans == IDFT_R2C)) < 0)
102  return ret;
103 
104  ff_init_ff_cos_tabs(nbits);
105  s->tcos = ff_cos_tabs[nbits];
106  s->tsin = ff_cos_tabs[nbits] + (n >> 2);
107  s->rdft_calc = rdft_calc_c;
108 
109  if (ARCH_ARM) ff_rdft_init_arm(s);
110 
111  return 0;
112 }
113 
115 {
116  ff_fft_end(&s->fft);
117 }
av_cold void ff_rdft_end(RDFTContext *s)
Definition: rdft.c:114
const char * s
Definition: avisynth_c.h:768
int nbits
Definition: rdft.h:29
Definition: avfft.h:75
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)
const char data[16]
Definition: mxf.c:90
RDFTransformType
Definition: avfft.h:71
void(* fft_permute)(struct FFTContext *s, FFTComplex *z)
Do the permutation needed BEFORE calling fft_calc().
Definition: fft.h:101
#define AVERROR(e)
Definition: error.h:43
Definition: avfft.h:73
float FFTSample
Definition: avfft.h:35
#define ARCH_ARM
Definition: config.h:19
#define ff_fft_init
Definition: fft.h:149
Definition: avfft.h:72
int n
Definition: avisynth_c.h:684
void(* rdft_calc)(struct RDFTContext *s, FFTSample *z)
Definition: rdft.h:38
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:88
const FFTSample * tcos
Definition: rdft.h:34