/*=============================================================================
|
| NAME
|
| testFFT.cpp
|
| DESCRIPTION
|
| Test program for FFT.cpp
|
| AUTHOR
|
| Sean O'Connor
|
| LEGAL
|
| FFT Version 1.5 - An FFT utility library in C++.
| Copyright (C) 2005-2024 by Sean Erik O'Connor. All Rights Reserved.
|
| This program is free software: you can redistribute it and/or modify
| it under the terms of the GNU General Public License as published by
| the Free Software Foundation, either version 3 of the License, or
| (at your option) any later version.
|
| This program is distributed in the hope that it will be useful,
| but WITHOUT ANY WARRANTY; without even the implied warranty of
| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
| GNU General Public License for more details.
|
| You should have received a copy of the GNU General Public License
| along with this program. If not, see .
|
| The author's address is seanerikoconnor!AT!gmail!DOT!com
| with !DOT! replaced by . and the !AT! replaced by @
|
+============================================================================*/
/*------------------------------------------------------------------------------
| Include Files |
------------------------------------------------------------------------------*/
#include // Basic stream I/O.
#include // Basic math functions.
#include // Complex data type and operations.
#include // STL vector class.
#include // File stream I/O.
#include // Iterators.
#include // STL string class.
#include // Exceptions.
#include // Numeric limits, e.g. floating point.
#include // Floating point formating.
using namespace std ; // Let us omit the std:: prefix for every STL class.
#include "FFT.h" // Our own FFT stuff.
string legalNotice
(
"\n"
"FFT Version 1.2 - An FFT utility in C++.\n"
"Copyright (C) 2005-2024 by Sean Erik O'Connor. All Rights Reserved.\n"
"\n"
"FFT comes with ABSOLUTELY NO WARRANTY; for details see the\n"
"GNU General Public License. This is free software, and you are welcome\n"
"to redistribute it under certain conditions; see the GNU General Public License\n"
"for details.\n\n"
) ;
/*=============================================================================
|
| NAME
|
| main
|
| DESCRIPTION
|
| This program tests the FFT implementation by performing forward and
| inverse transforms on the data and testing error conditions.
|
| INPUT
|
| File containing
|
| - Number of complex points in the transform.
| - The data points (complex numbers).
|
| OUTPUT
|
| - DFT of the input.
| - inverse DFT of the first DFT.
|
| EXAMPLE
|
| From the command line, call
|
| $ FastFourierTransform fftIn.txt fftOut.txt
|
| where the input file fftIn.txt contains
|
| 4
| (1, 1)
| (2, 2)
| (3, 3)
| (4, 4)
|
| and the output file fftOut.txt contains
|
| num_points = 4
| Point num. 0= (1,1)
| Point num. 1= (2,2)
| Point num. 2= (3,3)
| Point num. 3= (4,4)
| The transform of the input data is . . .
|
| Point num. 0 = (5,5)
| Point num. 1 = (-2,-1.11022e-016)
| Point num. 2 = (-1,-1)
| Point num. 3 = (0,-2)
| The inverse transform of the transform above is . . .
|
| Point num. 0 = (1,1)
| Point num. 1 = (2,2)
| Point num. 2 = (3,3)
| Point num. 3 = (4,4)
| Point num. 0 = (1,1)
| Point num. 1 = (2,2)
| Point num. 2 = (3,3)
| Point num. 3 = (4,4)
|
| Calling FFT with too many points.
| FFT has too many points for its sin/cos table
|
+============================================================================*/
int main( int argc, char * argv[] )
{
cout << legalNotice ;
int num_points = 0 ; // Number of points in the transform.
int i = 0 ;
// Open file streams for input and output getting file names
// from the command line.
if (argc != 3)
{
cerr << "Error: Missing the names of the input and output "
"files or too many files on the command line." ;
return 1 ;
}
char * inputFileName = argv[ 1 ] ;
char * outputFileName = argv[ 2 ] ;
fstream fin ( inputFileName, ios::in ) ;
fstream fout( outputFileName, ios::out ) ;
// Set a really high precision for our printouts.
// Of course some binary numbers will have repeating decimals, but we are checking to see if we get the same rounding results
// different computer architectures.
fout << scientific << setprecision( 70 ) << setw( 80 );
if (!(fin.is_open() && fout.is_open()))
{
cerr << "Error: Problem opening the input or output files. "
"Check the file names and permissions." ;
return 1 ;
}
try
{
//////////////////////////////////////////////////////////////////////////
/// Single precision FFT
//////////////////////////////////////////////////////////////////////////
// Create our DFT object.
FastFourierTransform dft_single ;
// Clear the fft object.
dft_single.clear() ;
// Read lines until end of file.
while (!fin.eof())
{
// Read the number of points in the transform.
fin >> num_points ;
// Check to see that the number of points in the FFT is valid.
if (num_points <= 0)
throw FastFourierTransformException( "The number of points is out of range.\n" ) ;
// Read in the complex data points from file into the vector.
for (i = 0 ; i < num_points ; ++i)
{
complex c ;
fin >> c ;
dft_single.push_back( c ) ;
}
}
fout << "* * * * * * Single precision FFT" << endl ;
// Print out the points just read in.
fout << "\n\n\nNumber of points = " << num_points << endl ;
for (i = 0 ; i <= num_points - 1 ; ++i)
fout << "Point num. " << i << "= " << dft_single[ i ] << endl ;
// Try a forward FFT.
dft_single.fft( true ) ;
// Print the transform results.
fout << "\n\nThe discrete Fourier transform of the input data is . . . \n\n" ;
for (i = 0 ; i <= num_points - 1 ; ++i)
fout << "Point num. " << i << " = " << dft_single[ i ] << endl ;
// Try an inverse transform to get the input data back, with roundoff error.
dft_single.fft( false ) ;
// Print the inverse transform.
fout << "\n\nThe inverse transform of the transform above is . . .\n\n" ;
for (i = 0 ; i <= num_points - 1 ; ++i)
fout << "Point num. " << i << " = " << dft_single[ i ] << endl ;
//////////////////////////////////////////////////////////////////////////
/// Double precision FFT
//////////////////////////////////////////////////////////////////////////
// Rewind the input file to the begining.
fin.clear() ; // Forget we hit the end of file.
fin.seekg( 0, ios::beg ) ; // Move to the start of the file.
// Create our DFT object.
FastFourierTransform dft_double ;
// Clear the fft object.
dft_double.clear() ;
// Read lines until end of file.
while (!fin.eof())
{
// Read the number of points in the transform.
fin >> num_points ;
// Check to see that the number of points in the FFT is valid.
if (num_points <= 0)
throw FastFourierTransformException( "The number of points is out of range.\n" ) ;
// Read in the complex data points from file into the vector.
for (i = 0 ; i < num_points ; ++i)
{
complex c ;
fin >> c ;
dft_double.push_back( c ) ;
}
}
fout << "* * * * * * Double precision FFT" << endl ;
// Print out the points just read in.
fout << "\n\n\nNumber of points = " << num_points << endl ;
for (i = 0 ; i <= num_points - 1 ; ++i)
fout << "Point num. " << i << "= " << dft_double[ i ] << endl ;
// Try a forward FFT.
dft_double.fft( true ) ;
// Print the transform results.
fout << "\n\nThe discrete Fourier transform of the input data is . . . \n\n" ;
for (i = 0 ; i <= num_points - 1 ; ++i)
fout << "Point num. " << i << " = " << dft_double[ i ] << endl ;
// Try an inverse transform to get the input data back, with roundoff error.
dft_double.fft( false ) ;
// Print the inverse transform.
fout << "\n\nThe inverse transform of the transform above is . . .\n\n" ;
for (i = 0 ; i <= num_points - 1 ; ++i)
fout << "Point num. " << i << " = " << dft_double[ i ] << endl ;
// =======================================
// Test the copy constructor.
// =======================================
fout << "\n\nTest the copy constructor" << endl ;
FastFourierTransform copy_of_dft_single( dft_single ) ;
for (i = 0 ; i <= num_points - 1 ; ++i)
fout << "Point num. " << i << " = " << copy_of_dft_single[ i ] << endl ;
// =======================================
// Floating point implementation details.
// =======================================
// Set a really high precision for our printouts.
fout << scientific << setprecision( 50 ) << setw( 60 ) << endl ;
fout << "= floating point epsilons =" << endl ;
fout << "single precision floating point epsilon = " << numeric_limits::epsilon() << endl ;
fout << "single precision test: (1.0f + 0.999999f * epsilon) - 1.0f = " << ((1.0f + 0.5f * numeric_limits::epsilon()) - 1.0f) << endl ;
fout << "single precision test: (1.0f + epsilon - 1.0f = " << ((1.0f + numeric_limits::epsilon()) - 1.0f) << endl << endl ;
fout << "double precision floating point epsilon = " << numeric_limits::epsilon() << endl ;
fout << "double precision test: (1.0f + 0.999999 * epsilon) - 1.0f = " << ((1.0f + 0.5f * numeric_limits::epsilon()) - 1.0f) << endl ;
fout << "double precision test: (1.0f + epsilon - 1.0f = " << ((1.0f + numeric_limits::epsilon()) - 1.0f) << endl;
// =======================================
// Test error conditions. Only the first will execute due to
// the throw, so comment out and recompile to test them all.
// =======================================
// Number of points is less than 2.
fout << "\n\nCalling FFT with 1 point.\n" ;
// Make sure we flush all the data to file, when we're debugging. Then close the files.
fout << flush ;
fin.close() ;
fout.close() ;
FastFourierTransform dft_single2 ;
dft_single2.push_back( complex(1.0, 2.0) ) ;
dft_single2.fft( true ) ;
}
// Catch exceptions and print their text explanations.
// Clean up resources.
catch( FastFourierTransformException & e )
{
fout << "FFT exception: " << e.what() << endl ;
fin.close() ;
fout.close() ;
return 1 ;
}
catch( runtime_error & e )
{
fout << "Unexpected runtime exception: Please email the author." << e.what() << endl ;
fin.close() ;
fout.close() ;
return 1 ;
}
catch( ... )
{
fout << "Unexpected exception: Please email the author." << endl ;
fin.close() ;
fout.close() ;
return 1 ;
}
return 0 ;
} /* ====================== end of function main ============================ */