1e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent/* 2e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent * http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html 3e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent * Copyright Takuya OOURA, 1996-2001 4e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent * 5e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent * You may use, copy, modify and distribute this code for any purpose (include 6e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent * commercial use) and without fee. Please refer to this package when you modify 7e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent * this code. 8e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent * 9e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent * Changes: 10e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent * Trivial type modifications by the WebRTC authors. 11e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent */ 12e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 13e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent/* 14e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric LaurentFast Fourier/Cosine/Sine Transform 15e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dimension :one 16e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent data length :power of 2 17e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent decimation :frequency 18e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent radix :4, 2 19e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent data :inplace 20e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent table :use 21e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurentfunctions 22e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cdft: Complex Discrete Fourier Transform 23e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rdft: Real Discrete Fourier Transform 24e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ddct: Discrete Cosine Transform 25e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ddst: Discrete Sine Transform 26e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dfct: Cosine Transform of RDFT (Real Symmetric DFT) 27e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dfst: Sine Transform of RDFT (Real Anti-symmetric DFT) 28e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurentfunction prototypes 29e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent void cdft(int, int, float *, int *, float *); 30e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent void rdft(int, int, float *, int *, float *); 31e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent void ddct(int, int, float *, int *, float *); 32e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent void ddst(int, int, float *, int *, float *); 33e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent void dfct(int, float *, float *, int *, float *); 34e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent void dfst(int, float *, float *, int *, float *); 35e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 36e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 37e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent-------- Complex DFT (Discrete Fourier Transform) -------- 38e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [definition] 39e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case1> 40e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n 41e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case2> 42e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n 43e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent (notes: sum_j=0^n-1 is a summation from j=0 to n-1) 44e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [usage] 45e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case1> 46e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; // first time only 47e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cdft(2*n, 1, a, ip, w); 48e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case2> 49e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; // first time only 50e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cdft(2*n, -1, a, ip, w); 51e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [parameters] 52e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 2*n :data length (int) 53e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n >= 1, n = power of 2 54e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0...2*n-1] :input/output data (float *) 55e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent input data 56e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2*j] = Re(x[j]), 57e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2*j+1] = Im(x[j]), 0<=j<n 58e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent output data 59e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2*k] = Re(X[k]), 60e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2*k+1] = Im(X[k]), 0<=k<n 61e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0...*] :work area for bit reversal (int *) 62e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 2+sqrt(n) 63e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent strictly, 64e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 65e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 2+(1<<(int)(log(n+0.5)/log(2))/2). 66e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0],ip[1] are pointers of the cos/sin table. 67e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[0...n/2-1] :cos/sin table (float *) 68e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[],ip[] are initialized if ip[0] == 0. 69e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [remark] 70e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent Inverse of 71e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cdft(2*n, -1, a, ip, w); 72e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent is 73e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cdft(2*n, 1, a, ip, w); 74e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j <= 2 * n - 1; j++) { 75e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] *= 1.0 / n; 76e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 77e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent . 78e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 79e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 80e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent-------- Real DFT / Inverse of Real DFT -------- 81e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [definition] 82e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case1> RDFT 83e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2 84e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2 85e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case2> IRDFT (excluding scale) 86e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 87e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 88e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n 89e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [usage] 90e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case1> 91e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; // first time only 92e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rdft(n, 1, a, ip, w); 93e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case2> 94e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; // first time only 95e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rdft(n, -1, a, ip, w); 96e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [parameters] 97e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n :data length (int) 98e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n >= 2, n = power of 2 99e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0...n-1] :input/output data (float *) 100e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case1> 101e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent output data 102e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2*k] = R[k], 0<=k<n/2 103e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2*k+1] = I[k], 0<k<n/2 104e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = R[n/2] 105e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case2> 106e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent input data 107e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2*j] = R[j], 0<=j<n/2 108e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2*j+1] = I[j], 0<j<n/2 109e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = R[n/2] 110e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0...*] :work area for bit reversal (int *) 111e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 2+sqrt(n/2) 112e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent strictly, 113e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 114e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 2+(1<<(int)(log(n/2+0.5)/log(2))/2). 115e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0],ip[1] are pointers of the cos/sin table. 116e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[0...n/2-1] :cos/sin table (float *) 117e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[],ip[] are initialized if ip[0] == 0. 118e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [remark] 119e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent Inverse of 120e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rdft(n, 1, a, ip, w); 121e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent is 122e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rdft(n, -1, a, ip, w); 123e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j <= n - 1; j++) { 124e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] *= 2.0 / n; 125e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 126e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent . 127e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 128e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 129e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent-------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- 130e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [definition] 131e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case1> IDCT (excluding scale) 132e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n 133e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case2> DCT 134e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n 135e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [usage] 136e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case1> 137e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; // first time only 138e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ddct(n, 1, a, ip, w); 139e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case2> 140e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; // first time only 141e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ddct(n, -1, a, ip, w); 142e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [parameters] 143e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n :data length (int) 144e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n >= 2, n = power of 2 145e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0...n-1] :input/output data (float *) 146e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent output data 147e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] = C[k], 0<=k<n 148e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0...*] :work area for bit reversal (int *) 149e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 2+sqrt(n/2) 150e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent strictly, 151e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 152e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 2+(1<<(int)(log(n/2+0.5)/log(2))/2). 153e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0],ip[1] are pointers of the cos/sin table. 154e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[0...n*5/4-1] :cos/sin table (float *) 155e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[],ip[] are initialized if ip[0] == 0. 156e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [remark] 157e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent Inverse of 158e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ddct(n, -1, a, ip, w); 159e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent is 160e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] *= 0.5; 161e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ddct(n, 1, a, ip, w); 162e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j <= n - 1; j++) { 163e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] *= 2.0 / n; 164e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 165e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent . 166e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 167e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 168e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent-------- DST (Discrete Sine Transform) / Inverse of DST -------- 169e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [definition] 170e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case1> IDST (excluding scale) 171e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n 172e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case2> DST 173e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n 174e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [usage] 175e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case1> 176e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; // first time only 177e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ddst(n, 1, a, ip, w); 178e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case2> 179e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; // first time only 180e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ddst(n, -1, a, ip, w); 181e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [parameters] 182e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n :data length (int) 183e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n >= 2, n = power of 2 184e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0...n-1] :input/output data (float *) 185e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case1> 186e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent input data 187e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = A[j], 0<j<n 188e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] = A[n] 189e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent output data 190e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] = S[k], 0<=k<n 191e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent <case2> 192e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent output data 193e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] = S[k], 0<k<n 194e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] = S[n] 195e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0...*] :work area for bit reversal (int *) 196e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 2+sqrt(n/2) 197e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent strictly, 198e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 199e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 2+(1<<(int)(log(n/2+0.5)/log(2))/2). 200e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0],ip[1] are pointers of the cos/sin table. 201e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[0...n*5/4-1] :cos/sin table (float *) 202e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[],ip[] are initialized if ip[0] == 0. 203e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [remark] 204e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent Inverse of 205e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ddst(n, -1, a, ip, w); 206e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent is 207e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] *= 0.5; 208e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ddst(n, 1, a, ip, w); 209e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j <= n - 1; j++) { 210e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] *= 2.0 / n; 211e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 212e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent . 213e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 214e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 215e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent-------- Cosine Transform of RDFT (Real Symmetric DFT) -------- 216e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [definition] 217e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n 218e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [usage] 219e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; // first time only 220e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dfct(n, a, t, ip, w); 221e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [parameters] 222e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n :data length - 1 (int) 223e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n >= 2, n = power of 2 224e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0...n] :input/output data (float *) 225e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent output data 226e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] = C[k], 0<=k<=n 227e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[0...n/2] :work area (float *) 228e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0...*] :work area for bit reversal (int *) 229e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 2+sqrt(n/4) 230e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent strictly, 231e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 232e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 2+(1<<(int)(log(n/4+0.5)/log(2))/2). 233e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0],ip[1] are pointers of the cos/sin table. 234e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[0...n*5/8-1] :cos/sin table (float *) 235e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[],ip[] are initialized if ip[0] == 0. 236e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [remark] 237e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent Inverse of 238e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] *= 0.5; 239e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[n] *= 0.5; 240e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dfct(n, a, t, ip, w); 241e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent is 242e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] *= 0.5; 243e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[n] *= 0.5; 244e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dfct(n, a, t, ip, w); 245e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j <= n; j++) { 246e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] *= 2.0 / n; 247e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 248e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent . 249e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 250e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 251e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent-------- Sine Transform of RDFT (Real Anti-symmetric DFT) -------- 252e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [definition] 253e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n 254e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [usage] 255e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; // first time only 256e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dfst(n, a, t, ip, w); 257e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [parameters] 258e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n :data length + 1 (int) 259e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent n >= 2, n = power of 2 260e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0...n-1] :input/output data (float *) 261e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent output data 262e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] = S[k], 0<k<n 263e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent (a[0] is used for work area) 264e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[0...n/2-1] :work area (float *) 265e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0...*] :work area for bit reversal (int *) 266e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 2+sqrt(n/4) 267e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent strictly, 268e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent length of ip >= 269e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 2+(1<<(int)(log(n/4+0.5)/log(2))/2). 270e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0],ip[1] are pointers of the cos/sin table. 271e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[0...n*5/8-1] :cos/sin table (float *) 272e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[],ip[] are initialized if ip[0] == 0. 273e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent [remark] 274e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent Inverse of 275e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dfst(n, a, t, ip, w); 276e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent is 277e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dfst(n, a, t, ip, w); 278e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 1; j <= n - 1; j++) { 279e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] *= 2.0 / n; 280e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 281e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent . 282e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 283e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 284e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric LaurentAppendix : 285e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent The cos/sin table is recalculated when the larger table required. 286e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[] and ip[] are compatible with all routines. 287e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent*/ 288e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 289c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void makewt(int nw, int *ip, float *w); 290c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void makect(int nc, int *ip, float *c); 291c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void bitrv2(int n, int *ip, float *a); 292c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void bitrv2conj(int n, int *ip, float *a); 293c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void cftfsub(int n, float *a, float *w); 294c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void cftbsub(int n, float *a, float *w); 295c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void cft1st(int n, float *a, float *w); 296c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void cftmdl(int n, int l, float *a, float *w); 297c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void rftfsub(int n, float *a, int nc, float *c); 298c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void rftbsub(int n, float *a, int nc, float *c); 299c55a96383497a772a307b346368133960b02ad03Eric Laurent#if 0 // Not used. 300c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void dctsub(int n, float *a, int nc, float *c) 301c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void dstsub(int n, float *a, int nc, float *c) 302c55a96383497a772a307b346368133960b02ad03Eric Laurent#endif 303c55a96383497a772a307b346368133960b02ad03Eric Laurent 304c55a96383497a772a307b346368133960b02ad03Eric Laurent 305c55a96383497a772a307b346368133960b02ad03Eric Laurentvoid WebRtc_cdft(int n, int isgn, float *a, int *ip, float *w) 306e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 307e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > (ip[0] << 2)) { 308e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makewt(n >> 2, ip, w); 309e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 310e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 4) { 311e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (isgn >= 0) { 312e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(n, ip + 2, a); 313e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 314e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else { 315e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2conj(n, ip + 2, a); 316e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftbsub(n, a, w); 317e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 318e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (n == 4) { 319e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 320e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 321e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 322e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 323e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 324c55a96383497a772a307b346368133960b02ad03Eric Laurentvoid WebRtc_rdft(int n, int isgn, float *a, int *ip, float *w) 325e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 326e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int nw, nc; 327e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float xi; 328e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 329e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nw = ip[0]; 330e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > (nw << 2)) { 331e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nw = n >> 2; 332e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makewt(nw, ip, w); 333e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 334e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nc = ip[1]; 335e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > (nc << 2)) { 336e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nc = n >> 2; 337e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makect(nc, ip, w + nw); 338e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 339e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (isgn >= 0) { 340e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 4) { 341e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(n, ip + 2, a); 342e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 343e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rftfsub(n, a, nc, w + nw); 344e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (n == 4) { 345e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 346e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 347e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[0] - a[1]; 348e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] += a[1]; 349e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = xi; 350e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else { 351e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = 0.5f * (a[0] - a[1]); 352e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] -= a[1]; 353e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 4) { 354e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rftbsub(n, a, nc, w + nw); 355e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(n, ip + 2, a); 356e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftbsub(n, a, w); 357e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (n == 4) { 358e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 359e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 360e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 361e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 362e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 363c55a96383497a772a307b346368133960b02ad03Eric Laurent#if 0 // Not used. 364c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void ddct(int n, int isgn, float *a, int *ip, float *w) 365e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 366e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, nw, nc; 367e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float xr; 368e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 369e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nw = ip[0]; 370e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > (nw << 2)) { 371e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nw = n >> 2; 372e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makewt(nw, ip, w); 373e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 374e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nc = ip[1]; 375e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > nc) { 376e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nc = n; 377e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makect(nc, ip, w + nw); 378e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 379e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (isgn < 0) { 380e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[n - 1]; 381e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = n - 2; j >= 2; j -= 2) { 382e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = a[j] - a[j - 1]; 383e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] += a[j - 1]; 384e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 385e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = a[0] - xr; 386e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] += xr; 387e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 4) { 388e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rftbsub(n, a, nc, w + nw); 389e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(n, ip + 2, a); 390e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftbsub(n, a, w); 391e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (n == 4) { 392e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 393e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 394e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 395e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dctsub(n, a, nc, w + nw); 396e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (isgn >= 0) { 397e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 4) { 398e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(n, ip + 2, a); 399e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 400e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rftfsub(n, a, nc, w + nw); 401e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (n == 4) { 402e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 403e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 404e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[0] - a[1]; 405e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] += a[1]; 406e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 2; j < n; j += 2) { 407e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j - 1] = a[j] - a[j + 1]; 408e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] += a[j + 1]; 409e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 410e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[n - 1] = xr; 411e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 412e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 413e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 414e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 415c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void ddst(int n, int isgn, float *a, int *ip, float *w) 416e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 417e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, nw, nc; 418e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float xr; 419e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 420e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nw = ip[0]; 421e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > (nw << 2)) { 422e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nw = n >> 2; 423e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makewt(nw, ip, w); 424e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 425e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nc = ip[1]; 426e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > nc) { 427e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nc = n; 428e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makect(nc, ip, w + nw); 429e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 430e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (isgn < 0) { 431e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[n - 1]; 432e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = n - 2; j >= 2; j -= 2) { 433e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = -a[j] - a[j - 1]; 434e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] -= a[j - 1]; 435e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 436e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = a[0] + xr; 437e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] -= xr; 438e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 4) { 439e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rftbsub(n, a, nc, w + nw); 440e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(n, ip + 2, a); 441e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftbsub(n, a, w); 442e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (n == 4) { 443e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 444e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 445e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 446e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dstsub(n, a, nc, w + nw); 447e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (isgn >= 0) { 448e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 4) { 449e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(n, ip + 2, a); 450e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 451e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rftfsub(n, a, nc, w + nw); 452e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (n == 4) { 453e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(n, a, w); 454e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 455e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[0] - a[1]; 456e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] += a[1]; 457e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 2; j < n; j += 2) { 458e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j - 1] = -a[j] - a[j + 1]; 459e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] -= a[j + 1]; 460e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 461e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[n - 1] = -xr; 462e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 463e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 464e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 465e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 466c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void dfct(int n, float *a, float *t, int *ip, float *w) 467e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 468e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, k, l, m, mh, nw, nc; 469e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float xr, xi, yr, yi; 470e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 471e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nw = ip[0]; 472e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > (nw << 3)) { 473e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nw = n >> 3; 474e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makewt(nw, ip, w); 475e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 476e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nc = ip[1]; 477e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > (nc << 1)) { 478e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nc = n >> 1; 479e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makect(nc, ip, w + nw); 480e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 481e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = n >> 1; 482e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = a[m]; 483e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[0] + a[n]; 484e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] -= a[n]; 485e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[0] = xi - yi; 486e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[m] = xi + yi; 487e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 2) { 488e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent mh = m >> 1; 489e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 1; j < mh; j++) { 490e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k = m - j; 491e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j] - a[n - j]; 492e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j] + a[n - j]; 493e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k] - a[n - k]; 494e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = a[k] + a[n - k]; 495e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = xr; 496e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] = yr; 497e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[j] = xi - yi; 498e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[k] = xi + yi; 499e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 500e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[mh] = a[mh] + a[n - mh]; 501e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[mh] -= a[n - mh]; 502e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dctsub(m, a, nc, w + nw); 503e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (m > 4) { 504e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(m, ip + 2, a); 505e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(m, a, w); 506e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rftfsub(m, a, nc, w + nw); 507e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (m == 4) { 508e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(m, a, w); 509e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 510e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[n - 1] = a[0] - a[1]; 511e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = a[0] + a[1]; 512e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = m - 2; j >= 2; j -= 2) { 513e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2 * j + 1] = a[j] + a[j + 1]; 514e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2 * j - 1] = a[j] - a[j + 1]; 515e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 516e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l = 2; 517e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = mh; 518e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent while (m >= 2) { 519e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dctsub(m, t, nc, w + nw); 520e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (m > 4) { 521e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(m, ip + 2, t); 522e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(m, t, w); 523e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rftfsub(m, t, nc, w + nw); 524e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (m == 4) { 525e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(m, t, w); 526e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 527e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[n - l] = t[0] - t[1]; 528e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[l] = t[0] + t[1]; 529e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k = 0; 530e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 2; j < m; j += 2) { 531e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k += l << 2; 532e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k - l] = t[j] - t[j + 1]; 533e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k + l] = t[j] + t[j + 1]; 534e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 535e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l <<= 1; 536e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent mh = m >> 1; 537e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < mh; j++) { 538e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k = m - j; 539e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[j] = t[m + k] - t[m + j]; 540e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[k] = t[m + k] + t[m + j]; 541e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 542e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[mh] = t[m + mh]; 543e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = mh; 544e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 545e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[l] = t[0]; 546e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[n] = t[2] - t[1]; 547e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] = t[2] + t[1]; 548e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else { 549e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = a[0]; 550e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2] = t[0]; 551e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] = t[1]; 552e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 553e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 554e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 555c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void dfst(int n, float *a, float *t, int *ip, float *w) 556e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 557e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, k, l, m, mh, nw, nc; 558e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float xr, xi, yr, yi; 559e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 560e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nw = ip[0]; 561e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > (nw << 3)) { 562e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nw = n >> 3; 563e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makewt(nw, ip, w); 564e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 565e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nc = ip[1]; 566e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > (nc << 1)) { 567e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nc = n >> 1; 568e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent makect(nc, ip, w + nw); 569e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 570e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 2) { 571e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = n >> 1; 572e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent mh = m >> 1; 573e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 1; j < mh; j++) { 574e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k = m - j; 575e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j] + a[n - j]; 576e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j] - a[n - j]; 577e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k] + a[n - k]; 578e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = a[k] - a[n - k]; 579e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = xr; 580e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] = yr; 581e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[j] = xi + yi; 582e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[k] = xi - yi; 583e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 584e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[0] = a[mh] - a[n - mh]; 585e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[mh] += a[n - mh]; 586e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] = a[m]; 587e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dstsub(m, a, nc, w + nw); 588e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (m > 4) { 589e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(m, ip + 2, a); 590e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(m, a, w); 591e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rftfsub(m, a, nc, w + nw); 592e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (m == 4) { 593e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(m, a, w); 594e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 595e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[n - 1] = a[1] - a[0]; 596e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = a[0] + a[1]; 597e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = m - 2; j >= 2; j -= 2) { 598e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2 * j + 1] = a[j] - a[j + 1]; 599e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2 * j - 1] = -a[j] - a[j + 1]; 600e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 601e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l = 2; 602e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = mh; 603e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent while (m >= 2) { 604e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent dstsub(m, t, nc, w + nw); 605e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (m > 4) { 606e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(m, ip + 2, t); 607e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(m, t, w); 608e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent rftfsub(m, t, nc, w + nw); 609e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else if (m == 4) { 610e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftfsub(m, t, w); 611e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 612e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[n - l] = t[1] - t[0]; 613e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[l] = t[0] + t[1]; 614e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k = 0; 615e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 2; j < m; j += 2) { 616e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k += l << 2; 617e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k - l] = -t[j] - t[j + 1]; 618e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k + l] = t[j] - t[j + 1]; 619e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 620e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l <<= 1; 621e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent mh = m >> 1; 622e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 1; j < mh; j++) { 623e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k = m - j; 624e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[j] = t[m + k] + t[m + j]; 625e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[k] = t[m + k] - t[m + j]; 626e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 627e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent t[0] = t[m + mh]; 628e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = mh; 629e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 630e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[l] = t[0]; 631e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 632e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] = 0; 633e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 634c55a96383497a772a307b346368133960b02ad03Eric Laurent#endif // Not used. 635e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 636e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 637e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent/* -------- initializing routines -------- */ 638e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 639e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 640e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent#include <math.h> 641e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 642c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void makewt(int nw, int *ip, float *w) 643e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 644e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, nwh; 645e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float delta, x, y; 646e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 647e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = nw; 648e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[1] = 1; 649e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (nw > 2) { 650e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nwh = nw >> 1; 651e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent delta = (float)atan(1.0f) / nwh; 652e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[0] = 1; 653e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[1] = 0; 654e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[nwh] = (float)cos(delta * nwh); 655e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[nwh + 1] = w[nwh]; 656e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (nwh > 2) { 657e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 2; j < nwh; j += 2) { 658e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x = (float)cos(delta * j); 659e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent y = (float)sin(delta * j); 660e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[j] = x; 661e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[j + 1] = y; 662e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[nw - j] = y; 663e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent w[nw - j + 1] = x; 664e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 665e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent bitrv2(nw, ip + 2, w); 666e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 667e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 668e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 669e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 670e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 671c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void makect(int nc, int *ip, float *c) 672e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 673e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, nch; 674e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float delta; 675e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 676e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[1] = nc; 677e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (nc > 1) { 678e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent nch = nc >> 1; 679e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent delta = (float)atan(1.0f) / nch; 680e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent c[0] = (float)cos(delta * nch); 681e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent c[nch] = 0.5f * c[0]; 682e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 1; j < nch; j++) { 683e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent c[j] = 0.5f * (float)cos(delta * j); 684e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent c[nc - j] = 0.5f * (float)sin(delta * j); 685e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 686e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 687e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 688e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 689e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 690e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent/* -------- child routines -------- */ 691e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 692e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 693c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void bitrv2(int n, int *ip, float *a) 694e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 695e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, j1, k, k1, l, m, m2; 696e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float xr, xi, yr, yi; 697e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 698e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; 699e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l = n; 700e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = 1; 701e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent while ((m << 3) < l) { 702e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l >>= 1; 703e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < m; j++) { 704e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[m + j] = ip[j] + l; 705e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 706e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m <<= 1; 707e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 708e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m2 = 2 * m; 709e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if ((m << 3) == l) { 710e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (k = 0; k < m; k++) { 711e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < k; j++) { 712e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = 2 * j + ip[k]; 713e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 = 2 * k + ip[j]; 714e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 715e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j1 + 1]; 716e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 717e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = a[k1 + 1]; 718e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 719e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 720e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 721e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 722e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 += m2; 723e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 += 2 * m2; 724e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 725e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j1 + 1]; 726e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 727e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = a[k1 + 1]; 728e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 729e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 730e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 731e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 732e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 += m2; 733e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 -= m2; 734e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 735e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j1 + 1]; 736e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 737e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = a[k1 + 1]; 738e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 739e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 740e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 741e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 742e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 += m2; 743e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 += 2 * m2; 744e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 745e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j1 + 1]; 746e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 747e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = a[k1 + 1]; 748e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 749e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 750e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 751e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 752e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 753e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = 2 * k + m2 + ip[k]; 754e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 = j1 + m2; 755e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 756e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j1 + 1]; 757e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 758e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = a[k1 + 1]; 759e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 760e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 761e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 762e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 763e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 764e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else { 765e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (k = 1; k < m; k++) { 766e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < k; j++) { 767e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = 2 * j + ip[k]; 768e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 = 2 * k + ip[j]; 769e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 770e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j1 + 1]; 771e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 772e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = a[k1 + 1]; 773e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 774e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 775e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 776e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 777e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 += m2; 778e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 += m2; 779e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 780e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j1 + 1]; 781e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 782e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = a[k1 + 1]; 783e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 784e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 785e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 786e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 787e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 788e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 789e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 790e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 791e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 792e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 793c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void bitrv2conj(int n, int *ip, float *a) 794e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 795e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, j1, k, k1, l, m, m2; 796e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float xr, xi, yr, yi; 797e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 798e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[0] = 0; 799e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l = n; 800e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = 1; 801e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent while ((m << 3) < l) { 802e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l >>= 1; 803e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < m; j++) { 804e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ip[m + j] = ip[j] + l; 805e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 806e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m <<= 1; 807e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 808e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m2 = 2 * m; 809e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if ((m << 3) == l) { 810e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (k = 0; k < m; k++) { 811e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < k; j++) { 812e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = 2 * j + ip[k]; 813e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 = 2 * k + ip[j]; 814e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 815e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = -a[j1 + 1]; 816e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 817e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = -a[k1 + 1]; 818e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 819e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 820e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 821e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 822e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 += m2; 823e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 += 2 * m2; 824e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 825e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = -a[j1 + 1]; 826e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 827e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = -a[k1 + 1]; 828e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 829e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 830e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 831e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 832e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 += m2; 833e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 -= m2; 834e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 835e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = -a[j1 + 1]; 836e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 837e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = -a[k1 + 1]; 838e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 839e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 840e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 841e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 842e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 += m2; 843e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 += 2 * m2; 844e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 845e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = -a[j1 + 1]; 846e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 847e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = -a[k1 + 1]; 848e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 849e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 850e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 851e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 852e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 853e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 = 2 * k + ip[k]; 854e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = -a[k1 + 1]; 855e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = k1 + m2; 856e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 = j1 + m2; 857e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 858e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = -a[j1 + 1]; 859e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 860e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = -a[k1 + 1]; 861e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 862e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 863e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 864e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 865e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 += m2; 866e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = -a[k1 + 1]; 867e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 868e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else { 869e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = -a[1]; 870e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[m2 + 1] = -a[m2 + 1]; 871e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (k = 1; k < m; k++) { 872e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < k; j++) { 873e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = 2 * j + ip[k]; 874e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 = 2 * k + ip[j]; 875e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 876e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = -a[j1 + 1]; 877e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 878e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = -a[k1 + 1]; 879e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 880e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 881e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 882e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 883e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 += m2; 884e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 += m2; 885e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j1]; 886e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = -a[j1 + 1]; 887e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = a[k1]; 888e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = -a[k1 + 1]; 889e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = yr; 890e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = yi; 891e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1] = xr; 892e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = xi; 893e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 894e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 = 2 * k + ip[k]; 895e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + 1] = -a[k1 + 1]; 896e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k1 + m2 + 1] = -a[k1 + m2 + 1]; 897e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 898e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 899e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 900e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 901e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 902c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void cftfsub(int n, float *a, float *w) 903e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 904e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, j1, j2, j3, l; 905e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 906e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 907e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l = 2; 908e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 8) { 909e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cft1st(n, a, w); 910e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l = 8; 911e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent while ((l << 2) < n) { 912e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftmdl(n, l, a, w); 913e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l <<= 2; 914e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 915e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 916e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if ((l << 2) == n) { 917e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < l; j += 2) { 918e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = j + l; 919e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j2 = j1 + l; 920e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j3 = j2 + l; 921e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[j] + a[j1]; 922e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = a[j + 1] + a[j1 + 1]; 923e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1r = a[j] - a[j1]; 924e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1i = a[j + 1] - a[j1 + 1]; 925e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2r = a[j2] + a[j3]; 926e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2i = a[j2 + 1] + a[j3 + 1]; 927e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3r = a[j2] - a[j3]; 928e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3i = a[j2 + 1] - a[j3 + 1]; 929e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = x0r + x2r; 930e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = x0i + x2i; 931e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2] = x0r - x2r; 932e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2 + 1] = x0i - x2i; 933e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = x1r - x3i; 934e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = x1i + x3r; 935e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3] = x1r + x3i; 936e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3 + 1] = x1i - x3r; 937e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 938e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else { 939e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < l; j += 2) { 940e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = j + l; 941e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[j] - a[j1]; 942e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = a[j + 1] - a[j1 + 1]; 943e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] += a[j1]; 944e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] += a[j1 + 1]; 945e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = x0r; 946e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = x0i; 947e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 948e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 949e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 950e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 951e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 952c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void cftbsub(int n, float *a, float *w) 953e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 954e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, j1, j2, j3, l; 955e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 956e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 957e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l = 2; 958e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if (n > 8) { 959e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cft1st(n, a, w); 960e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l = 8; 961e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent while ((l << 2) < n) { 962e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent cftmdl(n, l, a, w); 963e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent l <<= 2; 964e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 965e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 966e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent if ((l << 2) == n) { 967e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < l; j += 2) { 968e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = j + l; 969e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j2 = j1 + l; 970e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j3 = j2 + l; 971e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[j] + a[j1]; 972e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = -a[j + 1] - a[j1 + 1]; 973e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1r = a[j] - a[j1]; 974e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1i = -a[j + 1] + a[j1 + 1]; 975e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2r = a[j2] + a[j3]; 976e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2i = a[j2 + 1] + a[j3 + 1]; 977e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3r = a[j2] - a[j3]; 978e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3i = a[j2 + 1] - a[j3 + 1]; 979e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = x0r + x2r; 980e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = x0i - x2i; 981e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2] = x0r - x2r; 982e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2 + 1] = x0i + x2i; 983e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = x1r - x3i; 984e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = x1i - x3r; 985e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3] = x1r + x3i; 986e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3 + 1] = x1i + x3r; 987e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 988e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } else { 989e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < l; j += 2) { 990e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = j + l; 991e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[j] - a[j1]; 992e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = -a[j + 1] + a[j1 + 1]; 993e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] += a[j1]; 994e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = -a[j + 1] - a[j1 + 1]; 995e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = x0r; 996e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = x0i; 997e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 998e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 999e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 1000e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1001e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1002c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void cft1st(int n, float *a, float *w) 1003e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 1004e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, k1, k2; 1005e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; 1006e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 1007e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1008e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[0] + a[2]; 1009e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = a[1] + a[3]; 1010e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1r = a[0] - a[2]; 1011e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1i = a[1] - a[3]; 1012e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2r = a[4] + a[6]; 1013e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2i = a[5] + a[7]; 1014e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3r = a[4] - a[6]; 1015e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3i = a[5] - a[7]; 1016e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[0] = x0r + x2r; 1017e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = x0i + x2i; 1018e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[4] = x0r - x2r; 1019e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[5] = x0i - x2i; 1020e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[2] = x1r - x3i; 1021e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[3] = x1i + x3r; 1022e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[6] = x1r + x3i; 1023e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[7] = x1i - x3r; 1024e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk1r = w[2]; 1025e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[8] + a[10]; 1026e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = a[9] + a[11]; 1027e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1r = a[8] - a[10]; 1028e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1i = a[9] - a[11]; 1029e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2r = a[12] + a[14]; 1030e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2i = a[13] + a[15]; 1031e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3r = a[12] - a[14]; 1032e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3i = a[13] - a[15]; 1033e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[8] = x0r + x2r; 1034e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[9] = x0i + x2i; 1035e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[12] = x2i - x0i; 1036e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[13] = x0r - x2r; 1037e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x1r - x3i; 1038e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x1i + x3r; 1039e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[10] = wk1r * (x0r - x0i); 1040e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[11] = wk1r * (x0r + x0i); 1041e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x3i + x1r; 1042e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x3r - x1i; 1043e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[14] = wk1r * (x0i - x0r); 1044e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[15] = wk1r * (x0i + x0r); 1045e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 = 0; 1046e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 16; j < n; j += 16) { 1047e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 += 2; 1048e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k2 = 2 * k1; 1049e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk2r = w[k1]; 1050e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk2i = w[k1 + 1]; 1051e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk1r = w[k2]; 1052e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk1i = w[k2 + 1]; 1053e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk3r = wk1r - 2 * wk2i * wk1i; 1054e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk3i = 2 * wk2i * wk1r - wk1i; 1055e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[j] + a[j + 2]; 1056e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = a[j + 1] + a[j + 3]; 1057e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1r = a[j] - a[j + 2]; 1058e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1i = a[j + 1] - a[j + 3]; 1059e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2r = a[j + 4] + a[j + 6]; 1060e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2i = a[j + 5] + a[j + 7]; 1061e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3r = a[j + 4] - a[j + 6]; 1062e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3i = a[j + 5] - a[j + 7]; 1063e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = x0r + x2r; 1064e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = x0i + x2i; 1065e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r -= x2r; 1066e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i -= x2i; 1067e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 4] = wk2r * x0r - wk2i * x0i; 1068e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 5] = wk2r * x0i + wk2i * x0r; 1069e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x1r - x3i; 1070e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x1i + x3r; 1071e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 2] = wk1r * x0r - wk1i * x0i; 1072e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 3] = wk1r * x0i + wk1i * x0r; 1073e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x1r + x3i; 1074e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x1i - x3r; 1075e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 6] = wk3r * x0r - wk3i * x0i; 1076e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 7] = wk3r * x0i + wk3i * x0r; 1077e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk1r = w[k2 + 2]; 1078e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk1i = w[k2 + 3]; 1079e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk3r = wk1r - 2 * wk2r * wk1i; 1080e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk3i = 2 * wk2r * wk1r - wk1i; 1081e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[j + 8] + a[j + 10]; 1082e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = a[j + 9] + a[j + 11]; 1083e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1r = a[j + 8] - a[j + 10]; 1084e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1i = a[j + 9] - a[j + 11]; 1085e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2r = a[j + 12] + a[j + 14]; 1086e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2i = a[j + 13] + a[j + 15]; 1087e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3r = a[j + 12] - a[j + 14]; 1088e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3i = a[j + 13] - a[j + 15]; 1089e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 8] = x0r + x2r; 1090e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 9] = x0i + x2i; 1091e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r -= x2r; 1092e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i -= x2i; 1093e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 12] = -wk2i * x0r - wk2r * x0i; 1094e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 13] = -wk2i * x0i + wk2r * x0r; 1095e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x1r - x3i; 1096e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x1i + x3r; 1097e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 10] = wk1r * x0r - wk1i * x0i; 1098e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 11] = wk1r * x0i + wk1i * x0r; 1099e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x1r + x3i; 1100e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x1i - x3r; 1101e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 14] = wk3r * x0r - wk3i * x0i; 1102e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 15] = wk3r * x0i + wk3i * x0r; 1103e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 1104e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 1105e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1106e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1107c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void cftmdl(int n, int l, float *a, float *w) 1108e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 1109e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, j1, j2, j3, k, k1, k2, m, m2; 1110e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; 1111e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 1112e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1113e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = l << 2; 1114e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 0; j < l; j += 2) { 1115e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = j + l; 1116e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j2 = j1 + l; 1117e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j3 = j2 + l; 1118e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[j] + a[j1]; 1119e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = a[j + 1] + a[j1 + 1]; 1120e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1r = a[j] - a[j1]; 1121e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1i = a[j + 1] - a[j1 + 1]; 1122e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2r = a[j2] + a[j3]; 1123e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2i = a[j2 + 1] + a[j3 + 1]; 1124e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3r = a[j2] - a[j3]; 1125e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3i = a[j2 + 1] - a[j3 + 1]; 1126e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = x0r + x2r; 1127e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = x0i + x2i; 1128e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2] = x0r - x2r; 1129e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2 + 1] = x0i - x2i; 1130e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = x1r - x3i; 1131e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = x1i + x3r; 1132e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3] = x1r + x3i; 1133e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3 + 1] = x1i - x3r; 1134e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 1135e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk1r = w[2]; 1136e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = m; j < l + m; j += 2) { 1137e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = j + l; 1138e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j2 = j1 + l; 1139e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j3 = j2 + l; 1140e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[j] + a[j1]; 1141e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = a[j + 1] + a[j1 + 1]; 1142e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1r = a[j] - a[j1]; 1143e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1i = a[j + 1] - a[j1 + 1]; 1144e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2r = a[j2] + a[j3]; 1145e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2i = a[j2 + 1] + a[j3 + 1]; 1146e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3r = a[j2] - a[j3]; 1147e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3i = a[j2 + 1] - a[j3 + 1]; 1148e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = x0r + x2r; 1149e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = x0i + x2i; 1150e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2] = x2i - x0i; 1151e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2 + 1] = x0r - x2r; 1152e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x1r - x3i; 1153e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x1i + x3r; 1154e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = wk1r * (x0r - x0i); 1155e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = wk1r * (x0r + x0i); 1156e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x3i + x1r; 1157e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x3r - x1i; 1158e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3] = wk1r * (x0i - x0r); 1159e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3 + 1] = wk1r * (x0i + x0r); 1160e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 1161e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 = 0; 1162e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m2 = 2 * m; 1163e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (k = m2; k < n; k += m2) { 1164e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k1 += 2; 1165e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k2 = 2 * k1; 1166e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk2r = w[k1]; 1167e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk2i = w[k1 + 1]; 1168e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk1r = w[k2]; 1169e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk1i = w[k2 + 1]; 1170e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk3r = wk1r - 2 * wk2i * wk1i; 1171e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk3i = 2 * wk2i * wk1r - wk1i; 1172e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = k; j < l + k; j += 2) { 1173e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = j + l; 1174e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j2 = j1 + l; 1175e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j3 = j2 + l; 1176e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[j] + a[j1]; 1177e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = a[j + 1] + a[j1 + 1]; 1178e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1r = a[j] - a[j1]; 1179e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1i = a[j + 1] - a[j1 + 1]; 1180e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2r = a[j2] + a[j3]; 1181e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2i = a[j2 + 1] + a[j3 + 1]; 1182e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3r = a[j2] - a[j3]; 1183e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3i = a[j2 + 1] - a[j3 + 1]; 1184e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = x0r + x2r; 1185e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = x0i + x2i; 1186e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r -= x2r; 1187e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i -= x2i; 1188e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2] = wk2r * x0r - wk2i * x0i; 1189e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2 + 1] = wk2r * x0i + wk2i * x0r; 1190e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x1r - x3i; 1191e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x1i + x3r; 1192e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = wk1r * x0r - wk1i * x0i; 1193e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = wk1r * x0i + wk1i * x0r; 1194e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x1r + x3i; 1195e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x1i - x3r; 1196e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3] = wk3r * x0r - wk3i * x0i; 1197e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3 + 1] = wk3r * x0i + wk3i * x0r; 1198e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 1199e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk1r = w[k2 + 2]; 1200e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk1i = w[k2 + 3]; 1201e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk3r = wk1r - 2 * wk2r * wk1i; 1202e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wk3i = 2 * wk2r * wk1r - wk1i; 1203e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = k + m; j < l + (k + m); j += 2) { 1204e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j1 = j + l; 1205e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j2 = j1 + l; 1206e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent j3 = j2 + l; 1207e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = a[j] + a[j1]; 1208e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = a[j + 1] + a[j1 + 1]; 1209e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1r = a[j] - a[j1]; 1210e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x1i = a[j + 1] - a[j1 + 1]; 1211e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2r = a[j2] + a[j3]; 1212e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x2i = a[j2 + 1] + a[j3 + 1]; 1213e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3r = a[j2] - a[j3]; 1214e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x3i = a[j2 + 1] - a[j3 + 1]; 1215e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = x0r + x2r; 1216e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = x0i + x2i; 1217e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r -= x2r; 1218e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i -= x2i; 1219e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2] = -wk2i * x0r - wk2r * x0i; 1220e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j2 + 1] = -wk2i * x0i + wk2r * x0r; 1221e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x1r - x3i; 1222e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x1i + x3r; 1223e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1] = wk1r * x0r - wk1i * x0i; 1224e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j1 + 1] = wk1r * x0i + wk1i * x0r; 1225e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0r = x1r + x3i; 1226e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent x0i = x1i - x3r; 1227e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3] = wk3r * x0r - wk3i * x0i; 1228e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j3 + 1] = wk3r * x0i + wk3i * x0r; 1229e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 1230e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 1231e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 1232e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1233e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1234c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void rftfsub(int n, float *a, int nc, float *c) 1235e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 1236e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, k, kk, ks, m; 1237e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float wkr, wki, xr, xi, yr, yi; 1238e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1239e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = n >> 1; 1240e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ks = 2 * nc / m; 1241e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent kk = 0; 1242e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 2; j < m; j += 2) { 1243e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k = n - j; 1244e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent kk += ks; 1245e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wkr = 0.5f - c[nc - kk]; 1246e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wki = c[kk]; 1247e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j] - a[k]; 1248e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j + 1] + a[k + 1]; 1249e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = wkr * xr - wki * xi; 1250e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = wkr * xi + wki * xr; 1251e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] -= yr; 1252e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] -= yi; 1253e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] += yr; 1254e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k + 1] -= yi; 1255e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 1256e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 1257e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1258e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1259c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void rftbsub(int n, float *a, int nc, float *c) 1260e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 1261e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, k, kk, ks, m; 1262e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float wkr, wki, xr, xi, yr, yi; 1263e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1264e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[1] = -a[1]; 1265e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = n >> 1; 1266e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ks = 2 * nc / m; 1267e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent kk = 0; 1268e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 2; j < m; j += 2) { 1269e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k = n - j; 1270e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent kk += ks; 1271e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wkr = 0.5f - c[nc - kk]; 1272e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wki = c[kk]; 1273e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = a[j] - a[k]; 1274e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xi = a[j + 1] + a[k + 1]; 1275e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yr = wkr * xr + wki * xi; 1276e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent yi = wkr * xi - wki * xr; 1277e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] -= yr; 1278e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j + 1] = yi - a[j + 1]; 1279e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] += yr; 1280e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k + 1] = yi - a[k + 1]; 1281e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 1282e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[m + 1] = -a[m + 1]; 1283e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 1284e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1285c55a96383497a772a307b346368133960b02ad03Eric Laurent#if 0 // Not used. 1286c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void dctsub(int n, float *a, int nc, float *c) 1287e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 1288e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, k, kk, ks, m; 1289e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float wkr, wki, xr; 1290e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1291e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = n >> 1; 1292e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ks = nc / n; 1293e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent kk = 0; 1294e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 1; j < m; j++) { 1295e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k = n - j; 1296e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent kk += ks; 1297e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wkr = c[kk] - c[nc - kk]; 1298e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wki = c[kk] + c[nc - kk]; 1299e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = wki * a[j] - wkr * a[k]; 1300e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = wkr * a[j] + wki * a[k]; 1301e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] = xr; 1302e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 1303e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[m] *= c[0]; 1304e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 1305e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1306e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1307c55a96383497a772a307b346368133960b02ad03Eric Laurentstatic void dstsub(int n, float *a, int nc, float *c) 1308e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent{ 1309e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent int j, k, kk, ks, m; 1310e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent float wkr, wki, xr; 1311e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent 1312e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent m = n >> 1; 1313e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent ks = nc / n; 1314e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent kk = 0; 1315e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent for (j = 1; j < m; j++) { 1316e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent k = n - j; 1317e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent kk += ks; 1318e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wkr = c[kk] - c[nc - kk]; 1319e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent wki = c[kk] + c[nc - kk]; 1320e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent xr = wki * a[k] - wkr * a[j]; 1321e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[k] = wkr * a[k] + wki * a[j]; 1322e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[j] = xr; 1323e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent } 1324e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent a[m] *= c[0]; 1325e48d5845c8b35de2ab73ea055c18a61fa3a9f0beEric Laurent} 1326c55a96383497a772a307b346368133960b02ad03Eric Laurent#endif // Not used. 1327