example4.cpp

利用这个模板可以分析基因表达数据

    
    ofstream fout;
    fout.open("fft1D.out");
    // check that file can be opened
    if (!fout.good())
    {
        cerr << "\n\nError creating fft1D.out" << endl;
        return;
    }
    
    // Perform a shift so that the origin lies in the centre.
    unsigned int j;
    for (unsigned int iiii=0; iiii<dataSize; iiii++)
    {
    	if (iiii < dataSize/2)
    	    j = dataSize/2 + iiii;
    	else
    	    j = iiii - dataSize/2;
    	fout << powerSpectrum_sum[j] << '\n';
    }    
    
    fout.close();       
}

//..............................................................................

void Series_2Dtest()
{

    // Create containers to contain test data
    //
    
    typedef vector< complex_type > Container;
    unsigned int dim = 64, nrows = dim, ncols = dim, dataSize = nrows*ncols;
    Container sin_freq_1(dataSize),
    	      sin_freq_2(dataSize),
    	      sin_freq_3(dataSize),
    	      sin_freq_sum(dataSize);
 
    // Put sin series with 4 different frequencies into containers
    //
    double rowstep, colstep, stepSize = 2.0 * M_PI / dim, start=-M_PI;
    unsigned int row,col,iarr;
    for (row=0, iarr=0, rowstep=start; row<nrows; row++, rowstep+=stepSize)
    	for (col=0, colstep=start; col<ncols; col++, colstep+=stepSize, iarr++)
    	{
    	      sin_freq_1[iarr] = complex_type(sin(1.0*rowstep)*sin(1.0*colstep));
    	      sin_freq_2[iarr] = complex_type(sin(2.0*rowstep)*sin(2.0*colstep));
    	      sin_freq_3[iarr] = complex_type(sin(3.0*rowstep)*sin(3.0*colstep));
    	    sin_freq_sum[iarr] = sin_freq_1[iarr] + 
    	                         sin_freq_2[iarr] +
    	                         sin_freq_3[iarr];
    	}
    
    // Output input data to a file
    //    
    ofstream fin;
    fin.open("fft2D.in");

    // check that file can be opened
    if (!fin)
    {
        cerr << "\n\nError creating fft2D.in" << endl;
        return;
    }
    for (row=0, iarr=0; row<nrows; row++)
    {
    	for (col=0; col<ncols; col++, iarr++)
    	    fin << real(sin_freq_sum[iarr]) << '\n';
    	fin << '\n';
    }
    
    // Create containers to contain Fourier transforms of the input data
    //
    Container output_1(dataSize),
    	      output_2(dataSize),
    	      output_3(dataSize),
    	      output_sum(dataSize);

    // Do the FFT calculations
    //    	      

    fourier_transform(sin_freq_1.begin(),output_1.begin(), nrows, ncols);
    fourier_transform(sin_freq_2.begin(),output_2.begin(), nrows, ncols);
    fourier_transform(sin_freq_3.begin(),output_3.begin(), nrows, ncols);
    fourier_transform(sin_freq_sum.begin(),output_sum.begin(), nrows, ncols);
   

    // Calculate the power spectrum of the resulting complex numbers to 
    // remove the imaginary part. Compare the theoretical result with the
    // calculated result
    //
    vector<double> powerSpectrum_sum(dataSize);
    double theoretical_result, diff, maxdiff = 0.0;

    for (unsigned int i=0; i<dataSize; i++)
    {
    	theoretical_result  = norm(output_1[i]+output_2[i]+output_3[i]);
    	                      
    	powerSpectrum_sum[i] = norm(output_sum[i]);
    	diff = fabs(theoretical_result - powerSpectrum_sum[i]);
    	if (diff > maxdiff) maxdiff = diff;
    }
    
    // Output the result to a file
    //
    cout << "2D Sin series test: max difference from ideal result = " << maxdiff << endl;
    
    ofstream fout;
    fout.open("fft2D.out");
    // check that file can be opened
    if (!fout.good())
    {
        cerr << "\n\nError creating fft2D.out" << endl;
        return;
    }
    
    // Perform a shift so that the origin lies in the centre.
    unsigned int x,y;
    for (row=0; row<nrows; row++)
    {
    	if (row < nrows/2)
    	    x = nrows/2 + row;
    	else
    	    x = row - nrows/2;
    	for (col=0; col<ncols; col++)
    	{
    	    if (col < ncols/2) 
    	    	y = ncols/2 + col;
    	    else
    	    	y = col - ncols/2;
    	    iarr = ncols*x + y; 
    	    fout << powerSpectrum_sum[iarr] << '\n';
    	}
    	fout << '\n';
    }    
}

//..............................................................................

template <class InputIterator>
void ArrayToFile3D(char *fileName, InputIterator begin, InputIterator end, 
                  unsigned int nrows,
                  unsigned int ncols,
                  unsigned int nslabs,
                  bool takeLogs=false)
{
    ofstream fout;
    fout.open(fileName);
    // check that file can be opened
    if (!fout.good())
    {
        cerr << "\n\nError creating " << fileName << endl;
        return;
    }
    
    unsigned int dataSize = nrows*ncols*nslabs;
    
    fout << "# vtk DataFile Version 1.0\n"
        << "title\n"
        << "ASCII\n"
        << "DATASET STRUCTURED_POINTS\n"
        << "DIMENSIONS " << nrows << " " << ncols << " "  << nslabs << '\n'
        << "ORIGIN 0.0 0.0 0.0\n"
        << "ASPECT_RATIO 1.0 1.0 1.0\n\n"
    	<< "POINT_DATA " << dataSize << '\n'
    	<< "SCALARS scalars float\n"
    	<< "LOOKUP_TABLE default\n";
    
    if (takeLogs)
    	for (unsigned int i=0; begin != end; begin++, i++) 
    	{
    	    fout << log(1+*begin) << " ";
    	    if (i==nslabs) 
    	    {
    	    	i=0;
    	    	fout << '\n';
    	    }
    	}
    else
    	for (unsigned int i=0; begin != end; begin++, i++) 
    	{
    	    fout << *begin << " ";
    	    if (i==nslabs) 
    	    {
    	    	i=0;
    	    	fout << '\n';
    	    }
    	}
}    
    
// Perform a shift so that the origin lies in the centre. 
//
template <class Iterator>
void Shift(Iterator begin, 
           Iterator tempBegin,
    	   unsigned int nrows,
           unsigned int ncols,
	   unsigned int nslabs)
{
    unsigned int dataSize = nrows*ncols*nslabs;

    // Copy data to temp container
    //
    Iterator i=begin, j = tempBegin;   
    for (unsigned int k=0; k<dataSize; k++)
    	*j++ = *i++;

    
    unsigned int row,col,slab,x,y,z,iarr,jarr;
    for (row=0; row<nrows; row++)
    {
    	if (row < nrows/2)
    	    x = nrows/2 + row;
    	else
    	    x = row - nrows/2;
    	    
    	for (col=0; col<ncols; col++)
    	{
    	    if (col < ncols/2)
    		y = ncols/2 + col;
    	    else
    		y = col - ncols/2;
    		
    	    for (slab=0; slab<nslabs; slab++)
    	    {
    	    	if (slab < nslabs/2)
    	    	    z = nslabs/2 + slab;
    	    	else
    	    	    z = slab - nslabs/2;

    	    	iarr = ncols*nslabs*row + nslabs*col + slab;
    	    	jarr = ncols*nslabs*x + nslabs*y + z;

    	    	begin[iarr] = tempBegin[jarr];
    	    }
    	}
    }   	   
}

void Series_3Dtest()
{

    // Create containers to contain test data
    //
    
    typedef vector< complex_type > Container;
    unsigned int dim = 32, nrows = dim, ncols = dim, nslabs = dim,
                 dataSize = nrows*ncols*nslabs;
    Container sin_freq_1(dataSize),
    	      sin_freq_2(dataSize),
    	      sin_freq_3(dataSize),
    	      sin_freq_sum(dataSize);
 
    // Put sin series with 4 different frequencies into containers
    //
    double rowstep, colstep, slabstep, stepSize = 2.0 * M_PI / dim, start=-M_PI;
    unsigned int row,col,slab,iarr;
    for (row=0, iarr=0, rowstep=start; row<nrows; row++, rowstep+=stepSize)
    	for (col=0, colstep=start; col<ncols; col++, colstep+=stepSize)
    	    for (slab=0, slabstep=start; slab<nslabs; slab++,
    	         slabstep+=stepSize, iarr++)
    	    {
    	      	sin_freq_1[iarr] =complex_type(sin(1*rowstep)*
                                              sin(1*colstep)*sin(1*slabstep));
    	      	sin_freq_2[iarr] =complex_type(sin(2*rowstep)*
                                              sin(2*colstep)*sin(2*slabstep));
    	      	sin_freq_3[iarr] =complex_type(sin(3*rowstep)*
                                              sin(3*colstep)*sin(3*slabstep));
    	    	sin_freq_sum[iarr] = sin_freq_1[iarr] + 
    	                             sin_freq_2[iarr] +
    	                             sin_freq_3[iarr];
    	    }
    	    
    // Output input data to a file
    //    
    ofstream fin;
    fin.open("fft3D_in.vtk");

    // check that file can be opened
    if (!fin)
    {
        cerr << "\n\nError creating fft3D_in.vtk" << endl;
        return;
    }
    fin << "# vtk DataFile Version 1.0\n"
        << "Sum of three sin functions of different frequencies "
        << "- test input to FFTTest.cpp\n"
        << "ASCII\n"
        << "DATASET STRUCTURED_POINTS\n"
        << "DIMENSIONS " << nrows << " " << ncols << " "  << nslabs << '\n'
        << "ORIGIN 0.0 0.0 0.0\n"
        << "ASPECT_RATIO 1.0 1.0 1.0\n\n"
    	<< "POINT_DATA " << dataSize << '\n'
    	<< "SCALARS scalars float\n"
    	<< "LOOKUP_TABLE default\n";
    for (row=0, iarr=0; row<nrows; row++)
    	for (col=0; col<ncols; col++)
    	{
    	    for (slab=0; slab<nslabs; slab++, iarr++)
    	    	fin << real(sin_freq_sum[iarr]) << " ";
    	    fin << '\n';
    	}
    
    
    
    // Create containers to contain Fourier transforms of the input data
    //
    Container output_1(dataSize),
    	      output_2(dataSize),
    	      output_3(dataSize),
    	      output_sum(dataSize);

    // Do the FFT calculations
    //    	      

    fourier_transform(sin_freq_1.begin(),output_1.begin(), nrows, ncols, nslabs);
    fourier_transform(sin_freq_2.begin(),output_2.begin(), nrows, ncols, nslabs);
    fourier_transform(sin_freq_3.begin(),output_3.begin(), nrows, ncols, nslabs);
    fourier_transform(sin_freq_sum.begin(),output_sum.begin(), nrows,ncols,nslabs);
   

    // Calculate the power spectrum of the resulting complex numbers to 
    // remove the imaginary part. Compare the theoretical result with the
    // calculated result
    //
    vector<double> powerSpectrum_sum(dataSize);
    double theoretical_result, diff, maxdiff = 0.0;

    for (unsigned int i=0; i<dataSize; i++)
    {
    	theoretical_result  = norm(output_1[i]+output_2[i]+output_3[i]);
    	powerSpectrum_sum[i] = norm(output_sum[i]);
    	diff = fabs(theoretical_result - powerSpectrum_sum[i]);
    	if (diff > maxdiff) maxdiff = diff;
    }    
    // Output the error to the screen
    //
    cout << "3D Sin series test: max difference from ideal result = " << maxdiff << endl;
    
    vector<double> temp(dataSize);
    Shift(powerSpectrum_sum.begin(), temp.begin(), nrows, ncols, nslabs);
    ArrayToFile3D("fft3D_out.vtk",
                  powerSpectrum_sum.begin(),
                  powerSpectrum_sum.end(),
                  nrows, ncols, nslabs, true);
    
}


//..............................................................................

void SinglePoint_3Dtest()
{

    // Create containers to contain test data
    //
    
    typedef vector< complex_type > Container;
    unsigned int dim = 32, nrows = dim, ncols = dim, nslabs = dim,
                 dataSize = nrows*ncols*nslabs;

    // Create container of all zero's
    //
    Container testData(dataSize, (complex_type) 0.0);
    Container testResults(dataSize);
    
    // Make first element equal to 1
    //
    testData[0] = (complex_type) 1;

    // Do the FFT calculations
    //    	      
    fourier_transform(testData.begin(),testResults.begin(), nrows, ncols, nslabs);
    
    // The resulting transform should be evenly spread out.
    //

    
    vector<double> powerSpec(dataSize);
    double normal = norm(testResults[0]), max=normal, min=normal;
    for (unsigned int i=0; i<dataSize; i++)
    {
    	normal  = norm(testResults[i]);
    	if (normal < min) min = normal;
    	if (normal > max) max = normal;
    	powerSpec[i] = normal;
    }

    vector<double> temp(dataSize);
    Shift(powerSpec.begin(), temp.begin(), nrows, ncols, nslabs);
    ArrayToFile3D("fft3D_point.vtk",
                  powerSpec.begin(),
                  powerSpec.end(),
                  nrows, ncols, nslabs, true);
    
    complex_type error = max-min;
    
    TestResults("3D Single point transform",error);
}

//..............................................................................

int main(int argc, char* argv[])
{

    SinglePoint_3Dtest();
    Series_1Dtest();
    Series_2Dtest();
    Series_3Dtest();
    
    unsigned int dim1 = 30, dim2 = 20, dim3 = 7;
    if (argc > 1) dim1 = atoi(argv[1]);
    if (argc > 2) dim2 = atoi(argv[2]);
    if (argc > 3) dim3 = atoi(argv[3]);
    
    cout << "\nSize of 1D containers: " << dim1 << "\n\n";
    
    Check1DVectorInput(dim1);
    Check1DArrayInput(dim1);

    cout << "\nSize of 2D containers: " << dim1 << " x " << dim2 << "\n\n";
    
    Check2DVectorInput(dim1,dim2);
    Check2DArrayInput(dim1,dim2);

    cout << "\nSize of 3D containers: " 
    	 << dim1 << " x " 
    	 << dim2 << " x " 
    	 << dim3 << "\n\n";
    
    Check3DVectorInput(dim1,dim2,dim3);
    Check3DArrayInput(dim1,dim2,dim3);

    return 0;
}