Koichi Kurahashi / Mbed 2 deprecated FRDM_MAG3110

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との接続は以下の通りです：

AQM1248KL46
Vcc3.3v
CSD10
RESETD9
RSD8
SCLKD13
SDID11

fft4g.cpp@1:ad135c286d4d, 2015-05-04 (annotated)

Committer:
TareObjects
Date:
Mon May 04 08:53:50 2015 +0000
Revision:
1:ad135c286d4d
Parent:
0:47be4d9de4b9 
Faster; ; - Less Sampling datas (256->128); - Wait for status, and fetch data; // I tried to use INT1, but I could not get interrupt....

Who changed what in which revision?

UserRevisionLine numberNew 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