123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354 |
- struct cpx{
- cpx(){}
- cpx(double aa):a(aa) {}
- cpx(double aa, double bb):a(aa),b(bb) {}
- double a;
- double b;
- double modsq(void) const {
- return a * a + b * b;
- }
- cpx bar(void) const {
- return cpx(a, -b);
- }
- };
- cpx operator +(cpx a, cpx b) {
- return cpx(a.a + b.a, a.b + b.b);
- }
- cpx operator *(cpx a, cpx b) {
- return cpx(a.a * b.a - a.b * b.b, a.a * b.b + a.b * b.a);
- }
- cpx operator /(cpx a, cpx b) {
- cpx r = a * b.bar();
- return cpx(r.a / b.modsq(), r.b / b.modsq());
- }
- cpx EXP(double theta) {
- return cpx(cos(theta),sin(theta));
- }
- const double two_pi = 4 * acos(0);
- // in: input array
- // out: output array
- // step: {SET TO 1} (used internally)
- // size: length of the input/output {MUST BE A POWER OF 2} (use nextPowerOf2)
- // dir: either plus or minus one (direction of the FFT)
- // RESULT: out[k] = \sum_{j=0}^{size - 1} in[j] * exp(dir * 2pi * i * j * k / size)
- void FFT(cpx *in, cpx *out, int step, int size, int dir) {
- if(size < 1) return;
- if(size == 1) {
- out[0] = in[0];
- return;
- }
- FFT(in, out, step * 2, size / 2, dir);
- FFT(in + step, out + size / 2, step * 2, size / 2, dir);
- for(int i = 0 ; i < size / 2 ; i++) {
- cpx even = out[i];
- cpx odd = out[i + size / 2];
- out[i] = even + EXP(dir * two_pi * i / size) * odd;
- out[i + size / 2] = even + EXP(dir * two_pi * (i + size / 2) / size) * odd;
- }
- }
- unsigned int nextPowerOf2(unsigned int n) {
- unsigned int p = 1;
- if (n && !(n & (n - 1))) return n;
- while (p < n) p <<= 1;
- return p;
- }
|