1 /*
2 * FFT.c
3 * PRICE
4 *
5 * Created by Riccardo Mottola on Sat Sep 13 2003.
6 * Copyright (c) 2002-2005 Carduus. All rights reserved.
7 *
8 */
9
10 // This application is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
11 // This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12
13 #include "FFT.h"
14
15 double *cosinus;
16 double *sinus;
17 unsigned int *bitRevIndex;
18
initTrigonometrics(int num,unsigned int bitNumber)19 int initTrigonometrics(int num, unsigned int bitNumber)
20 /* we init here the trigonometric arrays
21 but we also init the bit reerse array */
22 {
23 int i;
24 double alpha;
25
26 alpha = (2 * PI) / num;
27 for (i = 0; i < num; i++)
28 {
29 cosinus[i] = cos(-alpha * i);
30 sinus[i] = sin(-alpha * i);
31 }
32
33 for (i = 0; i < num; i++)
34 bitRevIndex[i] = bitrev(i, bitNumber);
35
36 return 0;
37 }
38
binaryLog(unsigned int a)39 unsigned int binaryLog(unsigned int a)
40 /* returns the smallest 2 power that contains the a */
41 {
42 unsigned int i;
43 unsigned int j;
44
45 i = 1;
46 j = 1;
47 while (i < a)
48 {
49 i = i << 1;
50 j++;
51 }
52 return j - 1;
53 }
54
binpow(unsigned int n)55 unsigned int binpow(unsigned int n)
56 {
57 return 1 << n;
58 }
59
bitrev(unsigned int i,unsigned int len)60 unsigned int bitrev(unsigned int i, unsigned int len)
61 /* returns the bit-reverse of i given the bit number to work on */
62 {
63 unsigned int r;
64 unsigned int j;
65
66 r = 0;
67 for (j = 0; j < len-1; j++) /* we don't have to shift for the last position */
68 {
69 r = r | (i & 1);
70 r = r << 1;
71 i = i >> 1;
72 }
73 r = r | (i & 1);
74 return r;
75 }
76
fft(int num,unsigned int bitNumber,double Ar[],double Ai[],double yr[],double yi[])77 int fft(int num, unsigned int bitNumber, double Ar[], double Ai[], double yr[], double yi[])
78 {
79 unsigned int i;
80 unsigned int s;
81 unsigned int m;
82 unsigned int j;
83 unsigned int k;
84 register double omega_r, omega_i; /* current unity root */
85 register unsigned int gamma;
86 register double t_r, t_i;
87 register double u_r, u_i;
88
89 for (i = 0; i < num; i++) /* exponents are 0 to n-1 for n bits */
90 {
91 yr[bitRevIndex[i]] = Ar[i];
92 yi[bitRevIndex[i]] = Ai[i];
93 }
94 for (s = 1; s <= bitNumber; s++)
95 {
96 m = binpow(s);
97 for (j = 0; j < m/2; j++)
98 {
99 gamma = ((num / m) * j) % num;
100 omega_r = cosinus[gamma];
101 omega_i = sinus[gamma];
102 for (k = j; k < num; k += m)
103 {
104 t_r = omega_r * yr[k + m/2] - omega_i * yi[k + m/2];
105 t_i = omega_r * yi[k + m/2] + omega_i * yr[k + m/2];
106 u_r = yr[k];
107 u_i = yi[k];
108 yr[k] = u_r + t_r;
109 yi[k] = u_i + t_i;
110 yr[k + m/2] = u_r - t_r;
111 yi[k + m/2] = u_i - t_i;
112 }
113 }
114 }
115 return 0;
116 }
117
ifft(int num,unsigned int bitNumber,double Ar[],double Ai[],double yr[],double yi[])118 int ifft(int num, unsigned int bitNumber, double Ar[], double Ai[], double yr[], double yi[])
119 {
120 unsigned int i;
121 unsigned int s;
122 unsigned int m;
123 unsigned int j;
124 unsigned int k;
125 register double omega_r, omega_i; /* current unity root */
126 register unsigned int gamma;
127 register double t_r, t_i;
128 register double u_r, u_i;
129
130
131 for (i = 0; i < num; i++) /* exponents are 0 to n-1 for n bits */
132 {
133 yr[bitRevIndex[i]] = Ar[i];
134 yi[bitRevIndex[i]] = Ai[i];
135 }
136 for (s = 1; s <= bitNumber; s++)
137 {
138 m = binpow(s);
139 for (j = 0; j < m/2; j++)
140 {
141 /* we follow the unity roots in the opposite direction */
142 /* but since sin(-a) = -sin(a) but cos(-a)=cos(a) just a sign is needed */
143 gamma = ((num / m) * j) % num;
144 omega_r = cosinus[gamma];
145 omega_i = -sinus[gamma];
146 for (k = j; k < num; k += m)
147 {
148 t_r = omega_r * yr[k + m/2] - omega_i * yi[k + m/2];
149 t_i = omega_r * yi[k + m/2] + omega_i * yr[k + m/2];
150 u_r = yr[k];
151 u_i = yi[k];
152 yr[k] = u_r + t_r;
153 yi[k] = u_i + t_i;
154 yr[k + m/2] = u_r - t_r;
155 yi[k + m/2] = u_i - t_i;
156 }
157 }
158 }
159 for (i = 0; i < num; i++)
160 {
161 yr[i] = yr[i]/num;
162 yi[i] = yi[i]/num;
163 }
164 return 0;
165 }
166