Преобразование Фурье

#ifndef FFT_H
#define FFT_H
 
#include <complex>
#include <vector>
 
typedef std::complex<double> base;
 
class FFT
{
    static const int mod = 7340033;
    static const int root = 5;
    static const int root_1 = 4404020;
    static const int root_pw = 1 << 20;
public:
    static void fft(std::vector<base> & a, bool invert);
    static void multiply(const std::vector<int> & a, const std::vector<int> & b, std::vector<int> & res);
    static void fft(std::vector<int> & a, bool invert);
};
 
#endif // !FFT_H

#include "FFT.h"
 
#include <algorithm>
#include <cmath>
#define PI 3.1415926535897932384626433832795
 
using namespace std;
 
void FFT::fft(vector<base> & a, bool invert)
{
    int n = (int)a.size();
 
    for (int i = 1, j = 0; i < n; ++i)
    {
        int bit = n >> 1;
        for (; j >= bit; bit >>= 1)
            j -= bit;
        j += bit;
        if (i < j)
            swap(a[i], a[j]);
    }
 
    for (int len = 2; len <= n; len <<= 1)
    {
        double ang = 2 * PI / len * (invert ? -1 : 1);
        base wlen(cos(ang), sin(ang));
        for (int i = 0; i < n; i += len)
        {
            base w(1);
            for (int j = 0; j < len / 2; ++j)
            {
                base u = a[i + j], v = a[i + j + len / 2] * w;
                a[i + j] = u + v;
                a[i + j + len / 2] = u - v;
                w *= wlen;
            }
        }
    }
    if (invert)
        for (int i = 0; i < n; ++i)
            a[i] /= n;
}
 
void FFT::multiply(const vector<int> & a, const vector<int> & b, vector<int> & res) {
    vector<base> fa(a.begin(), a.end()), fb(b.begin(), b.end());
    size_t n = 1;
    while (n < max(a.size(), b.size()))
        n <<= 1;
    n <<= 1;
    fa.resize(n), fb.resize(n);
 
    fft(fa, false), fft(fb, false);
    for (size_t i = 0; i < n; ++i)
        fa[i] *= fb[i];
    fft(fa, true);
 
    res.resize(n);
    for (size_t i = 0; i < n; ++i)
        res[i] = int(fa[i].real() + 0.5);
}
 
void FFT::fft(vector<int> & a, bool invert) {
    int n = (int)a.size();
 
    for (int i = 1, j = 0; i < n; ++i) {
        int bit = n >> 1;
        for (; j >= bit; bit >>= 1)
            j -= bit;
        j += bit;
        if (i < j)
            swap(a[i], a[j]);
    }
 
    for (int len = 2; len <= n; len <<= 1) {
        int wlen = invert ? root_1 : root;
        for (int i = len; i < root_pw; i <<= 1)
            wlen = int(wlen * 1ll * wlen % mod);
        for (int i = 0; i < n; i += len) {
            int w = 1;
            for (int j = 0; j < len / 2; ++j) {
                int u = a[i + j], v = int(a[i + j + len / 2] * 1ll * w % mod);
                a[i + j] = u + v < mod ? u + v : u + v - mod;
                a[i + j + len / 2] = u - v >= 0 ? u - v : u - v + mod;
                w = int(w * 1ll * wlen % mod);
            }
        }
    }
    /*if (invert) {
        int nrev = reverse(n, mod);
        for (int i = 0; i < n; ++i)
            a[i] = int(a[i] * 1ll * nrev % mod);
    }*/
}

#include <iostream>
#include <vector>
 
#include "FFT.h"
 
using namespace std;
 
int main()
{
    vector< complex<double> > d;
    for (double i = 0; i < 1024.0; i += 0.5)
    {
        double temp = 5.35 * sin(2.0 * 3.1415926535897932384626433832795 * 5.0 * i / 1024.0) +
                      3.0 * cos(2.0 * 3.1415926535897932384626433832795 * 10.0 * i / 1024.0) +
                      7.0 * sin(2.0 * 3.1415926535897932384626433832795 * 12.0 * i / 1024.0 + 3.1415926535897932384626433832795 / 3.0);
        d.push_back((complex<double>)(temp));
    }
 
    vector< complex<double> > g(d);
 
    FFT::fft(g, false);
 
    double dt = 1.0 / (double)g.size();
    for (double i = 0; i < (double)g.size() / 2; i += 1)
    {
        double fi = double(i) / dt / (double)g.size();
        double ai = 2.0 / (double)g.size() / (double)g.size() / dt * sqrt(g[i].real() * g[i].real() + g[i].imag() * g[i].imag());
        if (ai > 0.001)
            cout << fi << "    " << ai << endl;
    }
    cout << endl << endl;
    cout << d.size() << endl << endl;
    system("PAUSE");
 
    FFT::fft(d, true);
 
    for (int i = 0; i < int(d.size()); i++)
    {
        if (d[i].real() - g[i].real() > 0.001)
            cout << "ERROR! " << abs(d[i].real() - g[i].real()) << endl;
    }
 
    system("PAUSE");
    return 0;
}

for (double i = 0; i < 1024.0; i += 0.5)
    {
        double temp = 5.35 * sin(2.0 * 3.1415926535897932384626433832795 * 5.2 * i / 1024.0) +
                      3.0 * cos(2.0 * 3.1415926535897932384626433832795 * 10.0 * i / 1024.0) +
                      7.0 * sin(2.0 * 3.1415926535897932384626433832795 * 12.0 * i / 1024.0 + 3.1415926535897932384626433832795 / 3.0);
        d.push_back((complex<double>)(temp));
    }

for (double i = 0; i < 1024.0; i += 0.1)

for (double i = 0; i < 1024.0; i += 0.5)
    {
        double temp = 5.35 * sin(2.0 * 3.1415926535897932384626433832795 * 5.0 * i / 1024.0) +
                      3.0 * cos(2.0 * 3.1415926535897932384626433832795 * 10.0 * i / 1024.0) +
                      7.0 * sin(2.0 * 3.1415926535897932384626433832795 * 12.0 * i / 1024.0 + 3.1415926535897932384626433832795 / 3.0);
        d.push_back((complex<double>)(temp));
    }

for (double i = 0; i < 1024.0; i += 0.5)
    {
        double temp = 5.35 * sin(2.0 * 3.1415926535897932384626433832795 * 5.5 * i / 1024.0) +
                      3.0 * cos(2.0 * 3.1415926535897932384626433832795 * 10.0 * i / 1024.0) +
                      7.0 * sin(2.0 * 3.1415926535897932384626433832795 * 12.0 * i / 1024.0 + 3.1415926535897932384626433832795 / 3.0);
        d.push_back((complex<double>)(temp));
    }

for (double i = 0; i < 2048.0; i += 0.5)
    {
        double temp = 5.35 * sin(2.0 * 3.1415926535897932384626433832795 * 5.5 * i / 2048.0) +
                      3.0 * cos(2.0 * 3.1415926535897932384626433832795 * 10.0 * i / 2048.0) +
                      7.0 * sin(2.0 * 3.1415926535897932384626433832795 * 12.0 * i / 2048.0 + 3.1415926535897932384626433832795 / 3.0);
        d.push_back((complex<double>)(temp));
    }