a shifted scaled fourier

gistfile1.txt

    // Do a shifted/scaled fourier transform of polarization in each direction.
    mat Fourier(const mat& pol, double zoom = 1.0/20.0, int axis=0, bool dering=true) const
    {
        cout << "Computing Fourier Transform ..." << endl;
        mat tore;
        double dt = pol(3,0)-pol(2,0);
        {
            int ndata = (pol.n_rows%2==0)? pol.n_rows : pol.n_rows-1; // Keep an even amount of data
            for (int r=1; r<ndata; ++r)
                if (pol(r,0)==0.0)
                    ndata =	(r%2==0)?r:r-1;  // Ignore trailing zeros.
            mat data(ndata,1);
            data.col(0) = pol.rows(0,ndata-1).col(1+axis);
            data = data - mean(data); // subtract any zero frequency component.

            if (dering)
            {
                double gam = log(0.00001/(ndata*ndata));
                mat de = linspace<mat>(0.0,ndata,ndata);
                de = exp(-1.0*de%de*gam);
                data = data%de;
            }

            cx_mat f(ndata,1); f.zeros();
            mat fr(ndata/2,1), fi(ndata/2,1); fr.zeros(); fi.zeros();

#pragma omp parallel for schedule(guided)
            for (int r=0; r<ndata; ++r)
                f(r,0) = exp(std::complex<double>(0.0,2.0*M_PI)*zoom*(r-1.0)/ndata);

#pragma omp parallel for schedule(guided)
            for (int i=0; i<ndata/2; ++i) // I only even generate the positive frequencies.
            {
                cx_mat ftmp = pow(f,std::complex<double>((double)i,0.0));
                fr(i,0) = real(accu(ftmp%data));
                fi(i,0) = imag(accu(ftmp%data));
            }

            tore.resize(ndata/2,3);
            for (int i=0; i<ndata/2; ++i)
                tore(i,0) = M_PI*(27.2113*DipolesEvery/dt)*zoom*i/ndata;
            tore.col(1) = -1.0*fi.col(0);
            tore.col(2) = fr.col(0);
        }
        return tore;
    }