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