C++ FFT::SetFlag方法代码示例

void test_scalar_generic(int nfft)
{
    typedef typename FFT<T>::Complex Complex;
    typedef typename FFT<T>::Scalar Scalar;
    typedef typename VectorType<Container, Scalar>::type ScalarVector;
    typedef typename VectorType<Container, Complex>::type ComplexVector;

    FFT<T>        fft;
    ScalarVector  tbuf(nfft);
    ComplexVector freqBuf;

    for (int k = 0; k < nfft; ++k)
        tbuf[k] = (T)(rand() / (double)RAND_MAX - .5);

    // make sure it DOESN'T give the right full spectrum answer
    // if we've asked for half-spectrum
    fft.SetFlag(fft.HalfSpectrum);
    fft.fwd(freqBuf, tbuf);
    VERIFY((size_t)freqBuf.size() == (size_t)((nfft >> 1) + 1));
    VERIFY(fft_rmse(freqBuf, tbuf) < test_precision<T>());  // gross check

    fft.ClearFlag(fft.HalfSpectrum);
    fft.fwd(freqBuf, tbuf);
    VERIFY((size_t)freqBuf.size() == (size_t)nfft);
    VERIFY(fft_rmse(freqBuf, tbuf) < test_precision<T>());  // gross check

    if (nfft & 1)
        return; // odd FFTs get the wrong size inverse FFT

    ScalarVector tbuf2;
    fft.inv(tbuf2, freqBuf);
    VERIFY(dif_rmse(tbuf, tbuf2) < test_precision<T>());  // gross check


    // verify that the Unscaled flag takes effect
    ScalarVector tbuf3;
    fft.SetFlag(fft.Unscaled);

    fft.inv(tbuf3, freqBuf);

    for (int k = 0; k < nfft; ++k)
        tbuf3[k] *= T(1. / nfft);


    // for (size_t i=0;i<(size_t) tbuf.size();++i)
    //    cout << "freqBuf=" << freqBuf[i] << " in2=" << tbuf3[i] << " -  in=" << tbuf[i] << " => " << (tbuf3[i] - tbuf[i] ) <<  endl;

    VERIFY(dif_rmse(tbuf, tbuf3) < test_precision<T>());  // gross check

    // verify that ClearFlag works
    fft.ClearFlag(fft.Unscaled);
    fft.inv(tbuf2, freqBuf);
    VERIFY(dif_rmse(tbuf, tbuf2) < test_precision<T>());  // gross check
}

void bench(int nfft,bool fwd,bool unscaled=false, bool halfspec=false)
{
    typedef typename NumTraits<T>::Real Scalar;
    typedef typename std::complex<Scalar> Complex;
    int nits = NDATA/nfft;
    vector<T> inbuf(nfft);
    vector<Complex > outbuf(nfft);
    FFT< Scalar > fft;

    if (unscaled) {
        fft.SetFlag(fft.Unscaled);
        cout << "unscaled ";
    }
    if (halfspec) {
        fft.SetFlag(fft.HalfSpectrum);
        cout << "halfspec ";
    }


    std::fill(inbuf.begin(),inbuf.end(),0);
    fft.fwd( outbuf , inbuf);

    BenchTimer timer;
    timer.reset();
    for (int k=0;k<8;++k) {
        timer.start();
        if (fwd)
            for(int i = 0; i < nits; i++)
                fft.fwd( outbuf , inbuf);
        else
            for(int i = 0; i < nits; i++)
                fft.inv(inbuf,outbuf);
        timer.stop();
    }

    cout << nameof<Scalar>() << " ";
    double mflops = 5.*nfft*log2((double)nfft) / (1e6 * timer.value() / (double)nits );
    if ( NumTraits<T>::IsComplex ) {
        cout << "complex";
    }else{
        cout << "real   ";
        mflops /= 2;
    }


    if (fwd)
        cout << " fwd";
    else
        cout << " inv";

    cout << " NFFT=" << nfft << "  " << (double(1e-6*nfft*nits)/timer.value()) << " MS/s  " << mflops << "MFLOPS\n";
}

void test_complex_generic(int nfft)
{
    typedef typename FFT<T>::Complex Complex;
    typedef typename VectorType<Container,Complex>::type ComplexVector;

    FFT<T> fft;

    ComplexVector inbuf(nfft);
    ComplexVector outbuf;
    ComplexVector buf3;
    for (int k=0;k<nfft;++k)
        inbuf[k]= Complex( (T)(rand()/(double)RAND_MAX - .5), (T)(rand()/(double)RAND_MAX - .5) );
    fft.fwd( outbuf , inbuf);

    VERIFY( fft_rmse(outbuf,inbuf) < test_precision<T>()  );// gross check
    fft.inv( buf3 , outbuf);

    VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>()  );// gross check

    // verify that the Unscaled flag takes effect
    ComplexVector buf4;
    fft.SetFlag(fft.Unscaled);
    fft.inv( buf4 , outbuf);
    for (int k=0;k<nfft;++k)
        buf4[k] *= T(1./nfft);
    VERIFY( dif_rmse(inbuf,buf4) < test_precision<T>()  );// gross check

    // verify that ClearFlag works
    fft.ClearFlag(fft.Unscaled);
    fft.inv( buf3 , outbuf);
    VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>()  );// gross check
}

void test_return_by_value(int len)
{
    VectorXf in;
    VectorXf in1;
    in.setRandom( len );
    VectorXcf out1,out2;
    FFT<float> fft;

    fft.SetFlag(fft.HalfSpectrum );

    fft.fwd(out1,in);
    out2 = fft.fwd(in);
    VERIFY( (out1-out2).norm() < test_precision<float>() );
    in1 = fft.inv(out1);
    VERIFY( (in1-in).norm() < test_precision<float>() );
}

void UnbiasedSquaredPhaseLagIndex::compute(ConnectivitySettings::IntermediateTrialData& inputData,
                                           QVector<QPair<int,MatrixXcd> >& vecPairCsdSum,
                                           QVector<QPair<int,MatrixXd> >& vecPairCsdImagSignSum,
                                           QMutex& mutex,
                                           int iNRows,
                                           int iNFreqs,
                                           int iNfft,
                                           const QPair<MatrixXd, VectorXd>& tapers)
{
    if(inputData.vecPairCsdImagSign.size() == iNRows) {
        //qDebug() << "UnbiasedSquaredPhaseLagIndex::compute - vecPairCsdImagSign was already computed for this trial.";
        return;
    }

    inputData.vecPairCsdImagSign.clear();

    int i,j;

    // Calculate tapered spectra if not available already
    // This code was copied and changed modified Utils/Spectra since we do not want to call the function due to time loss.
    if(inputData.vecTapSpectra.size() != iNRows) {
        inputData.vecTapSpectra.clear();

        RowVectorXd vecInputFFT, rowData;
        RowVectorXcd vecTmpFreq;

        MatrixXcd matTapSpectrum(tapers.first.rows(), iNFreqs);

        QVector<Eigen::MatrixXcd> vecTapSpectra;

        FFT<double> fft;
        fft.SetFlag(fft.HalfSpectrum);

        for (i = 0; i < iNRows; ++i) {
            // Substract mean
            rowData.array() = inputData.matData.row(i).array() - inputData.matData.row(i).mean();

            // Calculate tapered spectra if not available already
            for(j = 0; j < tapers.first.rows(); j++) {
                vecInputFFT = rowData.cwiseProduct(tapers.first.row(j));
                // FFT for freq domain returning the half spectrum and multiply taper weights
                fft.fwd(vecTmpFreq, vecInputFFT, iNfft);
                matTapSpectrum.row(j) = vecTmpFreq * tapers.second(j);
            }

            inputData.vecTapSpectra.append(matTapSpectrum);
        }
    }

    // Compute CSD
    if(inputData.vecPairCsd.isEmpty()) {
        double denomCSD = sqrt(tapers.second.cwiseAbs2().sum()) * sqrt(tapers.second.cwiseAbs2().sum()) / 2.0;

        bool bNfftEven = false;
        if (iNfft % 2 == 0){
            bNfftEven = true;
        }

        MatrixXcd matCsd = MatrixXcd(iNRows, iNFreqs);

        for (i = 0; i < iNRows; ++i) {
            for (j = i; j < iNRows; ++j) {
                // Compute CSD (average over tapers if necessary)
                matCsd.row(j) = inputData.vecTapSpectra.at(i).cwiseProduct(inputData.vecTapSpectra.at(j).conjugate()).colwise().sum() / denomCSD;

                // Divide first and last element by 2 due to half spectrum
                matCsd.row(j)(0) /= 2.0;
                if(bNfftEven) {
                    matCsd.row(j).tail(1) /= 2.0;
                }
            }

            inputData.vecPairCsd.append(QPair<int,MatrixXcd>(i,matCsd));
            inputData.vecPairCsdImagSign.append(QPair<int,MatrixXd>(i,matCsd.imag().cwiseSign()));
        }

        mutex.lock();

        if(vecPairCsdSum.isEmpty()) {
            vecPairCsdSum = inputData.vecPairCsd;
            vecPairCsdImagSignSum = inputData.vecPairCsdImagSign;
        } else {
            for (int j = 0; j < vecPairCsdSum.size(); ++j) {
                vecPairCsdSum[j].second += inputData.vecPairCsd.at(j).second;
                vecPairCsdImagSignSum[j].second += inputData.vecPairCsdImagSign.at(j).second;
            }
        }

        mutex.unlock();
    } else {
        if(inputData.vecPairCsdImagSign.isEmpty()) {
            for (i = 0; i < inputData.vecPairCsd.size(); ++i) {
                inputData.vecPairCsdImagSign.append(QPair<int,MatrixXd>(i,inputData.vecPairCsd.at(i).second.imag().cwiseSign()));
            }

            mutex.lock();

            if(vecPairCsdImagSignSum.isEmpty()) {
                vecPairCsdImagSignSum = inputData.vecPairCsdImagSign;
            } else {
//.........这里部分代码省略.........