rdft.c 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117
  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. * @file
  27. * (Inverse) Real Discrete Fourier Transforms.
  28. */
  29. /** Map one real FFT into two parallel real even and odd FFTs. Then interleave
  30. * the two real FFTs into one complex FFT. Unmangle the results.
  31. * ref: http://www.engineeringproductivitytools.com/stuff/T0001/PT10.HTM
  32. */
  33. static void rdft_calc_c(RDFTContext *s, FFTSample *data)
  34. {
  35. int i, i1, i2;
  36. FFTComplex ev, od, odsum;
  37. const int n = 1 << s->nbits;
  38. const float k1 = 0.5;
  39. const float k2 = 0.5 - s->inverse;
  40. const FFTSample *tcos = s->tcos;
  41. const FFTSample *tsin = s->tsin;
  42. if (!s->inverse) {
  43. s->fft.fft_permute(&s->fft, (FFTComplex*)data);
  44. s->fft.fft_calc(&s->fft, (FFTComplex*)data);
  45. }
  46. /* i=0 is a special case because of packing, the DC term is real, so we
  47. are going to throw the N/2 term (also real) in with it. */
  48. ev.re = data[0];
  49. data[0] = ev.re+data[1];
  50. data[1] = ev.re-data[1];
  51. #define RDFT_UNMANGLE(sign0, sign1) \
  52. for (i = 1; i < (n>>2); i++) { \
  53. i1 = 2*i; \
  54. i2 = n-i1; \
  55. /* Separate even and odd FFTs */ \
  56. ev.re = k1*(data[i1 ]+data[i2 ]); \
  57. od.im = k2*(data[i2 ]-data[i1 ]); \
  58. ev.im = k1*(data[i1+1]-data[i2+1]); \
  59. od.re = k2*(data[i1+1]+data[i2+1]); \
  60. /* Apply twiddle factors to the odd FFT and add to the even FFT */ \
  61. odsum.re = od.re*tcos[i] sign0 od.im*tsin[i]; \
  62. odsum.im = od.im*tcos[i] sign1 od.re*tsin[i]; \
  63. data[i1 ] = ev.re + odsum.re; \
  64. data[i1+1] = ev.im + odsum.im; \
  65. data[i2 ] = ev.re - odsum.re; \
  66. data[i2+1] = odsum.im - ev.im; \
  67. }
  68. if (s->negative_sin) {
  69. RDFT_UNMANGLE(+,-)
  70. } else {
  71. RDFT_UNMANGLE(-,+)
  72. }
  73. data[2*i+1]=s->sign_convention*data[2*i+1];
  74. if (s->inverse) {
  75. data[0] *= k1;
  76. data[1] *= k1;
  77. s->fft.fft_permute(&s->fft, (FFTComplex*)data);
  78. s->fft.fft_calc(&s->fft, (FFTComplex*)data);
  79. }
  80. }
  81. av_cold int ff_rdft_init(RDFTContext *s, int nbits, enum RDFTransformType trans)
  82. {
  83. int n = 1 << nbits;
  84. int ret;
  85. s->nbits = nbits;
  86. s->inverse = trans == IDFT_C2R || trans == DFT_C2R;
  87. s->sign_convention = trans == IDFT_R2C || trans == DFT_C2R ? 1 : -1;
  88. s->negative_sin = trans == DFT_C2R || trans == DFT_R2C;
  89. if (nbits < 4 || nbits > 16)
  90. return AVERROR(EINVAL);
  91. if ((ret = ff_fft_init(&s->fft, nbits-1, trans == IDFT_C2R || trans == IDFT_R2C)) < 0)
  92. return ret;
  93. ff_init_ff_cos_tabs(nbits);
  94. s->tcos = ff_cos_tabs[nbits];
  95. s->tsin = ff_cos_tabs[nbits] + (n >> 2);
  96. s->rdft_calc = rdft_calc_c;
  97. if (ARCH_ARM) ff_rdft_init_arm(s);
  98. return 0;
  99. }
  100. av_cold void ff_rdft_end(RDFTContext *s)
  101. {
  102. ff_fft_end(&s->fft);
  103. }