mdct.c 5.74 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
/* (C) 2008 Jean-Marc Valin, CSIRO
*/
/*
   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:
   
   - Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
   
   - Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
   
   - Neither the name of the Xiph.org Foundation nor the names of its
   contributors may be used to endorse or promote products derived from
   this software without specific prior written permission.
   
   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/* This is a simple MDCT implementation that uses a N/4 complex FFT
   to do most of the work. It should be relatively straightforward to
   plug in pretty much and FFT here.
   
   This replaces the Vorbis FFT (and uses the exact same API), which 
   was a bit too messy and that was ending up duplicating code 
   (might as well use the same FFT everywhere).
   
   The algorithm is similar to (and inspired from) Fabrice Bellard's
   MDCT implementation in FFMPEG, but has differences in signs, ordering
   and scaling in many places. 
*/
44

45
46
47
48
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

49
#include "mdct.h"
50
51
#include "kiss_fft.h"
#include <math.h>
52
#include "os_support.h"
Jean-Marc Valin's avatar
Jean-Marc Valin committed
53
#include "_kiss_fft_guts.h"
54
#include "mathops.h"
55

56
#ifndef M_PI
57
#define M_PI 3.141592653
58
59
#endif

60
61
62
63
64
65
66
void mdct_init(mdct_lookup *l,int N)
{
   int i;
   int N2, N4;
   l->n = N;
   N2 = N/2;
   N4 = N/4;
67
   l->kfft = kiss_fft_alloc(N4, NULL, NULL);
68
   l->trig = (kiss_twiddle_scalar*)celt_alloc(N2*sizeof(kiss_twiddle_scalar));
69
   /* We have enough points that sine isn't necessary */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
70
71
72
73
74
75
#if defined(FIXED_POINT)
#if defined(DOUBLE_PRECISION) & !defined(MIXED_PRECISION)
   for (i=0;i<N2;i++)
      l->trig[i] = SAMP_MAX*cos(2*M_PI*(i+1./8.)/N);
#else
   for (i=0;i<N2;i++)
76
      l->trig[i] = TRIG_UPSCALE*celt_cos_norm(DIV32(ADD32(SHL32(EXTEND32(i),17),16386),N));
Jean-Marc Valin's avatar
Jean-Marc Valin committed
77
78
#endif
#else
79
80
   for (i=0;i<N2;i++)
      l->trig[i] = cos(2*M_PI*(i+1./8.)/N);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
81
#endif
82
83
}

84
85
86
87
void mdct_clear(mdct_lookup *l)
{
   kiss_fft_free(l->kfft);
   celt_free(l->trig);
88
89
}

90
void mdct_forward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *out)
91
92
{
   int i;
93
   int N, N2, N4;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
94
   VARDECL(kiss_fft_scalar *f);
95
   SAVE_STACK;
96
97
98
   N = l->n;
   N2 = N/2;
   N4 = N/4;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
99
   ALLOC(f, N2, kiss_fft_scalar);
100
101
102
   
   /* Consider the input to be compused of four blocks: [a, b, c, d] */
   /* Shuffle, fold, pre-rotate (part 1) */
103
   for(i=0;i<N/8;i++)
104
   {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
105
      kiss_fft_scalar re, im;
106
      /* Real part arranged as -d-cR, Imag part arranged as -b+aR*/
107
108
      re = -HALF32(in[N2+N4+2*i] + in[N2+N4-2*i-1]);
      im = -HALF32(in[N4+2*i]    - in[N4-2*i-1]);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
109
110
      out[2*i]   = S_MUL(re,l->trig[i])  -  S_MUL(im,l->trig[i+N4]);
      out[2*i+1] = S_MUL(im,l->trig[i])  +  S_MUL(re,l->trig[i+N4]);
111
112
113
   }
   for(;i<N4;i++)
   {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
114
      kiss_fft_scalar re, im;
115
      /* Real part arranged as a-bR, Imag part arranged as -c-dR */
116
117
      re =  HALF32(in[2*i-N4] - in[N2+N4-2*i-1]);
      im = -HALF32(in[N4+2*i] + in[N+N4-2*i-1]);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
118
119
      out[2*i]   = S_MUL(re,l->trig[i])  -  S_MUL(im,l->trig[i+N4]);
      out[2*i+1] = S_MUL(im,l->trig[i])  +  S_MUL(re,l->trig[i+N4]);
120
121
122
123
124
125
126
127
   }

   /* N/4 complex FFT, which should normally down-scale by 4/N (but doesn't now) */
   kiss_fft(l->kfft, (const kiss_fft_cpx *)out, (kiss_fft_cpx *)f);

   /* Post-rotate and apply the scaling if the FFT doesn't to it itself */
   for(i=0;i<N4;i++)
   {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
128
129
      out[2*i]      = -S_MUL(f[2*i+1],l->trig[i+N4]) + S_MUL(f[2*i]  ,l->trig[i]);
      out[N2-1-2*i] = -S_MUL(f[2*i]  ,l->trig[i+N4]) - S_MUL(f[2*i+1],l->trig[i]);
130
   }
131
   RESTORE_STACK;
132
133
134
}


135
void mdct_backward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar *out)
136
137
{
   int i;
138
   int N, N2, N4;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
139
   VARDECL(kiss_fft_scalar *f);
140
   SAVE_STACK;
141
142
143
   N = l->n;
   N2 = N/2;
   N4 = N/4;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
144
   ALLOC(f, N2, kiss_fft_scalar);
145
146
147
148
   
   /* Pre-rotate */
   for(i=0;i<N4;i++) 
   {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
149
150
      out[2*i]   = -S_MUL(in[N2-2*i-1], l->trig[i])    - S_MUL(in[2*i],l->trig[i+N4]);
      out[2*i+1] =  S_MUL(in[N2-2*i-1], l->trig[i+N4]) - S_MUL(in[2*i],l->trig[i]);
151
152
153
   }

   /* Inverse N/4 complex FFT. This one should *not* downscale even in fixed-point */
154
   kiss_ifft(l->kfft, (const kiss_fft_cpx *)out, (kiss_fft_cpx *)f);
155
156
157
158
   
   /* Post-rotate */
   for(i=0;i<N4;i++)
   {
Jean-Marc Valin's avatar
Jean-Marc Valin committed
159
      kiss_fft_scalar re, im;
160
161
      re = f[2*i];
      im = f[2*i+1];
162
      /* We'd scale up by 2 here, but instead it's done when mixing the windows */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
163
164
      f[2*i]   = S_MUL(re,l->trig[i]) + S_MUL(im,l->trig[i+N4]);
      f[2*i+1] = S_MUL(im,l->trig[i]) - S_MUL(re,l->trig[i+N4]);
165
166
167
168
169
170
171
172
173
174
175
176
177
178
   }
   /* De-shuffle the components for the middle of the window only */
   for(i = 0; i < N4; i++)
   {
      out[N4+2*i]   =-f[2*i];
      out[N4+2*i+1] = f[N2-2*i-1];
   }

   /* Mirror on both sides for TDAC */
   for(i = 0; i < N4; i++)
   {
      out[i]     =-out[N2-i-1];
      out[N-i-1] = out[N2+i];
   }
179
   RESTORE_STACK;
180
181
182
}