[029bf4e] | 1 | // modifications made for aubio: |
---|
[5c6b264] | 2 | // - replace all 'double' with 'smpl_t' |
---|
| 3 | // - include "aubio_priv.h" (for config.h and types.h) |
---|
[029bf4e] | 4 | // - add missing prototypes |
---|
[5cb8abe] | 5 | // - use COS and SIN macros |
---|
[392ad1c] | 6 | // - declare initialization as static |
---|
| 7 | // - prefix public function with aubio_ooura_ |
---|
[5c6b264] | 8 | |
---|
| 9 | #include "aubio_priv.h" |
---|
| 10 | |
---|
[392ad1c] | 11 | void aubio_ooura_cdft(int n, int isgn, smpl_t *a, int *ip, smpl_t *w); |
---|
| 12 | void aubio_ooura_rdft(int n, int isgn, smpl_t *a, int *ip, smpl_t *w); |
---|
| 13 | void aubio_ooura_ddct(int n, int isgn, smpl_t *a, int *ip, smpl_t *w); |
---|
| 14 | void aubio_ooura_ddst(int n, int isgn, smpl_t *a, int *ip, smpl_t *w); |
---|
| 15 | void aubio_ooura_dfct(int n, smpl_t *a, smpl_t *t, int *ip, smpl_t *w); |
---|
| 16 | void aubio_ooura_dfst(int n, smpl_t *a, smpl_t *t, int *ip, smpl_t *w); |
---|
| 17 | static void makewt(int nw, int *ip, smpl_t *w); |
---|
| 18 | static void makect(int nc, int *ip, smpl_t *c); |
---|
| 19 | static void bitrv2(int n, int *ip, smpl_t *a); |
---|
| 20 | static void bitrv2conj(int n, int *ip, smpl_t *a); |
---|
| 21 | static void cftfsub(int n, smpl_t *a, smpl_t *w); |
---|
| 22 | static void cftbsub(int n, smpl_t *a, smpl_t *w); |
---|
| 23 | static void cft1st(int n, smpl_t *a, smpl_t *w); |
---|
| 24 | static void cftmdl(int n, int l, smpl_t *a, smpl_t *w); |
---|
| 25 | static void rftfsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
| 26 | static void rftbsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
| 27 | static void dctsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
| 28 | static void dstsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
[029bf4e] | 29 | |
---|
[cfaa3c4] | 30 | /* |
---|
| 31 | Fast Fourier/Cosine/Sine Transform |
---|
| 32 | dimension :one |
---|
| 33 | data length :power of 2 |
---|
| 34 | decimation :frequency |
---|
| 35 | radix :8, 4, 2 |
---|
| 36 | data :inplace |
---|
| 37 | table :use |
---|
| 38 | functions |
---|
| 39 | cdft: Complex Discrete Fourier Transform |
---|
| 40 | rdft: Real Discrete Fourier Transform |
---|
| 41 | ddct: Discrete Cosine Transform |
---|
| 42 | ddst: Discrete Sine Transform |
---|
| 43 | dfct: Cosine Transform of RDFT (Real Symmetric DFT) |
---|
| 44 | dfst: Sine Transform of RDFT (Real Anti-symmetric DFT) |
---|
| 45 | function prototypes |
---|
[5c6b264] | 46 | void cdft(int, int, smpl_t *, int *, smpl_t *); |
---|
| 47 | void rdft(int, int, smpl_t *, int *, smpl_t *); |
---|
| 48 | void ddct(int, int, smpl_t *, int *, smpl_t *); |
---|
| 49 | void ddst(int, int, smpl_t *, int *, smpl_t *); |
---|
| 50 | void dfct(int, smpl_t *, smpl_t *, int *, smpl_t *); |
---|
| 51 | void dfst(int, smpl_t *, smpl_t *, int *, smpl_t *); |
---|
[cfaa3c4] | 52 | |
---|
| 53 | |
---|
| 54 | -------- Complex DFT (Discrete Fourier Transform) -------- |
---|
| 55 | [definition] |
---|
| 56 | <case1> |
---|
| 57 | X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n |
---|
| 58 | <case2> |
---|
| 59 | X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n |
---|
| 60 | (notes: sum_j=0^n-1 is a summation from j=0 to n-1) |
---|
| 61 | [usage] |
---|
| 62 | <case1> |
---|
| 63 | ip[0] = 0; // first time only |
---|
| 64 | cdft(2*n, 1, a, ip, w); |
---|
| 65 | <case2> |
---|
| 66 | ip[0] = 0; // first time only |
---|
| 67 | cdft(2*n, -1, a, ip, w); |
---|
| 68 | [parameters] |
---|
| 69 | 2*n :data length (int) |
---|
| 70 | n >= 1, n = power of 2 |
---|
[5c6b264] | 71 | a[0...2*n-1] :input/output data (smpl_t *) |
---|
[cfaa3c4] | 72 | input data |
---|
| 73 | a[2*j] = Re(x[j]), |
---|
| 74 | a[2*j+1] = Im(x[j]), 0<=j<n |
---|
| 75 | output data |
---|
| 76 | a[2*k] = Re(X[k]), |
---|
| 77 | a[2*k+1] = Im(X[k]), 0<=k<n |
---|
| 78 | ip[0...*] :work area for bit reversal (int *) |
---|
| 79 | length of ip >= 2+sqrt(n) |
---|
| 80 | strictly, |
---|
| 81 | length of ip >= |
---|
| 82 | 2+(1<<(int)(log(n+0.5)/log(2))/2). |
---|
| 83 | ip[0],ip[1] are pointers of the cos/sin table. |
---|
[5c6b264] | 84 | w[0...n/2-1] :cos/sin table (smpl_t *) |
---|
[cfaa3c4] | 85 | w[],ip[] are initialized if ip[0] == 0. |
---|
| 86 | [remark] |
---|
| 87 | Inverse of |
---|
| 88 | cdft(2*n, -1, a, ip, w); |
---|
| 89 | is |
---|
| 90 | cdft(2*n, 1, a, ip, w); |
---|
| 91 | for (j = 0; j <= 2 * n - 1; j++) { |
---|
| 92 | a[j] *= 1.0 / n; |
---|
| 93 | } |
---|
| 94 | . |
---|
| 95 | |
---|
| 96 | |
---|
| 97 | -------- Real DFT / Inverse of Real DFT -------- |
---|
| 98 | [definition] |
---|
| 99 | <case1> RDFT |
---|
| 100 | R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2 |
---|
| 101 | I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2 |
---|
| 102 | <case2> IRDFT (excluding scale) |
---|
| 103 | a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + |
---|
| 104 | sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + |
---|
| 105 | sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n |
---|
| 106 | [usage] |
---|
| 107 | <case1> |
---|
| 108 | ip[0] = 0; // first time only |
---|
| 109 | rdft(n, 1, a, ip, w); |
---|
| 110 | <case2> |
---|
| 111 | ip[0] = 0; // first time only |
---|
| 112 | rdft(n, -1, a, ip, w); |
---|
| 113 | [parameters] |
---|
| 114 | n :data length (int) |
---|
| 115 | n >= 2, n = power of 2 |
---|
[5c6b264] | 116 | a[0...n-1] :input/output data (smpl_t *) |
---|
[cfaa3c4] | 117 | <case1> |
---|
| 118 | output data |
---|
| 119 | a[2*k] = R[k], 0<=k<n/2 |
---|
| 120 | a[2*k+1] = I[k], 0<k<n/2 |
---|
| 121 | a[1] = R[n/2] |
---|
| 122 | <case2> |
---|
| 123 | input data |
---|
| 124 | a[2*j] = R[j], 0<=j<n/2 |
---|
| 125 | a[2*j+1] = I[j], 0<j<n/2 |
---|
| 126 | a[1] = R[n/2] |
---|
| 127 | ip[0...*] :work area for bit reversal (int *) |
---|
| 128 | length of ip >= 2+sqrt(n/2) |
---|
| 129 | strictly, |
---|
| 130 | length of ip >= |
---|
| 131 | 2+(1<<(int)(log(n/2+0.5)/log(2))/2). |
---|
| 132 | ip[0],ip[1] are pointers of the cos/sin table. |
---|
[5c6b264] | 133 | w[0...n/2-1] :cos/sin table (smpl_t *) |
---|
[cfaa3c4] | 134 | w[],ip[] are initialized if ip[0] == 0. |
---|
| 135 | [remark] |
---|
| 136 | Inverse of |
---|
| 137 | rdft(n, 1, a, ip, w); |
---|
| 138 | is |
---|
| 139 | rdft(n, -1, a, ip, w); |
---|
| 140 | for (j = 0; j <= n - 1; j++) { |
---|
| 141 | a[j] *= 2.0 / n; |
---|
| 142 | } |
---|
| 143 | . |
---|
| 144 | |
---|
| 145 | |
---|
| 146 | -------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- |
---|
| 147 | [definition] |
---|
| 148 | <case1> IDCT (excluding scale) |
---|
| 149 | C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n |
---|
| 150 | <case2> DCT |
---|
| 151 | C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n |
---|
| 152 | [usage] |
---|
| 153 | <case1> |
---|
| 154 | ip[0] = 0; // first time only |
---|
| 155 | ddct(n, 1, a, ip, w); |
---|
| 156 | <case2> |
---|
| 157 | ip[0] = 0; // first time only |
---|
| 158 | ddct(n, -1, a, ip, w); |
---|
| 159 | [parameters] |
---|
| 160 | n :data length (int) |
---|
| 161 | n >= 2, n = power of 2 |
---|
[5c6b264] | 162 | a[0...n-1] :input/output data (smpl_t *) |
---|
[cfaa3c4] | 163 | output data |
---|
| 164 | a[k] = C[k], 0<=k<n |
---|
| 165 | ip[0...*] :work area for bit reversal (int *) |
---|
| 166 | length of ip >= 2+sqrt(n/2) |
---|
| 167 | strictly, |
---|
| 168 | length of ip >= |
---|
| 169 | 2+(1<<(int)(log(n/2+0.5)/log(2))/2). |
---|
| 170 | ip[0],ip[1] are pointers of the cos/sin table. |
---|
[5c6b264] | 171 | w[0...n*5/4-1] :cos/sin table (smpl_t *) |
---|
[cfaa3c4] | 172 | w[],ip[] are initialized if ip[0] == 0. |
---|
| 173 | [remark] |
---|
| 174 | Inverse of |
---|
| 175 | ddct(n, -1, a, ip, w); |
---|
| 176 | is |
---|
| 177 | a[0] *= 0.5; |
---|
| 178 | ddct(n, 1, a, ip, w); |
---|
| 179 | for (j = 0; j <= n - 1; j++) { |
---|
| 180 | a[j] *= 2.0 / n; |
---|
| 181 | } |
---|
| 182 | . |
---|
| 183 | |
---|
| 184 | |
---|
| 185 | -------- DST (Discrete Sine Transform) / Inverse of DST -------- |
---|
| 186 | [definition] |
---|
| 187 | <case1> IDST (excluding scale) |
---|
| 188 | S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n |
---|
| 189 | <case2> DST |
---|
| 190 | S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n |
---|
| 191 | [usage] |
---|
| 192 | <case1> |
---|
| 193 | ip[0] = 0; // first time only |
---|
| 194 | ddst(n, 1, a, ip, w); |
---|
| 195 | <case2> |
---|
| 196 | ip[0] = 0; // first time only |
---|
| 197 | ddst(n, -1, a, ip, w); |
---|
| 198 | [parameters] |
---|
| 199 | n :data length (int) |
---|
| 200 | n >= 2, n = power of 2 |
---|
[5c6b264] | 201 | a[0...n-1] :input/output data (smpl_t *) |
---|
[cfaa3c4] | 202 | <case1> |
---|
| 203 | input data |
---|
| 204 | a[j] = A[j], 0<j<n |
---|
| 205 | a[0] = A[n] |
---|
| 206 | output data |
---|
| 207 | a[k] = S[k], 0<=k<n |
---|
| 208 | <case2> |
---|
| 209 | output data |
---|
| 210 | a[k] = S[k], 0<k<n |
---|
| 211 | a[0] = S[n] |
---|
| 212 | ip[0...*] :work area for bit reversal (int *) |
---|
| 213 | length of ip >= 2+sqrt(n/2) |
---|
| 214 | strictly, |
---|
| 215 | length of ip >= |
---|
| 216 | 2+(1<<(int)(log(n/2+0.5)/log(2))/2). |
---|
| 217 | ip[0],ip[1] are pointers of the cos/sin table. |
---|
[5c6b264] | 218 | w[0...n*5/4-1] :cos/sin table (smpl_t *) |
---|
[cfaa3c4] | 219 | w[],ip[] are initialized if ip[0] == 0. |
---|
| 220 | [remark] |
---|
| 221 | Inverse of |
---|
| 222 | ddst(n, -1, a, ip, w); |
---|
| 223 | is |
---|
| 224 | a[0] *= 0.5; |
---|
| 225 | ddst(n, 1, a, ip, w); |
---|
| 226 | for (j = 0; j <= n - 1; j++) { |
---|
| 227 | a[j] *= 2.0 / n; |
---|
| 228 | } |
---|
| 229 | . |
---|
| 230 | |
---|
| 231 | |
---|
| 232 | -------- Cosine Transform of RDFT (Real Symmetric DFT) -------- |
---|
| 233 | [definition] |
---|
| 234 | C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n |
---|
| 235 | [usage] |
---|
| 236 | ip[0] = 0; // first time only |
---|
| 237 | dfct(n, a, t, ip, w); |
---|
| 238 | [parameters] |
---|
| 239 | n :data length - 1 (int) |
---|
| 240 | n >= 2, n = power of 2 |
---|
[5c6b264] | 241 | a[0...n] :input/output data (smpl_t *) |
---|
[cfaa3c4] | 242 | output data |
---|
| 243 | a[k] = C[k], 0<=k<=n |
---|
[5c6b264] | 244 | t[0...n/2] :work area (smpl_t *) |
---|
[cfaa3c4] | 245 | ip[0...*] :work area for bit reversal (int *) |
---|
| 246 | length of ip >= 2+sqrt(n/4) |
---|
| 247 | strictly, |
---|
| 248 | length of ip >= |
---|
| 249 | 2+(1<<(int)(log(n/4+0.5)/log(2))/2). |
---|
| 250 | ip[0],ip[1] are pointers of the cos/sin table. |
---|
[5c6b264] | 251 | w[0...n*5/8-1] :cos/sin table (smpl_t *) |
---|
[cfaa3c4] | 252 | w[],ip[] are initialized if ip[0] == 0. |
---|
| 253 | [remark] |
---|
| 254 | Inverse of |
---|
| 255 | a[0] *= 0.5; |
---|
| 256 | a[n] *= 0.5; |
---|
| 257 | dfct(n, a, t, ip, w); |
---|
| 258 | is |
---|
| 259 | a[0] *= 0.5; |
---|
| 260 | a[n] *= 0.5; |
---|
| 261 | dfct(n, a, t, ip, w); |
---|
| 262 | for (j = 0; j <= n; j++) { |
---|
| 263 | a[j] *= 2.0 / n; |
---|
| 264 | } |
---|
| 265 | . |
---|
| 266 | |
---|
| 267 | |
---|
| 268 | -------- Sine Transform of RDFT (Real Anti-symmetric DFT) -------- |
---|
| 269 | [definition] |
---|
| 270 | S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n |
---|
| 271 | [usage] |
---|
| 272 | ip[0] = 0; // first time only |
---|
| 273 | dfst(n, a, t, ip, w); |
---|
| 274 | [parameters] |
---|
| 275 | n :data length + 1 (int) |
---|
| 276 | n >= 2, n = power of 2 |
---|
[5c6b264] | 277 | a[0...n-1] :input/output data (smpl_t *) |
---|
[cfaa3c4] | 278 | output data |
---|
| 279 | a[k] = S[k], 0<k<n |
---|
| 280 | (a[0] is used for work area) |
---|
[5c6b264] | 281 | t[0...n/2-1] :work area (smpl_t *) |
---|
[cfaa3c4] | 282 | ip[0...*] :work area for bit reversal (int *) |
---|
| 283 | length of ip >= 2+sqrt(n/4) |
---|
| 284 | strictly, |
---|
| 285 | length of ip >= |
---|
| 286 | 2+(1<<(int)(log(n/4+0.5)/log(2))/2). |
---|
| 287 | ip[0],ip[1] are pointers of the cos/sin table. |
---|
[5c6b264] | 288 | w[0...n*5/8-1] :cos/sin table (smpl_t *) |
---|
[cfaa3c4] | 289 | w[],ip[] are initialized if ip[0] == 0. |
---|
| 290 | [remark] |
---|
| 291 | Inverse of |
---|
| 292 | dfst(n, a, t, ip, w); |
---|
| 293 | is |
---|
| 294 | dfst(n, a, t, ip, w); |
---|
| 295 | for (j = 1; j <= n - 1; j++) { |
---|
| 296 | a[j] *= 2.0 / n; |
---|
| 297 | } |
---|
| 298 | . |
---|
| 299 | |
---|
| 300 | |
---|
| 301 | Appendix : |
---|
| 302 | The cos/sin table is recalculated when the larger table required. |
---|
| 303 | w[] and ip[] are compatible with all routines. |
---|
| 304 | */ |
---|
| 305 | |
---|
| 306 | |
---|
[392ad1c] | 307 | void aubio_ooura_cdft(int n, int isgn, smpl_t *a, int *ip, smpl_t *w) |
---|
[cfaa3c4] | 308 | { |
---|
[5c6b264] | 309 | void makewt(int nw, int *ip, smpl_t *w); |
---|
| 310 | void bitrv2(int n, int *ip, smpl_t *a); |
---|
| 311 | void bitrv2conj(int n, int *ip, smpl_t *a); |
---|
| 312 | void cftfsub(int n, smpl_t *a, smpl_t *w); |
---|
| 313 | void cftbsub(int n, smpl_t *a, smpl_t *w); |
---|
[cfaa3c4] | 314 | |
---|
| 315 | if (n > (ip[0] << 2)) { |
---|
| 316 | makewt(n >> 2, ip, w); |
---|
| 317 | } |
---|
| 318 | if (n > 4) { |
---|
| 319 | if (isgn >= 0) { |
---|
| 320 | bitrv2(n, ip + 2, a); |
---|
| 321 | cftfsub(n, a, w); |
---|
| 322 | } else { |
---|
| 323 | bitrv2conj(n, ip + 2, a); |
---|
| 324 | cftbsub(n, a, w); |
---|
| 325 | } |
---|
| 326 | } else if (n == 4) { |
---|
| 327 | cftfsub(n, a, w); |
---|
| 328 | } |
---|
| 329 | } |
---|
| 330 | |
---|
| 331 | |
---|
[392ad1c] | 332 | void aubio_ooura_rdft(int n, int isgn, smpl_t *a, int *ip, smpl_t *w) |
---|
[cfaa3c4] | 333 | { |
---|
[5c6b264] | 334 | void makewt(int nw, int *ip, smpl_t *w); |
---|
| 335 | void makect(int nc, int *ip, smpl_t *c); |
---|
| 336 | void bitrv2(int n, int *ip, smpl_t *a); |
---|
| 337 | void cftfsub(int n, smpl_t *a, smpl_t *w); |
---|
| 338 | void cftbsub(int n, smpl_t *a, smpl_t *w); |
---|
| 339 | void rftfsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
| 340 | void rftbsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
[cfaa3c4] | 341 | int nw, nc; |
---|
[5c6b264] | 342 | smpl_t xi; |
---|
[cfaa3c4] | 343 | |
---|
| 344 | nw = ip[0]; |
---|
| 345 | if (n > (nw << 2)) { |
---|
| 346 | nw = n >> 2; |
---|
| 347 | makewt(nw, ip, w); |
---|
| 348 | } |
---|
| 349 | nc = ip[1]; |
---|
| 350 | if (n > (nc << 2)) { |
---|
| 351 | nc = n >> 2; |
---|
| 352 | makect(nc, ip, w + nw); |
---|
| 353 | } |
---|
| 354 | if (isgn >= 0) { |
---|
| 355 | if (n > 4) { |
---|
| 356 | bitrv2(n, ip + 2, a); |
---|
| 357 | cftfsub(n, a, w); |
---|
| 358 | rftfsub(n, a, nc, w + nw); |
---|
| 359 | } else if (n == 4) { |
---|
| 360 | cftfsub(n, a, w); |
---|
| 361 | } |
---|
| 362 | xi = a[0] - a[1]; |
---|
| 363 | a[0] += a[1]; |
---|
| 364 | a[1] = xi; |
---|
| 365 | } else { |
---|
| 366 | a[1] = 0.5 * (a[0] - a[1]); |
---|
| 367 | a[0] -= a[1]; |
---|
| 368 | if (n > 4) { |
---|
| 369 | rftbsub(n, a, nc, w + nw); |
---|
| 370 | bitrv2(n, ip + 2, a); |
---|
| 371 | cftbsub(n, a, w); |
---|
| 372 | } else if (n == 4) { |
---|
| 373 | cftfsub(n, a, w); |
---|
| 374 | } |
---|
| 375 | } |
---|
| 376 | } |
---|
| 377 | |
---|
| 378 | |
---|
[392ad1c] | 379 | void aubio_ooura_ddct(int n, int isgn, smpl_t *a, int *ip, smpl_t *w) |
---|
[cfaa3c4] | 380 | { |
---|
[5c6b264] | 381 | void makewt(int nw, int *ip, smpl_t *w); |
---|
| 382 | void makect(int nc, int *ip, smpl_t *c); |
---|
| 383 | void bitrv2(int n, int *ip, smpl_t *a); |
---|
| 384 | void cftfsub(int n, smpl_t *a, smpl_t *w); |
---|
| 385 | void cftbsub(int n, smpl_t *a, smpl_t *w); |
---|
| 386 | void rftfsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
| 387 | void rftbsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
| 388 | void dctsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
[cfaa3c4] | 389 | int j, nw, nc; |
---|
[5c6b264] | 390 | smpl_t xr; |
---|
[cfaa3c4] | 391 | |
---|
| 392 | nw = ip[0]; |
---|
| 393 | if (n > (nw << 2)) { |
---|
| 394 | nw = n >> 2; |
---|
| 395 | makewt(nw, ip, w); |
---|
| 396 | } |
---|
| 397 | nc = ip[1]; |
---|
| 398 | if (n > nc) { |
---|
| 399 | nc = n; |
---|
| 400 | makect(nc, ip, w + nw); |
---|
| 401 | } |
---|
| 402 | if (isgn < 0) { |
---|
| 403 | xr = a[n - 1]; |
---|
| 404 | for (j = n - 2; j >= 2; j -= 2) { |
---|
| 405 | a[j + 1] = a[j] - a[j - 1]; |
---|
| 406 | a[j] += a[j - 1]; |
---|
| 407 | } |
---|
| 408 | a[1] = a[0] - xr; |
---|
| 409 | a[0] += xr; |
---|
| 410 | if (n > 4) { |
---|
| 411 | rftbsub(n, a, nc, w + nw); |
---|
| 412 | bitrv2(n, ip + 2, a); |
---|
| 413 | cftbsub(n, a, w); |
---|
| 414 | } else if (n == 4) { |
---|
| 415 | cftfsub(n, a, w); |
---|
| 416 | } |
---|
| 417 | } |
---|
| 418 | dctsub(n, a, nc, w + nw); |
---|
| 419 | if (isgn >= 0) { |
---|
| 420 | if (n > 4) { |
---|
| 421 | bitrv2(n, ip + 2, a); |
---|
| 422 | cftfsub(n, a, w); |
---|
| 423 | rftfsub(n, a, nc, w + nw); |
---|
| 424 | } else if (n == 4) { |
---|
| 425 | cftfsub(n, a, w); |
---|
| 426 | } |
---|
| 427 | xr = a[0] - a[1]; |
---|
| 428 | a[0] += a[1]; |
---|
| 429 | for (j = 2; j < n; j += 2) { |
---|
| 430 | a[j - 1] = a[j] - a[j + 1]; |
---|
| 431 | a[j] += a[j + 1]; |
---|
| 432 | } |
---|
| 433 | a[n - 1] = xr; |
---|
| 434 | } |
---|
| 435 | } |
---|
| 436 | |
---|
| 437 | |
---|
[392ad1c] | 438 | void aubio_ooura_ddst(int n, int isgn, smpl_t *a, int *ip, smpl_t *w) |
---|
[cfaa3c4] | 439 | { |
---|
[5c6b264] | 440 | void makewt(int nw, int *ip, smpl_t *w); |
---|
| 441 | void makect(int nc, int *ip, smpl_t *c); |
---|
| 442 | void bitrv2(int n, int *ip, smpl_t *a); |
---|
| 443 | void cftfsub(int n, smpl_t *a, smpl_t *w); |
---|
| 444 | void cftbsub(int n, smpl_t *a, smpl_t *w); |
---|
| 445 | void rftfsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
| 446 | void rftbsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
| 447 | void dstsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
[cfaa3c4] | 448 | int j, nw, nc; |
---|
[5c6b264] | 449 | smpl_t xr; |
---|
[cfaa3c4] | 450 | |
---|
| 451 | nw = ip[0]; |
---|
| 452 | if (n > (nw << 2)) { |
---|
| 453 | nw = n >> 2; |
---|
| 454 | makewt(nw, ip, w); |
---|
| 455 | } |
---|
| 456 | nc = ip[1]; |
---|
| 457 | if (n > nc) { |
---|
| 458 | nc = n; |
---|
| 459 | makect(nc, ip, w + nw); |
---|
| 460 | } |
---|
| 461 | if (isgn < 0) { |
---|
| 462 | xr = a[n - 1]; |
---|
| 463 | for (j = n - 2; j >= 2; j -= 2) { |
---|
| 464 | a[j + 1] = -a[j] - a[j - 1]; |
---|
| 465 | a[j] -= a[j - 1]; |
---|
| 466 | } |
---|
| 467 | a[1] = a[0] + xr; |
---|
| 468 | a[0] -= xr; |
---|
| 469 | if (n > 4) { |
---|
| 470 | rftbsub(n, a, nc, w + nw); |
---|
| 471 | bitrv2(n, ip + 2, a); |
---|
| 472 | cftbsub(n, a, w); |
---|
| 473 | } else if (n == 4) { |
---|
| 474 | cftfsub(n, a, w); |
---|
| 475 | } |
---|
| 476 | } |
---|
| 477 | dstsub(n, a, nc, w + nw); |
---|
| 478 | if (isgn >= 0) { |
---|
| 479 | if (n > 4) { |
---|
| 480 | bitrv2(n, ip + 2, a); |
---|
| 481 | cftfsub(n, a, w); |
---|
| 482 | rftfsub(n, a, nc, w + nw); |
---|
| 483 | } else if (n == 4) { |
---|
| 484 | cftfsub(n, a, w); |
---|
| 485 | } |
---|
| 486 | xr = a[0] - a[1]; |
---|
| 487 | a[0] += a[1]; |
---|
| 488 | for (j = 2; j < n; j += 2) { |
---|
| 489 | a[j - 1] = -a[j] - a[j + 1]; |
---|
| 490 | a[j] -= a[j + 1]; |
---|
| 491 | } |
---|
| 492 | a[n - 1] = -xr; |
---|
| 493 | } |
---|
| 494 | } |
---|
| 495 | |
---|
| 496 | |
---|
[392ad1c] | 497 | void aubio_ooura_dfct(int n, smpl_t *a, smpl_t *t, int *ip, smpl_t *w) |
---|
[cfaa3c4] | 498 | { |
---|
[5c6b264] | 499 | void makewt(int nw, int *ip, smpl_t *w); |
---|
| 500 | void makect(int nc, int *ip, smpl_t *c); |
---|
| 501 | void bitrv2(int n, int *ip, smpl_t *a); |
---|
| 502 | void cftfsub(int n, smpl_t *a, smpl_t *w); |
---|
| 503 | void rftfsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
| 504 | void dctsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
[cfaa3c4] | 505 | int j, k, l, m, mh, nw, nc; |
---|
[5c6b264] | 506 | smpl_t xr, xi, yr, yi; |
---|
[cfaa3c4] | 507 | |
---|
| 508 | nw = ip[0]; |
---|
| 509 | if (n > (nw << 3)) { |
---|
| 510 | nw = n >> 3; |
---|
| 511 | makewt(nw, ip, w); |
---|
| 512 | } |
---|
| 513 | nc = ip[1]; |
---|
| 514 | if (n > (nc << 1)) { |
---|
| 515 | nc = n >> 1; |
---|
| 516 | makect(nc, ip, w + nw); |
---|
| 517 | } |
---|
| 518 | m = n >> 1; |
---|
| 519 | yi = a[m]; |
---|
| 520 | xi = a[0] + a[n]; |
---|
| 521 | a[0] -= a[n]; |
---|
| 522 | t[0] = xi - yi; |
---|
| 523 | t[m] = xi + yi; |
---|
| 524 | if (n > 2) { |
---|
| 525 | mh = m >> 1; |
---|
| 526 | for (j = 1; j < mh; j++) { |
---|
| 527 | k = m - j; |
---|
| 528 | xr = a[j] - a[n - j]; |
---|
| 529 | xi = a[j] + a[n - j]; |
---|
| 530 | yr = a[k] - a[n - k]; |
---|
| 531 | yi = a[k] + a[n - k]; |
---|
| 532 | a[j] = xr; |
---|
| 533 | a[k] = yr; |
---|
| 534 | t[j] = xi - yi; |
---|
| 535 | t[k] = xi + yi; |
---|
| 536 | } |
---|
| 537 | t[mh] = a[mh] + a[n - mh]; |
---|
| 538 | a[mh] -= a[n - mh]; |
---|
| 539 | dctsub(m, a, nc, w + nw); |
---|
| 540 | if (m > 4) { |
---|
| 541 | bitrv2(m, ip + 2, a); |
---|
| 542 | cftfsub(m, a, w); |
---|
| 543 | rftfsub(m, a, nc, w + nw); |
---|
| 544 | } else if (m == 4) { |
---|
| 545 | cftfsub(m, a, w); |
---|
| 546 | } |
---|
| 547 | a[n - 1] = a[0] - a[1]; |
---|
| 548 | a[1] = a[0] + a[1]; |
---|
| 549 | for (j = m - 2; j >= 2; j -= 2) { |
---|
| 550 | a[2 * j + 1] = a[j] + a[j + 1]; |
---|
| 551 | a[2 * j - 1] = a[j] - a[j + 1]; |
---|
| 552 | } |
---|
| 553 | l = 2; |
---|
| 554 | m = mh; |
---|
| 555 | while (m >= 2) { |
---|
| 556 | dctsub(m, t, nc, w + nw); |
---|
| 557 | if (m > 4) { |
---|
| 558 | bitrv2(m, ip + 2, t); |
---|
| 559 | cftfsub(m, t, w); |
---|
| 560 | rftfsub(m, t, nc, w + nw); |
---|
| 561 | } else if (m == 4) { |
---|
| 562 | cftfsub(m, t, w); |
---|
| 563 | } |
---|
| 564 | a[n - l] = t[0] - t[1]; |
---|
| 565 | a[l] = t[0] + t[1]; |
---|
| 566 | k = 0; |
---|
| 567 | for (j = 2; j < m; j += 2) { |
---|
| 568 | k += l << 2; |
---|
| 569 | a[k - l] = t[j] - t[j + 1]; |
---|
| 570 | a[k + l] = t[j] + t[j + 1]; |
---|
| 571 | } |
---|
| 572 | l <<= 1; |
---|
| 573 | mh = m >> 1; |
---|
| 574 | for (j = 0; j < mh; j++) { |
---|
| 575 | k = m - j; |
---|
| 576 | t[j] = t[m + k] - t[m + j]; |
---|
| 577 | t[k] = t[m + k] + t[m + j]; |
---|
| 578 | } |
---|
| 579 | t[mh] = t[m + mh]; |
---|
| 580 | m = mh; |
---|
| 581 | } |
---|
| 582 | a[l] = t[0]; |
---|
| 583 | a[n] = t[2] - t[1]; |
---|
| 584 | a[0] = t[2] + t[1]; |
---|
| 585 | } else { |
---|
| 586 | a[1] = a[0]; |
---|
| 587 | a[2] = t[0]; |
---|
| 588 | a[0] = t[1]; |
---|
| 589 | } |
---|
| 590 | } |
---|
| 591 | |
---|
| 592 | |
---|
[392ad1c] | 593 | void aubio_ooura_dfst(int n, smpl_t *a, smpl_t *t, int *ip, smpl_t *w) |
---|
[cfaa3c4] | 594 | { |
---|
[5c6b264] | 595 | void makewt(int nw, int *ip, smpl_t *w); |
---|
| 596 | void makect(int nc, int *ip, smpl_t *c); |
---|
| 597 | void bitrv2(int n, int *ip, smpl_t *a); |
---|
| 598 | void cftfsub(int n, smpl_t *a, smpl_t *w); |
---|
| 599 | void rftfsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
| 600 | void dstsub(int n, smpl_t *a, int nc, smpl_t *c); |
---|
[cfaa3c4] | 601 | int j, k, l, m, mh, nw, nc; |
---|
[5c6b264] | 602 | smpl_t xr, xi, yr, yi; |
---|
[cfaa3c4] | 603 | |
---|
| 604 | nw = ip[0]; |
---|
| 605 | if (n > (nw << 3)) { |
---|
| 606 | nw = n >> 3; |
---|
| 607 | makewt(nw, ip, w); |
---|
| 608 | } |
---|
| 609 | nc = ip[1]; |
---|
| 610 | if (n > (nc << 1)) { |
---|
| 611 | nc = n >> 1; |
---|
| 612 | makect(nc, ip, w + nw); |
---|
| 613 | } |
---|
| 614 | if (n > 2) { |
---|
| 615 | m = n >> 1; |
---|
| 616 | mh = m >> 1; |
---|
| 617 | for (j = 1; j < mh; j++) { |
---|
| 618 | k = m - j; |
---|
| 619 | xr = a[j] + a[n - j]; |
---|
| 620 | xi = a[j] - a[n - j]; |
---|
| 621 | yr = a[k] + a[n - k]; |
---|
| 622 | yi = a[k] - a[n - k]; |
---|
| 623 | a[j] = xr; |
---|
| 624 | a[k] = yr; |
---|
| 625 | t[j] = xi + yi; |
---|
| 626 | t[k] = xi - yi; |
---|
| 627 | } |
---|
| 628 | t[0] = a[mh] - a[n - mh]; |
---|
| 629 | a[mh] += a[n - mh]; |
---|
| 630 | a[0] = a[m]; |
---|
| 631 | dstsub(m, a, nc, w + nw); |
---|
| 632 | if (m > 4) { |
---|
| 633 | bitrv2(m, ip + 2, a); |
---|
| 634 | cftfsub(m, a, w); |
---|
| 635 | rftfsub(m, a, nc, w + nw); |
---|
| 636 | } else if (m == 4) { |
---|
| 637 | cftfsub(m, a, w); |
---|
| 638 | } |
---|
| 639 | a[n - 1] = a[1] - a[0]; |
---|
| 640 | a[1] = a[0] + a[1]; |
---|
| 641 | for (j = m - 2; j >= 2; j -= 2) { |
---|
| 642 | a[2 * j + 1] = a[j] - a[j + 1]; |
---|
| 643 | a[2 * j - 1] = -a[j] - a[j + 1]; |
---|
| 644 | } |
---|
| 645 | l = 2; |
---|
| 646 | m = mh; |
---|
| 647 | while (m >= 2) { |
---|
| 648 | dstsub(m, t, nc, w + nw); |
---|
| 649 | if (m > 4) { |
---|
| 650 | bitrv2(m, ip + 2, t); |
---|
| 651 | cftfsub(m, t, w); |
---|
| 652 | rftfsub(m, t, nc, w + nw); |
---|
| 653 | } else if (m == 4) { |
---|
| 654 | cftfsub(m, t, w); |
---|
| 655 | } |
---|
| 656 | a[n - l] = t[1] - t[0]; |
---|
| 657 | a[l] = t[0] + t[1]; |
---|
| 658 | k = 0; |
---|
| 659 | for (j = 2; j < m; j += 2) { |
---|
| 660 | k += l << 2; |
---|
| 661 | a[k - l] = -t[j] - t[j + 1]; |
---|
| 662 | a[k + l] = t[j] - t[j + 1]; |
---|
| 663 | } |
---|
| 664 | l <<= 1; |
---|
| 665 | mh = m >> 1; |
---|
| 666 | for (j = 1; j < mh; j++) { |
---|
| 667 | k = m - j; |
---|
| 668 | t[j] = t[m + k] + t[m + j]; |
---|
| 669 | t[k] = t[m + k] - t[m + j]; |
---|
| 670 | } |
---|
| 671 | t[0] = t[m + mh]; |
---|
| 672 | m = mh; |
---|
| 673 | } |
---|
| 674 | a[l] = t[0]; |
---|
| 675 | } |
---|
| 676 | a[0] = 0; |
---|
| 677 | } |
---|
| 678 | |
---|
| 679 | |
---|
| 680 | /* -------- initializing routines -------- */ |
---|
| 681 | |
---|
| 682 | |
---|
| 683 | #include <math.h> |
---|
| 684 | |
---|
[5c6b264] | 685 | void makewt(int nw, int *ip, smpl_t *w) |
---|
[cfaa3c4] | 686 | { |
---|
[5c6b264] | 687 | void bitrv2(int n, int *ip, smpl_t *a); |
---|
[cfaa3c4] | 688 | int j, nwh; |
---|
[5c6b264] | 689 | smpl_t delta, x, y; |
---|
[cfaa3c4] | 690 | |
---|
| 691 | ip[0] = nw; |
---|
| 692 | ip[1] = 1; |
---|
| 693 | if (nw > 2) { |
---|
| 694 | nwh = nw >> 1; |
---|
| 695 | delta = atan(1.0) / nwh; |
---|
| 696 | w[0] = 1; |
---|
| 697 | w[1] = 0; |
---|
[5cb8abe] | 698 | w[nwh] = COS(delta * nwh); |
---|
[cfaa3c4] | 699 | w[nwh + 1] = w[nwh]; |
---|
| 700 | if (nwh > 2) { |
---|
| 701 | for (j = 2; j < nwh; j += 2) { |
---|
[5cb8abe] | 702 | x = COS(delta * j); |
---|
| 703 | y = SIN(delta * j); |
---|
[cfaa3c4] | 704 | w[j] = x; |
---|
| 705 | w[j + 1] = y; |
---|
| 706 | w[nw - j] = y; |
---|
| 707 | w[nw - j + 1] = x; |
---|
| 708 | } |
---|
| 709 | for (j = nwh - 2; j >= 2; j -= 2) { |
---|
| 710 | x = w[2 * j]; |
---|
| 711 | y = w[2 * j + 1]; |
---|
| 712 | w[nwh + j] = x; |
---|
| 713 | w[nwh + j + 1] = y; |
---|
| 714 | } |
---|
| 715 | bitrv2(nw, ip + 2, w); |
---|
| 716 | } |
---|
| 717 | } |
---|
| 718 | } |
---|
| 719 | |
---|
| 720 | |
---|
[5c6b264] | 721 | void makect(int nc, int *ip, smpl_t *c) |
---|
[cfaa3c4] | 722 | { |
---|
| 723 | int j, nch; |
---|
[5c6b264] | 724 | smpl_t delta; |
---|
[cfaa3c4] | 725 | |
---|
| 726 | ip[1] = nc; |
---|
| 727 | if (nc > 1) { |
---|
| 728 | nch = nc >> 1; |
---|
| 729 | delta = atan(1.0) / nch; |
---|
| 730 | c[0] = cos(delta * nch); |
---|
| 731 | c[nch] = 0.5 * c[0]; |
---|
| 732 | for (j = 1; j < nch; j++) { |
---|
| 733 | c[j] = 0.5 * cos(delta * j); |
---|
| 734 | c[nc - j] = 0.5 * sin(delta * j); |
---|
| 735 | } |
---|
| 736 | } |
---|
| 737 | } |
---|
| 738 | |
---|
| 739 | |
---|
| 740 | /* -------- child routines -------- */ |
---|
| 741 | |
---|
| 742 | |
---|
[5c6b264] | 743 | void bitrv2(int n, int *ip, smpl_t *a) |
---|
[cfaa3c4] | 744 | { |
---|
| 745 | int j, j1, k, k1, l, m, m2; |
---|
[5c6b264] | 746 | smpl_t xr, xi, yr, yi; |
---|
[cfaa3c4] | 747 | |
---|
| 748 | ip[0] = 0; |
---|
| 749 | l = n; |
---|
| 750 | m = 1; |
---|
| 751 | while ((m << 3) < l) { |
---|
| 752 | l >>= 1; |
---|
| 753 | for (j = 0; j < m; j++) { |
---|
| 754 | ip[m + j] = ip[j] + l; |
---|
| 755 | } |
---|
| 756 | m <<= 1; |
---|
| 757 | } |
---|
| 758 | m2 = 2 * m; |
---|
| 759 | if ((m << 3) == l) { |
---|
| 760 | for (k = 0; k < m; k++) { |
---|
| 761 | for (j = 0; j < k; j++) { |
---|
| 762 | j1 = 2 * j + ip[k]; |
---|
| 763 | k1 = 2 * k + ip[j]; |
---|
| 764 | xr = a[j1]; |
---|
| 765 | xi = a[j1 + 1]; |
---|
| 766 | yr = a[k1]; |
---|
| 767 | yi = a[k1 + 1]; |
---|
| 768 | a[j1] = yr; |
---|
| 769 | a[j1 + 1] = yi; |
---|
| 770 | a[k1] = xr; |
---|
| 771 | a[k1 + 1] = xi; |
---|
| 772 | j1 += m2; |
---|
| 773 | k1 += 2 * m2; |
---|
| 774 | xr = a[j1]; |
---|
| 775 | xi = a[j1 + 1]; |
---|
| 776 | yr = a[k1]; |
---|
| 777 | yi = a[k1 + 1]; |
---|
| 778 | a[j1] = yr; |
---|
| 779 | a[j1 + 1] = yi; |
---|
| 780 | a[k1] = xr; |
---|
| 781 | a[k1 + 1] = xi; |
---|
| 782 | j1 += m2; |
---|
| 783 | k1 -= m2; |
---|
| 784 | xr = a[j1]; |
---|
| 785 | xi = a[j1 + 1]; |
---|
| 786 | yr = a[k1]; |
---|
| 787 | yi = a[k1 + 1]; |
---|
| 788 | a[j1] = yr; |
---|
| 789 | a[j1 + 1] = yi; |
---|
| 790 | a[k1] = xr; |
---|
| 791 | a[k1 + 1] = xi; |
---|
| 792 | j1 += m2; |
---|
| 793 | k1 += 2 * m2; |
---|
| 794 | xr = a[j1]; |
---|
| 795 | xi = a[j1 + 1]; |
---|
| 796 | yr = a[k1]; |
---|
| 797 | yi = a[k1 + 1]; |
---|
| 798 | a[j1] = yr; |
---|
| 799 | a[j1 + 1] = yi; |
---|
| 800 | a[k1] = xr; |
---|
| 801 | a[k1 + 1] = xi; |
---|
| 802 | } |
---|
| 803 | j1 = 2 * k + m2 + ip[k]; |
---|
| 804 | k1 = j1 + m2; |
---|
| 805 | xr = a[j1]; |
---|
| 806 | xi = a[j1 + 1]; |
---|
| 807 | yr = a[k1]; |
---|
| 808 | yi = a[k1 + 1]; |
---|
| 809 | a[j1] = yr; |
---|
| 810 | a[j1 + 1] = yi; |
---|
| 811 | a[k1] = xr; |
---|
| 812 | a[k1 + 1] = xi; |
---|
| 813 | } |
---|
| 814 | } else { |
---|
| 815 | for (k = 1; k < m; k++) { |
---|
| 816 | for (j = 0; j < k; j++) { |
---|
| 817 | j1 = 2 * j + ip[k]; |
---|
| 818 | k1 = 2 * k + ip[j]; |
---|
| 819 | xr = a[j1]; |
---|
| 820 | xi = a[j1 + 1]; |
---|
| 821 | yr = a[k1]; |
---|
| 822 | yi = a[k1 + 1]; |
---|
| 823 | a[j1] = yr; |
---|
| 824 | a[j1 + 1] = yi; |
---|
| 825 | a[k1] = xr; |
---|
| 826 | a[k1 + 1] = xi; |
---|
| 827 | j1 += m2; |
---|
| 828 | k1 += m2; |
---|
| 829 | xr = a[j1]; |
---|
| 830 | xi = a[j1 + 1]; |
---|
| 831 | yr = a[k1]; |
---|
| 832 | yi = a[k1 + 1]; |
---|
| 833 | a[j1] = yr; |
---|
| 834 | a[j1 + 1] = yi; |
---|
| 835 | a[k1] = xr; |
---|
| 836 | a[k1 + 1] = xi; |
---|
| 837 | } |
---|
| 838 | } |
---|
| 839 | } |
---|
| 840 | } |
---|
| 841 | |
---|
| 842 | |
---|
[5c6b264] | 843 | void bitrv2conj(int n, int *ip, smpl_t *a) |
---|
[cfaa3c4] | 844 | { |
---|
| 845 | int j, j1, k, k1, l, m, m2; |
---|
[5c6b264] | 846 | smpl_t xr, xi, yr, yi; |
---|
[cfaa3c4] | 847 | |
---|
| 848 | ip[0] = 0; |
---|
| 849 | l = n; |
---|
| 850 | m = 1; |
---|
| 851 | while ((m << 3) < l) { |
---|
| 852 | l >>= 1; |
---|
| 853 | for (j = 0; j < m; j++) { |
---|
| 854 | ip[m + j] = ip[j] + l; |
---|
| 855 | } |
---|
| 856 | m <<= 1; |
---|
| 857 | } |
---|
| 858 | m2 = 2 * m; |
---|
| 859 | if ((m << 3) == l) { |
---|
| 860 | for (k = 0; k < m; k++) { |
---|
| 861 | for (j = 0; j < k; j++) { |
---|
| 862 | j1 = 2 * j + ip[k]; |
---|
| 863 | k1 = 2 * k + ip[j]; |
---|
| 864 | xr = a[j1]; |
---|
| 865 | xi = -a[j1 + 1]; |
---|
| 866 | yr = a[k1]; |
---|
| 867 | yi = -a[k1 + 1]; |
---|
| 868 | a[j1] = yr; |
---|
| 869 | a[j1 + 1] = yi; |
---|
| 870 | a[k1] = xr; |
---|
| 871 | a[k1 + 1] = xi; |
---|
| 872 | j1 += m2; |
---|
| 873 | k1 += 2 * m2; |
---|
| 874 | xr = a[j1]; |
---|
| 875 | xi = -a[j1 + 1]; |
---|
| 876 | yr = a[k1]; |
---|
| 877 | yi = -a[k1 + 1]; |
---|
| 878 | a[j1] = yr; |
---|
| 879 | a[j1 + 1] = yi; |
---|
| 880 | a[k1] = xr; |
---|
| 881 | a[k1 + 1] = xi; |
---|
| 882 | j1 += m2; |
---|
| 883 | k1 -= m2; |
---|
| 884 | xr = a[j1]; |
---|
| 885 | xi = -a[j1 + 1]; |
---|
| 886 | yr = a[k1]; |
---|
| 887 | yi = -a[k1 + 1]; |
---|
| 888 | a[j1] = yr; |
---|
| 889 | a[j1 + 1] = yi; |
---|
| 890 | a[k1] = xr; |
---|
| 891 | a[k1 + 1] = xi; |
---|
| 892 | j1 += m2; |
---|
| 893 | k1 += 2 * m2; |
---|
| 894 | xr = a[j1]; |
---|
| 895 | xi = -a[j1 + 1]; |
---|
| 896 | yr = a[k1]; |
---|
| 897 | yi = -a[k1 + 1]; |
---|
| 898 | a[j1] = yr; |
---|
| 899 | a[j1 + 1] = yi; |
---|
| 900 | a[k1] = xr; |
---|
| 901 | a[k1 + 1] = xi; |
---|
| 902 | } |
---|
| 903 | k1 = 2 * k + ip[k]; |
---|
| 904 | a[k1 + 1] = -a[k1 + 1]; |
---|
| 905 | j1 = k1 + m2; |
---|
| 906 | k1 = j1 + m2; |
---|
| 907 | xr = a[j1]; |
---|
| 908 | xi = -a[j1 + 1]; |
---|
| 909 | yr = a[k1]; |
---|
| 910 | yi = -a[k1 + 1]; |
---|
| 911 | a[j1] = yr; |
---|
| 912 | a[j1 + 1] = yi; |
---|
| 913 | a[k1] = xr; |
---|
| 914 | a[k1 + 1] = xi; |
---|
| 915 | k1 += m2; |
---|
| 916 | a[k1 + 1] = -a[k1 + 1]; |
---|
| 917 | } |
---|
| 918 | } else { |
---|
| 919 | a[1] = -a[1]; |
---|
| 920 | a[m2 + 1] = -a[m2 + 1]; |
---|
| 921 | for (k = 1; k < m; k++) { |
---|
| 922 | for (j = 0; j < k; j++) { |
---|
| 923 | j1 = 2 * j + ip[k]; |
---|
| 924 | k1 = 2 * k + ip[j]; |
---|
| 925 | xr = a[j1]; |
---|
| 926 | xi = -a[j1 + 1]; |
---|
| 927 | yr = a[k1]; |
---|
| 928 | yi = -a[k1 + 1]; |
---|
| 929 | a[j1] = yr; |
---|
| 930 | a[j1 + 1] = yi; |
---|
| 931 | a[k1] = xr; |
---|
| 932 | a[k1 + 1] = xi; |
---|
| 933 | j1 += m2; |
---|
| 934 | k1 += m2; |
---|
| 935 | xr = a[j1]; |
---|
| 936 | xi = -a[j1 + 1]; |
---|
| 937 | yr = a[k1]; |
---|
| 938 | yi = -a[k1 + 1]; |
---|
| 939 | a[j1] = yr; |
---|
| 940 | a[j1 + 1] = yi; |
---|
| 941 | a[k1] = xr; |
---|
| 942 | a[k1 + 1] = xi; |
---|
| 943 | } |
---|
| 944 | k1 = 2 * k + ip[k]; |
---|
| 945 | a[k1 + 1] = -a[k1 + 1]; |
---|
| 946 | a[k1 + m2 + 1] = -a[k1 + m2 + 1]; |
---|
| 947 | } |
---|
| 948 | } |
---|
| 949 | } |
---|
| 950 | |
---|
| 951 | |
---|
[5c6b264] | 952 | void cftfsub(int n, smpl_t *a, smpl_t *w) |
---|
[cfaa3c4] | 953 | { |
---|
[5c6b264] | 954 | void cft1st(int n, smpl_t *a, smpl_t *w); |
---|
| 955 | void cftmdl(int n, int l, smpl_t *a, smpl_t *w); |
---|
[cfaa3c4] | 956 | int j, j1, j2, j3, l; |
---|
[5c6b264] | 957 | smpl_t x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; |
---|
[cfaa3c4] | 958 | |
---|
| 959 | l = 2; |
---|
| 960 | if (n >= 16) { |
---|
| 961 | cft1st(n, a, w); |
---|
| 962 | l = 16; |
---|
| 963 | while ((l << 3) <= n) { |
---|
| 964 | cftmdl(n, l, a, w); |
---|
| 965 | l <<= 3; |
---|
| 966 | } |
---|
| 967 | } |
---|
| 968 | if ((l << 1) < n) { |
---|
| 969 | for (j = 0; j < l; j += 2) { |
---|
| 970 | j1 = j + l; |
---|
| 971 | j2 = j1 + l; |
---|
| 972 | j3 = j2 + l; |
---|
| 973 | x0r = a[j] + a[j1]; |
---|
| 974 | x0i = a[j + 1] + a[j1 + 1]; |
---|
| 975 | x1r = a[j] - a[j1]; |
---|
| 976 | x1i = a[j + 1] - a[j1 + 1]; |
---|
| 977 | x2r = a[j2] + a[j3]; |
---|
| 978 | x2i = a[j2 + 1] + a[j3 + 1]; |
---|
| 979 | x3r = a[j2] - a[j3]; |
---|
| 980 | x3i = a[j2 + 1] - a[j3 + 1]; |
---|
| 981 | a[j] = x0r + x2r; |
---|
| 982 | a[j + 1] = x0i + x2i; |
---|
| 983 | a[j2] = x0r - x2r; |
---|
| 984 | a[j2 + 1] = x0i - x2i; |
---|
| 985 | a[j1] = x1r - x3i; |
---|
| 986 | a[j1 + 1] = x1i + x3r; |
---|
| 987 | a[j3] = x1r + x3i; |
---|
| 988 | a[j3 + 1] = x1i - x3r; |
---|
| 989 | } |
---|
| 990 | } else if ((l << 1) == n) { |
---|
| 991 | for (j = 0; j < l; j += 2) { |
---|
| 992 | j1 = j + l; |
---|
| 993 | x0r = a[j] - a[j1]; |
---|
| 994 | x0i = a[j + 1] - a[j1 + 1]; |
---|
| 995 | a[j] += a[j1]; |
---|
| 996 | a[j + 1] += a[j1 + 1]; |
---|
| 997 | a[j1] = x0r; |
---|
| 998 | a[j1 + 1] = x0i; |
---|
| 999 | } |
---|
| 1000 | } |
---|
| 1001 | } |
---|
| 1002 | |
---|
| 1003 | |
---|
[5c6b264] | 1004 | void cftbsub(int n, smpl_t *a, smpl_t *w) |
---|
[cfaa3c4] | 1005 | { |
---|
[5c6b264] | 1006 | void cft1st(int n, smpl_t *a, smpl_t *w); |
---|
| 1007 | void cftmdl(int n, int l, smpl_t *a, smpl_t *w); |
---|
[cfaa3c4] | 1008 | int j, j1, j2, j3, j4, j5, j6, j7, l; |
---|
[5c6b264] | 1009 | smpl_t wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, |
---|
[cfaa3c4] | 1010 | y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, |
---|
| 1011 | y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i; |
---|
| 1012 | |
---|
| 1013 | l = 2; |
---|
| 1014 | if (n > 16) { |
---|
| 1015 | cft1st(n, a, w); |
---|
| 1016 | l = 16; |
---|
| 1017 | while ((l << 3) < n) { |
---|
| 1018 | cftmdl(n, l, a, w); |
---|
| 1019 | l <<= 3; |
---|
| 1020 | } |
---|
| 1021 | } |
---|
| 1022 | if ((l << 2) < n) { |
---|
| 1023 | wn4r = w[2]; |
---|
| 1024 | for (j = 0; j < l; j += 2) { |
---|
| 1025 | j1 = j + l; |
---|
| 1026 | j2 = j1 + l; |
---|
| 1027 | j3 = j2 + l; |
---|
| 1028 | j4 = j3 + l; |
---|
| 1029 | j5 = j4 + l; |
---|
| 1030 | j6 = j5 + l; |
---|
| 1031 | j7 = j6 + l; |
---|
| 1032 | x0r = a[j] + a[j1]; |
---|
| 1033 | x0i = -a[j + 1] - a[j1 + 1]; |
---|
| 1034 | x1r = a[j] - a[j1]; |
---|
| 1035 | x1i = -a[j + 1] + a[j1 + 1]; |
---|
| 1036 | x2r = a[j2] + a[j3]; |
---|
| 1037 | x2i = a[j2 + 1] + a[j3 + 1]; |
---|
| 1038 | x3r = a[j2] - a[j3]; |
---|
| 1039 | x3i = a[j2 + 1] - a[j3 + 1]; |
---|
| 1040 | y0r = x0r + x2r; |
---|
| 1041 | y0i = x0i - x2i; |
---|
| 1042 | y2r = x0r - x2r; |
---|
| 1043 | y2i = x0i + x2i; |
---|
| 1044 | y1r = x1r - x3i; |
---|
| 1045 | y1i = x1i - x3r; |
---|
| 1046 | y3r = x1r + x3i; |
---|
| 1047 | y3i = x1i + x3r; |
---|
| 1048 | x0r = a[j4] + a[j5]; |
---|
| 1049 | x0i = a[j4 + 1] + a[j5 + 1]; |
---|
| 1050 | x1r = a[j4] - a[j5]; |
---|
| 1051 | x1i = a[j4 + 1] - a[j5 + 1]; |
---|
| 1052 | x2r = a[j6] + a[j7]; |
---|
| 1053 | x2i = a[j6 + 1] + a[j7 + 1]; |
---|
| 1054 | x3r = a[j6] - a[j7]; |
---|
| 1055 | x3i = a[j6 + 1] - a[j7 + 1]; |
---|
| 1056 | y4r = x0r + x2r; |
---|
| 1057 | y4i = x0i + x2i; |
---|
| 1058 | y6r = x0r - x2r; |
---|
| 1059 | y6i = x0i - x2i; |
---|
| 1060 | x0r = x1r - x3i; |
---|
| 1061 | x0i = x1i + x3r; |
---|
| 1062 | x2r = x1r + x3i; |
---|
| 1063 | x2i = x1i - x3r; |
---|
| 1064 | y5r = wn4r * (x0r - x0i); |
---|
| 1065 | y5i = wn4r * (x0r + x0i); |
---|
| 1066 | y7r = wn4r * (x2r - x2i); |
---|
| 1067 | y7i = wn4r * (x2r + x2i); |
---|
| 1068 | a[j1] = y1r + y5r; |
---|
| 1069 | a[j1 + 1] = y1i - y5i; |
---|
| 1070 | a[j5] = y1r - y5r; |
---|
| 1071 | a[j5 + 1] = y1i + y5i; |
---|
| 1072 | a[j3] = y3r - y7i; |
---|
| 1073 | a[j3 + 1] = y3i - y7r; |
---|
| 1074 | a[j7] = y3r + y7i; |
---|
| 1075 | a[j7 + 1] = y3i + y7r; |
---|
| 1076 | a[j] = y0r + y4r; |
---|
| 1077 | a[j + 1] = y0i - y4i; |
---|
| 1078 | a[j4] = y0r - y4r; |
---|
| 1079 | a[j4 + 1] = y0i + y4i; |
---|
| 1080 | a[j2] = y2r - y6i; |
---|
| 1081 | a[j2 + 1] = y2i - y6r; |
---|
| 1082 | a[j6] = y2r + y6i; |
---|
| 1083 | a[j6 + 1] = y2i + y6r; |
---|
| 1084 | } |
---|
| 1085 | } else if ((l << 2) == n) { |
---|
| 1086 | for (j = 0; j < l; j += 2) { |
---|
| 1087 | j1 = j + l; |
---|
| 1088 | j2 = j1 + l; |
---|
| 1089 | j3 = j2 + l; |
---|
| 1090 | x0r = a[j] + a[j1]; |
---|
| 1091 | x0i = -a[j + 1] - a[j1 + 1]; |
---|
| 1092 | x1r = a[j] - a[j1]; |
---|
| 1093 | x1i = -a[j + 1] + a[j1 + 1]; |
---|
| 1094 | x2r = a[j2] + a[j3]; |
---|
| 1095 | x2i = a[j2 + 1] + a[j3 + 1]; |
---|
| 1096 | x3r = a[j2] - a[j3]; |
---|
| 1097 | x3i = a[j2 + 1] - a[j3 + 1]; |
---|
| 1098 | a[j] = x0r + x2r; |
---|
| 1099 | a[j + 1] = x0i - x2i; |
---|
| 1100 | a[j2] = x0r - x2r; |
---|
| 1101 | a[j2 + 1] = x0i + x2i; |
---|
| 1102 | a[j1] = x1r - x3i; |
---|
| 1103 | a[j1 + 1] = x1i - x3r; |
---|
| 1104 | a[j3] = x1r + x3i; |
---|
| 1105 | a[j3 + 1] = x1i + x3r; |
---|
| 1106 | } |
---|
| 1107 | } else { |
---|
| 1108 | for (j = 0; j < l; j += 2) { |
---|
| 1109 | j1 = j + l; |
---|
| 1110 | x0r = a[j] - a[j1]; |
---|
| 1111 | x0i = -a[j + 1] + a[j1 + 1]; |
---|
| 1112 | a[j] += a[j1]; |
---|
| 1113 | a[j + 1] = -a[j + 1] - a[j1 + 1]; |
---|
| 1114 | a[j1] = x0r; |
---|
| 1115 | a[j1 + 1] = x0i; |
---|
| 1116 | } |
---|
| 1117 | } |
---|
| 1118 | } |
---|
| 1119 | |
---|
| 1120 | |
---|
[5c6b264] | 1121 | void cft1st(int n, smpl_t *a, smpl_t *w) |
---|
[cfaa3c4] | 1122 | { |
---|
| 1123 | int j, k1; |
---|
[5c6b264] | 1124 | smpl_t wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, |
---|
[cfaa3c4] | 1125 | wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i; |
---|
[5c6b264] | 1126 | smpl_t x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, |
---|
[cfaa3c4] | 1127 | y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, |
---|
| 1128 | y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i; |
---|
| 1129 | |
---|
| 1130 | wn4r = w[2]; |
---|
| 1131 | x0r = a[0] + a[2]; |
---|
| 1132 | x0i = a[1] + a[3]; |
---|
| 1133 | x1r = a[0] - a[2]; |
---|
| 1134 | x1i = a[1] - a[3]; |
---|
| 1135 | x2r = a[4] + a[6]; |
---|
| 1136 | x2i = a[5] + a[7]; |
---|
| 1137 | x3r = a[4] - a[6]; |
---|
| 1138 | x3i = a[5] - a[7]; |
---|
| 1139 | y0r = x0r + x2r; |
---|
| 1140 | y0i = x0i + x2i; |
---|
| 1141 | y2r = x0r - x2r; |
---|
| 1142 | y2i = x0i - x2i; |
---|
| 1143 | y1r = x1r - x3i; |
---|
| 1144 | y1i = x1i + x3r; |
---|
| 1145 | y3r = x1r + x3i; |
---|
| 1146 | y3i = x1i - x3r; |
---|
| 1147 | x0r = a[8] + a[10]; |
---|
| 1148 | x0i = a[9] + a[11]; |
---|
| 1149 | x1r = a[8] - a[10]; |
---|
| 1150 | x1i = a[9] - a[11]; |
---|
| 1151 | x2r = a[12] + a[14]; |
---|
| 1152 | x2i = a[13] + a[15]; |
---|
| 1153 | x3r = a[12] - a[14]; |
---|
| 1154 | x3i = a[13] - a[15]; |
---|
| 1155 | y4r = x0r + x2r; |
---|
| 1156 | y4i = x0i + x2i; |
---|
| 1157 | y6r = x0r - x2r; |
---|
| 1158 | y6i = x0i - x2i; |
---|
| 1159 | x0r = x1r - x3i; |
---|
| 1160 | x0i = x1i + x3r; |
---|
| 1161 | x2r = x1r + x3i; |
---|
| 1162 | x2i = x1i - x3r; |
---|
| 1163 | y5r = wn4r * (x0r - x0i); |
---|
| 1164 | y5i = wn4r * (x0r + x0i); |
---|
| 1165 | y7r = wn4r * (x2r - x2i); |
---|
| 1166 | y7i = wn4r * (x2r + x2i); |
---|
| 1167 | a[2] = y1r + y5r; |
---|
| 1168 | a[3] = y1i + y5i; |
---|
| 1169 | a[10] = y1r - y5r; |
---|
| 1170 | a[11] = y1i - y5i; |
---|
| 1171 | a[6] = y3r - y7i; |
---|
| 1172 | a[7] = y3i + y7r; |
---|
| 1173 | a[14] = y3r + y7i; |
---|
| 1174 | a[15] = y3i - y7r; |
---|
| 1175 | a[0] = y0r + y4r; |
---|
| 1176 | a[1] = y0i + y4i; |
---|
| 1177 | a[8] = y0r - y4r; |
---|
| 1178 | a[9] = y0i - y4i; |
---|
| 1179 | a[4] = y2r - y6i; |
---|
| 1180 | a[5] = y2i + y6r; |
---|
| 1181 | a[12] = y2r + y6i; |
---|
| 1182 | a[13] = y2i - y6r; |
---|
| 1183 | if (n > 16) { |
---|
| 1184 | wk1r = w[4]; |
---|
| 1185 | wk1i = w[5]; |
---|
| 1186 | x0r = a[16] + a[18]; |
---|
| 1187 | x0i = a[17] + a[19]; |
---|
| 1188 | x1r = a[16] - a[18]; |
---|
| 1189 | x1i = a[17] - a[19]; |
---|
| 1190 | x2r = a[20] + a[22]; |
---|
| 1191 | x2i = a[21] + a[23]; |
---|
| 1192 | x3r = a[20] - a[22]; |
---|
| 1193 | x3i = a[21] - a[23]; |
---|
| 1194 | y0r = x0r + x2r; |
---|
| 1195 | y0i = x0i + x2i; |
---|
| 1196 | y2r = x0r - x2r; |
---|
| 1197 | y2i = x0i - x2i; |
---|
| 1198 | y1r = x1r - x3i; |
---|
| 1199 | y1i = x1i + x3r; |
---|
| 1200 | y3r = x1r + x3i; |
---|
| 1201 | y3i = x1i - x3r; |
---|
| 1202 | x0r = a[24] + a[26]; |
---|
| 1203 | x0i = a[25] + a[27]; |
---|
| 1204 | x1r = a[24] - a[26]; |
---|
| 1205 | x1i = a[25] - a[27]; |
---|
| 1206 | x2r = a[28] + a[30]; |
---|
| 1207 | x2i = a[29] + a[31]; |
---|
| 1208 | x3r = a[28] - a[30]; |
---|
| 1209 | x3i = a[29] - a[31]; |
---|
| 1210 | y4r = x0r + x2r; |
---|
| 1211 | y4i = x0i + x2i; |
---|
| 1212 | y6r = x0r - x2r; |
---|
| 1213 | y6i = x0i - x2i; |
---|
| 1214 | x0r = x1r - x3i; |
---|
| 1215 | x0i = x1i + x3r; |
---|
| 1216 | x2r = x1r + x3i; |
---|
| 1217 | x2i = x3r - x1i; |
---|
| 1218 | y5r = wk1i * x0r - wk1r * x0i; |
---|
| 1219 | y5i = wk1i * x0i + wk1r * x0r; |
---|
| 1220 | y7r = wk1r * x2r + wk1i * x2i; |
---|
| 1221 | y7i = wk1r * x2i - wk1i * x2r; |
---|
| 1222 | x0r = wk1r * y1r - wk1i * y1i; |
---|
| 1223 | x0i = wk1r * y1i + wk1i * y1r; |
---|
| 1224 | a[18] = x0r + y5r; |
---|
| 1225 | a[19] = x0i + y5i; |
---|
| 1226 | a[26] = y5i - x0i; |
---|
| 1227 | a[27] = x0r - y5r; |
---|
| 1228 | x0r = wk1i * y3r - wk1r * y3i; |
---|
| 1229 | x0i = wk1i * y3i + wk1r * y3r; |
---|
| 1230 | a[22] = x0r - y7r; |
---|
| 1231 | a[23] = x0i + y7i; |
---|
| 1232 | a[30] = y7i - x0i; |
---|
| 1233 | a[31] = x0r + y7r; |
---|
| 1234 | a[16] = y0r + y4r; |
---|
| 1235 | a[17] = y0i + y4i; |
---|
| 1236 | a[24] = y4i - y0i; |
---|
| 1237 | a[25] = y0r - y4r; |
---|
| 1238 | x0r = y2r - y6i; |
---|
| 1239 | x0i = y2i + y6r; |
---|
| 1240 | a[20] = wn4r * (x0r - x0i); |
---|
| 1241 | a[21] = wn4r * (x0i + x0r); |
---|
| 1242 | x0r = y6r - y2i; |
---|
| 1243 | x0i = y2r + y6i; |
---|
| 1244 | a[28] = wn4r * (x0r - x0i); |
---|
| 1245 | a[29] = wn4r * (x0i + x0r); |
---|
| 1246 | k1 = 4; |
---|
| 1247 | for (j = 32; j < n; j += 16) { |
---|
| 1248 | k1 += 4; |
---|
| 1249 | wk1r = w[k1]; |
---|
| 1250 | wk1i = w[k1 + 1]; |
---|
| 1251 | wk2r = w[k1 + 2]; |
---|
| 1252 | wk2i = w[k1 + 3]; |
---|
| 1253 | wtmp = 2 * wk2i; |
---|
| 1254 | wk3r = wk1r - wtmp * wk1i; |
---|
| 1255 | wk3i = wtmp * wk1r - wk1i; |
---|
| 1256 | wk4r = 1 - wtmp * wk2i; |
---|
| 1257 | wk4i = wtmp * wk2r; |
---|
| 1258 | wtmp = 2 * wk4i; |
---|
| 1259 | wk5r = wk3r - wtmp * wk1i; |
---|
| 1260 | wk5i = wtmp * wk1r - wk3i; |
---|
| 1261 | wk6r = wk2r - wtmp * wk2i; |
---|
| 1262 | wk6i = wtmp * wk2r - wk2i; |
---|
| 1263 | wk7r = wk1r - wtmp * wk3i; |
---|
| 1264 | wk7i = wtmp * wk3r - wk1i; |
---|
| 1265 | x0r = a[j] + a[j + 2]; |
---|
| 1266 | x0i = a[j + 1] + a[j + 3]; |
---|
| 1267 | x1r = a[j] - a[j + 2]; |
---|
| 1268 | x1i = a[j + 1] - a[j + 3]; |
---|
| 1269 | x2r = a[j + 4] + a[j + 6]; |
---|
| 1270 | x2i = a[j + 5] + a[j + 7]; |
---|
| 1271 | x3r = a[j + 4] - a[j + 6]; |
---|
| 1272 | x3i = a[j + 5] - a[j + 7]; |
---|
| 1273 | y0r = x0r + x2r; |
---|
| 1274 | y0i = x0i + x2i; |
---|
| 1275 | y2r = x0r - x2r; |
---|
| 1276 | y2i = x0i - x2i; |
---|
| 1277 | y1r = x1r - x3i; |
---|
| 1278 | y1i = x1i + x3r; |
---|
| 1279 | y3r = x1r + x3i; |
---|
| 1280 | y3i = x1i - x3r; |
---|
| 1281 | x0r = a[j + 8] + a[j + 10]; |
---|
| 1282 | x0i = a[j + 9] + a[j + 11]; |
---|
| 1283 | x1r = a[j + 8] - a[j + 10]; |
---|
| 1284 | x1i = a[j + 9] - a[j + 11]; |
---|
| 1285 | x2r = a[j + 12] + a[j + 14]; |
---|
| 1286 | x2i = a[j + 13] + a[j + 15]; |
---|
| 1287 | x3r = a[j + 12] - a[j + 14]; |
---|
| 1288 | x3i = a[j + 13] - a[j + 15]; |
---|
| 1289 | y4r = x0r + x2r; |
---|
| 1290 | y4i = x0i + x2i; |
---|
| 1291 | y6r = x0r - x2r; |
---|
| 1292 | y6i = x0i - x2i; |
---|
| 1293 | x0r = x1r - x3i; |
---|
| 1294 | x0i = x1i + x3r; |
---|
| 1295 | x2r = x1r + x3i; |
---|
| 1296 | x2i = x1i - x3r; |
---|
| 1297 | y5r = wn4r * (x0r - x0i); |
---|
| 1298 | y5i = wn4r * (x0r + x0i); |
---|
| 1299 | y7r = wn4r * (x2r - x2i); |
---|
| 1300 | y7i = wn4r * (x2r + x2i); |
---|
| 1301 | x0r = y1r + y5r; |
---|
| 1302 | x0i = y1i + y5i; |
---|
| 1303 | a[j + 2] = wk1r * x0r - wk1i * x0i; |
---|
| 1304 | a[j + 3] = wk1r * x0i + wk1i * x0r; |
---|
| 1305 | x0r = y1r - y5r; |
---|
| 1306 | x0i = y1i - y5i; |
---|
| 1307 | a[j + 10] = wk5r * x0r - wk5i * x0i; |
---|
| 1308 | a[j + 11] = wk5r * x0i + wk5i * x0r; |
---|
| 1309 | x0r = y3r - y7i; |
---|
| 1310 | x0i = y3i + y7r; |
---|
| 1311 | a[j + 6] = wk3r * x0r - wk3i * x0i; |
---|
| 1312 | a[j + 7] = wk3r * x0i + wk3i * x0r; |
---|
| 1313 | x0r = y3r + y7i; |
---|
| 1314 | x0i = y3i - y7r; |
---|
| 1315 | a[j + 14] = wk7r * x0r - wk7i * x0i; |
---|
| 1316 | a[j + 15] = wk7r * x0i + wk7i * x0r; |
---|
| 1317 | a[j] = y0r + y4r; |
---|
| 1318 | a[j + 1] = y0i + y4i; |
---|
| 1319 | x0r = y0r - y4r; |
---|
| 1320 | x0i = y0i - y4i; |
---|
| 1321 | a[j + 8] = wk4r * x0r - wk4i * x0i; |
---|
| 1322 | a[j + 9] = wk4r * x0i + wk4i * x0r; |
---|
| 1323 | x0r = y2r - y6i; |
---|
| 1324 | x0i = y2i + y6r; |
---|
| 1325 | a[j + 4] = wk2r * x0r - wk2i * x0i; |
---|
| 1326 | a[j + 5] = wk2r * x0i + wk2i * x0r; |
---|
| 1327 | x0r = y2r + y6i; |
---|
| 1328 | x0i = y2i - y6r; |
---|
| 1329 | a[j + 12] = wk6r * x0r - wk6i * x0i; |
---|
| 1330 | a[j + 13] = wk6r * x0i + wk6i * x0r; |
---|
| 1331 | } |
---|
| 1332 | } |
---|
| 1333 | } |
---|
| 1334 | |
---|
| 1335 | |
---|
[5c6b264] | 1336 | void cftmdl(int n, int l, smpl_t *a, smpl_t *w) |
---|
[cfaa3c4] | 1337 | { |
---|
| 1338 | int j, j1, j2, j3, j4, j5, j6, j7, k, k1, m; |
---|
[5c6b264] | 1339 | smpl_t wn4r, wtmp, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, |
---|
[cfaa3c4] | 1340 | wk4r, wk4i, wk5r, wk5i, wk6r, wk6i, wk7r, wk7i; |
---|
[5c6b264] | 1341 | smpl_t x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, |
---|
[cfaa3c4] | 1342 | y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, |
---|
| 1343 | y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i; |
---|
| 1344 | |
---|
| 1345 | m = l << 3; |
---|
| 1346 | wn4r = w[2]; |
---|
| 1347 | for (j = 0; j < l; j += 2) { |
---|
| 1348 | j1 = j + l; |
---|
| 1349 | j2 = j1 + l; |
---|
| 1350 | j3 = j2 + l; |
---|
| 1351 | j4 = j3 + l; |
---|
| 1352 | j5 = j4 + l; |
---|
| 1353 | j6 = j5 + l; |
---|
| 1354 | j7 = j6 + l; |
---|
| 1355 | x0r = a[j] + a[j1]; |
---|
| 1356 | x0i = a[j + 1] + a[j1 + 1]; |
---|
| 1357 | x1r = a[j] - a[j1]; |
---|
| 1358 | x1i = a[j + 1] - a[j1 + 1]; |
---|
| 1359 | x2r = a[j2] + a[j3]; |
---|
| 1360 | x2i = a[j2 + 1] + a[j3 + 1]; |
---|
| 1361 | x3r = a[j2] - a[j3]; |
---|
| 1362 | x3i = a[j2 + 1] - a[j3 + 1]; |
---|
| 1363 | y0r = x0r + x2r; |
---|
| 1364 | y0i = x0i + x2i; |
---|
| 1365 | y2r = x0r - x2r; |
---|
| 1366 | y2i = x0i - x2i; |
---|
| 1367 | y1r = x1r - x3i; |
---|
| 1368 | y1i = x1i + x3r; |
---|
| 1369 | y3r = x1r + x3i; |
---|
| 1370 | y3i = x1i - x3r; |
---|
| 1371 | x0r = a[j4] + a[j5]; |
---|
| 1372 | x0i = a[j4 + 1] + a[j5 + 1]; |
---|
| 1373 | x1r = a[j4] - a[j5]; |
---|
| 1374 | x1i = a[j4 + 1] - a[j5 + 1]; |
---|
| 1375 | x2r = a[j6] + a[j7]; |
---|
| 1376 | x2i = a[j6 + 1] + a[j7 + 1]; |
---|
| 1377 | x3r = a[j6] - a[j7]; |
---|
| 1378 | x3i = a[j6 + 1] - a[j7 + 1]; |
---|
| 1379 | y4r = x0r + x2r; |
---|
| 1380 | y4i = x0i + x2i; |
---|
| 1381 | y6r = x0r - x2r; |
---|
| 1382 | y6i = x0i - x2i; |
---|
| 1383 | x0r = x1r - x3i; |
---|
| 1384 | x0i = x1i + x3r; |
---|
| 1385 | x2r = x1r + x3i; |
---|
| 1386 | x2i = x1i - x3r; |
---|
| 1387 | y5r = wn4r * (x0r - x0i); |
---|
| 1388 | y5i = wn4r * (x0r + x0i); |
---|
| 1389 | y7r = wn4r * (x2r - x2i); |
---|
| 1390 | y7i = wn4r * (x2r + x2i); |
---|
| 1391 | a[j1] = y1r + y5r; |
---|
| 1392 | a[j1 + 1] = y1i + y5i; |
---|
| 1393 | a[j5] = y1r - y5r; |
---|
| 1394 | a[j5 + 1] = y1i - y5i; |
---|
| 1395 | a[j3] = y3r - y7i; |
---|
| 1396 | a[j3 + 1] = y3i + y7r; |
---|
| 1397 | a[j7] = y3r + y7i; |
---|
| 1398 | a[j7 + 1] = y3i - y7r; |
---|
| 1399 | a[j] = y0r + y4r; |
---|
| 1400 | a[j + 1] = y0i + y4i; |
---|
| 1401 | a[j4] = y0r - y4r; |
---|
| 1402 | a[j4 + 1] = y0i - y4i; |
---|
| 1403 | a[j2] = y2r - y6i; |
---|
| 1404 | a[j2 + 1] = y2i + y6r; |
---|
| 1405 | a[j6] = y2r + y6i; |
---|
| 1406 | a[j6 + 1] = y2i - y6r; |
---|
| 1407 | } |
---|
| 1408 | if (m < n) { |
---|
| 1409 | wk1r = w[4]; |
---|
| 1410 | wk1i = w[5]; |
---|
| 1411 | for (j = m; j < l + m; j += 2) { |
---|
| 1412 | j1 = j + l; |
---|
| 1413 | j2 = j1 + l; |
---|
| 1414 | j3 = j2 + l; |
---|
| 1415 | j4 = j3 + l; |
---|
| 1416 | j5 = j4 + l; |
---|
| 1417 | j6 = j5 + l; |
---|
| 1418 | j7 = j6 + l; |
---|
| 1419 | x0r = a[j] + a[j1]; |
---|
| 1420 | x0i = a[j + 1] + a[j1 + 1]; |
---|
| 1421 | x1r = a[j] - a[j1]; |
---|
| 1422 | x1i = a[j + 1] - a[j1 + 1]; |
---|
| 1423 | x2r = a[j2] + a[j3]; |
---|
| 1424 | x2i = a[j2 + 1] + a[j3 + 1]; |
---|
| 1425 | x3r = a[j2] - a[j3]; |
---|
| 1426 | x3i = a[j2 + 1] - a[j3 + 1]; |
---|
| 1427 | y0r = x0r + x2r; |
---|
| 1428 | y0i = x0i + x2i; |
---|
| 1429 | y2r = x0r - x2r; |
---|
| 1430 | y2i = x0i - x2i; |
---|
| 1431 | y1r = x1r - x3i; |
---|
| 1432 | y1i = x1i + x3r; |
---|
| 1433 | y3r = x1r + x3i; |
---|
| 1434 | y3i = x1i - x3r; |
---|
| 1435 | x0r = a[j4] + a[j5]; |
---|
| 1436 | x0i = a[j4 + 1] + a[j5 + 1]; |
---|
| 1437 | x1r = a[j4] - a[j5]; |
---|
| 1438 | x1i = a[j4 + 1] - a[j5 + 1]; |
---|
| 1439 | x2r = a[j6] + a[j7]; |
---|
| 1440 | x2i = a[j6 + 1] + a[j7 + 1]; |
---|
| 1441 | x3r = a[j6] - a[j7]; |
---|
| 1442 | x3i = a[j6 + 1] - a[j7 + 1]; |
---|
| 1443 | y4r = x0r + x2r; |
---|
| 1444 | y4i = x0i + x2i; |
---|
| 1445 | y6r = x0r - x2r; |
---|
| 1446 | y6i = x0i - x2i; |
---|
| 1447 | x0r = x1r - x3i; |
---|
| 1448 | x0i = x1i + x3r; |
---|
| 1449 | x2r = x1r + x3i; |
---|
| 1450 | x2i = x3r - x1i; |
---|
| 1451 | y5r = wk1i * x0r - wk1r * x0i; |
---|
| 1452 | y5i = wk1i * x0i + wk1r * x0r; |
---|
| 1453 | y7r = wk1r * x2r + wk1i * x2i; |
---|
| 1454 | y7i = wk1r * x2i - wk1i * x2r; |
---|
| 1455 | x0r = wk1r * y1r - wk1i * y1i; |
---|
| 1456 | x0i = wk1r * y1i + wk1i * y1r; |
---|
| 1457 | a[j1] = x0r + y5r; |
---|
| 1458 | a[j1 + 1] = x0i + y5i; |
---|
| 1459 | a[j5] = y5i - x0i; |
---|
| 1460 | a[j5 + 1] = x0r - y5r; |
---|
| 1461 | x0r = wk1i * y3r - wk1r * y3i; |
---|
| 1462 | x0i = wk1i * y3i + wk1r * y3r; |
---|
| 1463 | a[j3] = x0r - y7r; |
---|
| 1464 | a[j3 + 1] = x0i + y7i; |
---|
| 1465 | a[j7] = y7i - x0i; |
---|
| 1466 | a[j7 + 1] = x0r + y7r; |
---|
| 1467 | a[j] = y0r + y4r; |
---|
| 1468 | a[j + 1] = y0i + y4i; |
---|
| 1469 | a[j4] = y4i - y0i; |
---|
| 1470 | a[j4 + 1] = y0r - y4r; |
---|
| 1471 | x0r = y2r - y6i; |
---|
| 1472 | x0i = y2i + y6r; |
---|
| 1473 | a[j2] = wn4r * (x0r - x0i); |
---|
| 1474 | a[j2 + 1] = wn4r * (x0i + x0r); |
---|
| 1475 | x0r = y6r - y2i; |
---|
| 1476 | x0i = y2r + y6i; |
---|
| 1477 | a[j6] = wn4r * (x0r - x0i); |
---|
| 1478 | a[j6 + 1] = wn4r * (x0i + x0r); |
---|
| 1479 | } |
---|
| 1480 | k1 = 4; |
---|
| 1481 | for (k = 2 * m; k < n; k += m) { |
---|
| 1482 | k1 += 4; |
---|
| 1483 | wk1r = w[k1]; |
---|
| 1484 | wk1i = w[k1 + 1]; |
---|
| 1485 | wk2r = w[k1 + 2]; |
---|
| 1486 | wk2i = w[k1 + 3]; |
---|
| 1487 | wtmp = 2 * wk2i; |
---|
| 1488 | wk3r = wk1r - wtmp * wk1i; |
---|
| 1489 | wk3i = wtmp * wk1r - wk1i; |
---|
| 1490 | wk4r = 1 - wtmp * wk2i; |
---|
| 1491 | wk4i = wtmp * wk2r; |
---|
| 1492 | wtmp = 2 * wk4i; |
---|
| 1493 | wk5r = wk3r - wtmp * wk1i; |
---|
| 1494 | wk5i = wtmp * wk1r - wk3i; |
---|
| 1495 | wk6r = wk2r - wtmp * wk2i; |
---|
| 1496 | wk6i = wtmp * wk2r - wk2i; |
---|
| 1497 | wk7r = wk1r - wtmp * wk3i; |
---|
| 1498 | wk7i = wtmp * wk3r - wk1i; |
---|
| 1499 | for (j = k; j < l + k; j += 2) { |
---|
| 1500 | j1 = j + l; |
---|
| 1501 | j2 = j1 + l; |
---|
| 1502 | j3 = j2 + l; |
---|
| 1503 | j4 = j3 + l; |
---|
| 1504 | j5 = j4 + l; |
---|
| 1505 | j6 = j5 + l; |
---|
| 1506 | j7 = j6 + l; |
---|
| 1507 | x0r = a[j] + a[j1]; |
---|
| 1508 | x0i = a[j + 1] + a[j1 + 1]; |
---|
| 1509 | x1r = a[j] - a[j1]; |
---|
| 1510 | x1i = a[j + 1] - a[j1 + 1]; |
---|
| 1511 | x2r = a[j2] + a[j3]; |
---|
| 1512 | x2i = a[j2 + 1] + a[j3 + 1]; |
---|
| 1513 | x3r = a[j2] - a[j3]; |
---|
| 1514 | x3i = a[j2 + 1] - a[j3 + 1]; |
---|
| 1515 | y0r = x0r + x2r; |
---|
| 1516 | y0i = x0i + x2i; |
---|
| 1517 | y2r = x0r - x2r; |
---|
| 1518 | y2i = x0i - x2i; |
---|
| 1519 | y1r = x1r - x3i; |
---|
| 1520 | y1i = x1i + x3r; |
---|
| 1521 | y3r = x1r + x3i; |
---|
| 1522 | y3i = x1i - x3r; |
---|
| 1523 | x0r = a[j4] + a[j5]; |
---|
| 1524 | x0i = a[j4 + 1] + a[j5 + 1]; |
---|
| 1525 | x1r = a[j4] - a[j5]; |
---|
| 1526 | x1i = a[j4 + 1] - a[j5 + 1]; |
---|
| 1527 | x2r = a[j6] + a[j7]; |
---|
| 1528 | x2i = a[j6 + 1] + a[j7 + 1]; |
---|
| 1529 | x3r = a[j6] - a[j7]; |
---|
| 1530 | x3i = a[j6 + 1] - a[j7 + 1]; |
---|
| 1531 | y4r = x0r + x2r; |
---|
| 1532 | y4i = x0i + x2i; |
---|
| 1533 | y6r = x0r - x2r; |
---|
| 1534 | y6i = x0i - x2i; |
---|
| 1535 | x0r = x1r - x3i; |
---|
| 1536 | x0i = x1i + x3r; |
---|
| 1537 | x2r = x1r + x3i; |
---|
| 1538 | x2i = x1i - x3r; |
---|
| 1539 | y5r = wn4r * (x0r - x0i); |
---|
| 1540 | y5i = wn4r * (x0r + x0i); |
---|
| 1541 | y7r = wn4r * (x2r - x2i); |
---|
| 1542 | y7i = wn4r * (x2r + x2i); |
---|
| 1543 | x0r = y1r + y5r; |
---|
| 1544 | x0i = y1i + y5i; |
---|
| 1545 | a[j1] = wk1r * x0r - wk1i * x0i; |
---|
| 1546 | a[j1 + 1] = wk1r * x0i + wk1i * x0r; |
---|
| 1547 | x0r = y1r - y5r; |
---|
| 1548 | x0i = y1i - y5i; |
---|
| 1549 | a[j5] = wk5r * x0r - wk5i * x0i; |
---|
| 1550 | a[j5 + 1] = wk5r * x0i + wk5i * x0r; |
---|
| 1551 | x0r = y3r - y7i; |
---|
| 1552 | x0i = y3i + y7r; |
---|
| 1553 | a[j3] = wk3r * x0r - wk3i * x0i; |
---|
| 1554 | a[j3 + 1] = wk3r * x0i + wk3i * x0r; |
---|
| 1555 | x0r = y3r + y7i; |
---|
| 1556 | x0i = y3i - y7r; |
---|
| 1557 | a[j7] = wk7r * x0r - wk7i * x0i; |
---|
| 1558 | a[j7 + 1] = wk7r * x0i + wk7i * x0r; |
---|
| 1559 | a[j] = y0r + y4r; |
---|
| 1560 | a[j + 1] = y0i + y4i; |
---|
| 1561 | x0r = y0r - y4r; |
---|
| 1562 | x0i = y0i - y4i; |
---|
| 1563 | a[j4] = wk4r * x0r - wk4i * x0i; |
---|
| 1564 | a[j4 + 1] = wk4r * x0i + wk4i * x0r; |
---|
| 1565 | x0r = y2r - y6i; |
---|
| 1566 | x0i = y2i + y6r; |
---|
| 1567 | a[j2] = wk2r * x0r - wk2i * x0i; |
---|
| 1568 | a[j2 + 1] = wk2r * x0i + wk2i * x0r; |
---|
| 1569 | x0r = y2r + y6i; |
---|
| 1570 | x0i = y2i - y6r; |
---|
| 1571 | a[j6] = wk6r * x0r - wk6i * x0i; |
---|
| 1572 | a[j6 + 1] = wk6r * x0i + wk6i * x0r; |
---|
| 1573 | } |
---|
| 1574 | } |
---|
| 1575 | } |
---|
| 1576 | } |
---|
| 1577 | |
---|
| 1578 | |
---|
[5c6b264] | 1579 | void rftfsub(int n, smpl_t *a, int nc, smpl_t *c) |
---|
[cfaa3c4] | 1580 | { |
---|
| 1581 | int j, k, kk, ks, m; |
---|
[5c6b264] | 1582 | smpl_t wkr, wki, xr, xi, yr, yi; |
---|
[cfaa3c4] | 1583 | |
---|
| 1584 | m = n >> 1; |
---|
| 1585 | ks = 2 * nc / m; |
---|
| 1586 | kk = 0; |
---|
| 1587 | for (j = 2; j < m; j += 2) { |
---|
| 1588 | k = n - j; |
---|
| 1589 | kk += ks; |
---|
| 1590 | wkr = 0.5 - c[nc - kk]; |
---|
| 1591 | wki = c[kk]; |
---|
| 1592 | xr = a[j] - a[k]; |
---|
| 1593 | xi = a[j + 1] + a[k + 1]; |
---|
| 1594 | yr = wkr * xr - wki * xi; |
---|
| 1595 | yi = wkr * xi + wki * xr; |
---|
| 1596 | a[j] -= yr; |
---|
| 1597 | a[j + 1] -= yi; |
---|
| 1598 | a[k] += yr; |
---|
| 1599 | a[k + 1] -= yi; |
---|
| 1600 | } |
---|
| 1601 | } |
---|
| 1602 | |
---|
| 1603 | |
---|
[5c6b264] | 1604 | void rftbsub(int n, smpl_t *a, int nc, smpl_t *c) |
---|
[cfaa3c4] | 1605 | { |
---|
| 1606 | int j, k, kk, ks, m; |
---|
[5c6b264] | 1607 | smpl_t wkr, wki, xr, xi, yr, yi; |
---|
[cfaa3c4] | 1608 | |
---|
| 1609 | a[1] = -a[1]; |
---|
| 1610 | m = n >> 1; |
---|
| 1611 | ks = 2 * nc / m; |
---|
| 1612 | kk = 0; |
---|
| 1613 | for (j = 2; j < m; j += 2) { |
---|
| 1614 | k = n - j; |
---|
| 1615 | kk += ks; |
---|
| 1616 | wkr = 0.5 - c[nc - kk]; |
---|
| 1617 | wki = c[kk]; |
---|
| 1618 | xr = a[j] - a[k]; |
---|
| 1619 | xi = a[j + 1] + a[k + 1]; |
---|
| 1620 | yr = wkr * xr + wki * xi; |
---|
| 1621 | yi = wkr * xi - wki * xr; |
---|
| 1622 | a[j] -= yr; |
---|
| 1623 | a[j + 1] = yi - a[j + 1]; |
---|
| 1624 | a[k] += yr; |
---|
| 1625 | a[k + 1] = yi - a[k + 1]; |
---|
| 1626 | } |
---|
| 1627 | a[m + 1] = -a[m + 1]; |
---|
| 1628 | } |
---|
| 1629 | |
---|
| 1630 | |
---|
[5c6b264] | 1631 | void dctsub(int n, smpl_t *a, int nc, smpl_t *c) |
---|
[cfaa3c4] | 1632 | { |
---|
| 1633 | int j, k, kk, ks, m; |
---|
[5c6b264] | 1634 | smpl_t wkr, wki, xr; |
---|
[cfaa3c4] | 1635 | |
---|
| 1636 | m = n >> 1; |
---|
| 1637 | ks = nc / n; |
---|
| 1638 | kk = 0; |
---|
| 1639 | for (j = 1; j < m; j++) { |
---|
| 1640 | k = n - j; |
---|
| 1641 | kk += ks; |
---|
| 1642 | wkr = c[kk] - c[nc - kk]; |
---|
| 1643 | wki = c[kk] + c[nc - kk]; |
---|
| 1644 | xr = wki * a[j] - wkr * a[k]; |
---|
| 1645 | a[j] = wkr * a[j] + wki * a[k]; |
---|
| 1646 | a[k] = xr; |
---|
| 1647 | } |
---|
| 1648 | a[m] *= c[0]; |
---|
| 1649 | } |
---|
| 1650 | |
---|
| 1651 | |
---|
[5c6b264] | 1652 | void dstsub(int n, smpl_t *a, int nc, smpl_t *c) |
---|
[cfaa3c4] | 1653 | { |
---|
| 1654 | int j, k, kk, ks, m; |
---|
[5c6b264] | 1655 | smpl_t wkr, wki, xr; |
---|
[cfaa3c4] | 1656 | |
---|
| 1657 | m = n >> 1; |
---|
| 1658 | ks = nc / n; |
---|
| 1659 | kk = 0; |
---|
| 1660 | for (j = 1; j < m; j++) { |
---|
| 1661 | k = n - j; |
---|
| 1662 | kk += ks; |
---|
| 1663 | wkr = c[kk] - c[nc - kk]; |
---|
| 1664 | wki = c[kk] + c[nc - kk]; |
---|
| 1665 | xr = wki * a[k] - wkr * a[j]; |
---|
| 1666 | a[k] = wkr * a[k] + wki * a[j]; |
---|
| 1667 | a[j] = xr; |
---|
| 1668 | } |
---|
| 1669 | a[m] *= c[0]; |
---|
| 1670 | } |
---|
| 1671 | |
---|