/*=============================================================================
|
| 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 <http://www.gnu.org/licenses/>.
|     
|     The author's address is seanerikoconnor!AT!gmail!DOT!com
|     with !DOT! replaced by . and the !AT! replaced by @
|     
+============================================================================*/

/*------------------------------------------------------------------------------
|                                Include Files                                 |
------------------------------------------------------------------------------*/

#include <iostream>     // Basic stream I/O.
#include <cmath>        // Basic math functions.
#include <complex>      // Complex data type and operations.
#include <vector>       // STL vector class.
#include <fstream>      // File stream I/O.
#include <algorithm>    // Iterators.
#include <string>       // STL string class.
#include <stdexcept>    // Exceptions.
#include <limits>       // Numeric limits, e.g. floating point.
#include <iomanip>      // 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<float> 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<float> 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<double> 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<double> 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<float>::epsilon() << endl ;
        fout << "single precision test: (1.0f + 0.999999f * epsilon) - 1.0f = " << ((1.0f + 0.5f * numeric_limits<float>::epsilon()) - 1.0f) << endl ;
        fout << "single precision test: (1.0f + epsilon - 1.0f              = " << ((1.0f + numeric_limits<float>::epsilon()) - 1.0f) << endl << endl ;

        fout << "double precision floating point epsilon = " << numeric_limits<double>::epsilon() << endl ;
        fout << "double precision test: (1.0f + 0.999999 * epsilon) - 1.0f = " << ((1.0f + 0.5f * numeric_limits<double>::epsilon()) - 1.0f) << endl ;
        fout << "double precision test: (1.0f + epsilon - 1.0f             = " << ((1.0f + numeric_limits<double>::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<double> dft_single2 ;
        dft_single2.push_back( complex<double>(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 ============================ */