C++ FFT类代码示例

void Test::test_PolyMult_works() {
  srand(time(0));

  int n = 4;

  int randMin = -100, randMax = 100;

  Vector<Complex> P(n), Q(n);

  for(int i = 0; i < n; ++i) {

    double real = rand()%(randMax-randMin)+randMin;
    double imaginary = rand()%(randMax-randMin)+randMin;

    P[i] = Complex(real, imaginary);
    Q[i] = Complex(real, imaginary);
  }

  FFT fft;

  fft.preComputeOmega(n);
  Vector<Complex> resultsR = fft.polyMultR(P, Q, n);
  fft.preComputeRBS(n);
  Vector<Complex> resultsD = fft.polyMultD(P, Q, n);

  for(int i = 0; i < n; ++i) {
    if(resultsD[i] != resultsR[i])
      std::cout << "*** Failed ***" << std::endl;

    std::cout << resultsR[i] << " " << resultsD[i] << std::endl;
  }
}

void ofApp::audioIn(float * input, int bufferSize, int nChannels){

    if( !processAudio ) return;

    for (int i=0; i<bufferSize; i++) {
        trainingSample[i][0] = input[i];
    }

    if( record ){
        trainingData.addSample( trainingClassLabel, trainingSample );
    }

    if( pipeline.getTrained() ){

        //Run the prediction using the matrix of audio data
        pipeline.predict( trainingSample );

        //Update the FFT plot
        FFT *fft = pipeline.getFeatureExtractionModule< FFT >( 0 );
        if( fft ){
            vector< FastFourierTransform > &results =  fft->getFFTResultsPtr();
            magnitudePlot.setData( results[0].getMagnitudeData() );
        }

        //Update the likelihood plot
        classLikelihoodsPlot.update( pipeline.getClassLikelihoods() );
    }
}

bool test3() {
    FFT fft;

    vector<complex<double> > in, out, expected, temp;
    in.push_back(complex<double>(0.0,0.0));
    in.push_back(complex<double>(1.0,0.0));
    in.push_back(complex<double>(0.0,0.0));
    in.push_back(complex<double>(0.0,0.0));
    in.push_back(complex<double>(0.0,0.0));

    expected.push_back(complex<double>( 0.125,  0.0));
    expected.push_back(complex<double>( 0.088, -0.088));
    expected.push_back(complex<double>( 0.000, -0.125));
    expected.push_back(complex<double>(-0.088, -0.088));
    expected.push_back(complex<double>(-0.125,  0.0));
    expected.push_back(complex<double>(-0.088,  0.088));
    expected.push_back(complex<double>( 0.000,  0.125));
    expected.push_back(complex<double>( 0.088,  0.088));

    fft.FT(in,out);
    
    if (!checkEqual(expected,out)) {
	printError(expected,out);
	cout << __FILE__ << ": " << __LINE__ << endl;
	return false;
    } 
    
    fft.inverseFT(out, temp);

    return true;

}

//--------------------------------------------------------------
void ofApp::setup(){
    
    ofSetFrameRate(60);
    
    //Setup the FFT
    FFT fft;
    fft.init(FFT_WINDOW_SIZE,FFT_HOP_SIZE,1,FFT::RECTANGULAR_WINDOW,true,false,DATA_TYPE_MATRIX);

    //Setup the classifier
    RandomForests forest;
    forest.setForestSize( 10 );
    forest.setNumRandomSplits( 100 );
    forest.setMaxDepth( 10 );
    forest.setMinNumSamplesPerNode( 10 );

    //Add the feature extraction and classifier to the pipeline
    pipeline.addFeatureExtractionModule( fft );
    pipeline.setClassifier( forest );

    trainingClassLabel = 1;
    record = false;
    processAudio = true;
    trainingData.setNumDimensions( 1 ); //We are only going to use the data from one microphone channel, so the dimensions are 1
    trainingSample.resize( AUDIO_BUFFER_SIZE, 1 ); //We will set the training matrix to match the audio buffer size

    //Setup the audio card
    ofSoundStreamSetup(2, 1, this, AUDIO_SAMPLE_RATE, AUDIO_BUFFER_SIZE, 4);
}

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 fft_Y(MultiArray<complex<float>, 3> &a) {
    FFT<float> fft;
    int ny = a.dims()[1];
    for (int z = 0; z < a.dims()[2]; z++) {
        for (int x = 0; x < a.dims()[0]; x++) {
            VectorXcf fft_in = a.slice<1>({x,0,z},{0,ny,0}).asArray();
            VectorXcf fft_out(ny);
            fft.fwd(fft_out, fft_in);
            a.slice<1>({x,0,z},{0,ny,0}).asArray() = fft_out;
        }
    }
}

void fft_Z(MultiArray<complex<float>, 3> &a) {
    FFT<float> fft;
    int nz = a.dims()[2];
    for (int x = 0; x < a.dims()[0]; x++) {
        for (int y = 0; y < a.dims()[1]; y++) {
            VectorXcf fft_in = a.slice<1>({x,y,0},{0,0,nz}).asArray();
            VectorXcf fft_out(nz);
            fft.fwd(fft_out, fft_in);
            a.slice<1>({x,y,0},{0,0,nz}).asArray() = fft_out;
        }
    }
}

void fft_X(MultiArray<complex<float>, 3> &a) {
    FFT<float> fft;
    int nx = a.dims()[0];
    for (int z = 0; z < a.dims()[2]; z++) {
        for (int y = 0; y < a.dims()[1]; y++) {
            VectorXcf fft_in = a.slice<1>({0,y,z},{nx,0,0}).asArray();
            VectorXcf fft_out(nx);
            fft.fwd(fft_out, fft_in);
            a.slice<1>({0,y,z},{nx,0,0}).asArray() = fft_out;
        }
    }
}

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 fwd_inv(size_t nfft)
{
    typedef typename NumTraits<T_freq>::Real Scalar;
    vector<T_time> timebuf(nfft);
    RandomFill(timebuf);

    vector<T_freq> freqbuf;
    static FFT<Scalar> fft;
    fft.fwd(freqbuf,timebuf);

    vector<T_time> timebuf2;
    fft.inv(timebuf2,freqbuf);

    long double rmse = mag2(timebuf - timebuf2) / mag2(timebuf);
    cout << "roundtrip rmse: " << rmse << endl;
}

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
}

int main()
{
	int n;
	scanf("%d", &n);
	while(n--)
	{
		string a, b;
		cin>>a>>b;
		if(a == "0" || b == "0")
		{
			cout<<"0\n";
			continue;
		}
		vector<int> ai;
		vector<int> bi;
		vector<int> res;
		for(int i = a.size()-1; i >= 0; --i)
		{
			ai.pb((a[i]-'0'));
		}
		for(int i = b.size()-1; i >= 0; --i)
		{
			bi.pb((b[i]-'0'));
		}
		FFT f;
		f.multiply(ai, bi, res);
		int carry=0;
		for(size_t i = 0; i < res.size(); i++)
		{
			res[i] += carry;
			carry = res[i]/10;
			res[i] %= 10;
		}
		int i = res.size()-1;
		while(res[i] == 0)
			i--;
		while(i >= 0)
		{
			cout<<res[i];
			--i;
		}
		cout<<endl;
	}
	return 0;
}

void Test::test_PolyMult_recursive() {
  int n = 2;
  Vector<Complex> P(n),Q(n);

  P[0] = Complex(2,0);
  P[1] = Complex(2,0);

  Q[0] = Complex(2,0);
  Q[1] = Complex(2,0);

  FFT fft;

  fft.preComputeOmega(n);
  Vector<Complex> results = fft.polyMultR(P, Q, n);

  for(Complex c : results) {
    std::cout << c << std::endl;
  }
}

vector<vector<double> > compute(vector<double> sound)
{
    FFT fft;
    vector<vector<double> > res;
    vector<complex<double> > tmp;
    while(sound.size()%SMPLS_PR_BLCK)
        sound.push_back(0);
    for(int i=0;i<sound.size();i+=SMPLS_PR_BLCK)
    {
        vector<complex<double> > t;
        for(int j=0;j<SMPLS_PR_BLCK;j++)
            t.push_back(complex<double>(sound[i+j],0));
        tmp = fft.compute(t);
        vector<double> tmp2;
        for(int i=0;i<tmp.size();i++)
            tmp2.push_back(tmp[i].real());
        res.push_back(transform(tmp2));
    }
    return res;
}

int main(int argc, char*argv[])
{
        // create FFT class implementation
        // each implementation has its own plan map.
        FFT<double> fft;
        // create some data
        size_t N = 1024 * 1024;
        std::vector<complex<double>,fftalloc<complex<double> > > data(N);
        std::vector<complex<double>,fftalloc<complex<double> > > dataFourier(N);
        std::vector<complex<double>,fftalloc<complex<double> > > dataCalc(N);
        std::cout << "1) --- original data ---" << endl;
        for(std::vector<complex<double> >::size_type i = 0; i != data.size(); i++) {
            data[i] = polar(sin(double(i)/N * M_PI * 2), 0.0);
//            std::cout << i << data[i] << endl;
        }

        // fft
        std::cout << "2) --- fft data ---" << endl;
        fft.fwd(dataFourier, data);
        std::cout << "speed / ms: " << fft.speed()/1000 << endl;
//        for(std::vector<complex<double> >::size_type i = 0; i != dataFourier.size(); i++) {
//            std::cout << i << dataFourier[i] << endl;
//        }

        // ifft
        std::cout << "3) --- ifft data ---" << endl;
        fft.inv(dataCalc, dataFourier);
        std::cout << "speed / ms: " << fft.speed()/1000 << endl;
//        for(std::vector<complex<double> >::size_type i = 0; i != dataCalc.size(); i++) {
//            std::cout << i << dataCalc[i] << endl;
//        }
//        std::cout << "4) --- comparison data ---" << endl;
//        for(std::vector<complex<double> >::size_type i = 0; i != dataCalc.size(); i++) {
//            std::cout << i << data[i] - dataCalc[i] << endl;
//        }

    std::getchar();
    return 0;
}