1 /*============================================================================= 2 | 3 | NAME 4 | 5 | testFFT.cpp 6 | 7 | DESCRIPTION 8 | 9 | Test program for FFT.cpp 10 | 11 | AUTHOR 12 | 13 | Sean O'Connor 14 | 15 | LEGAL 16 | 17 | FFT Version 1.4 - An FFT utility library in C++. 18 | Copyright (C) 2005-2024 by Sean Erik O'Connor. All Rights Reserved. 19 | 20 | This program is free software: you can redistribute it and/or modify 21 | it under the terms of the GNU General Public License as published by 22 | the Free Software Foundation, either version 3 of the License, or 23 | (at your option) any later version. 24 | 25 | This program is distributed in the hope that it will be useful, 26 | but WITHOUT ANY WARRANTY; without even the implied warranty of 27 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 28 | GNU General Public License for more details. 29 | 30 | You should have received a copy of the GNU General Public License 31 | along with this program. If not, see <http://www.gnu.org/licenses/>. 32 | 33 | The author's address is seanerikoconnor!AT!gmail!DOT!com 34 | with !DOT! replaced by . and the !AT! replaced by @ 35 | 36 +============================================================================*/ 37 38 /*------------------------------------------------------------------------------ 39 | Include Files | 40 ------------------------------------------------------------------------------*/ 41 42 #include <iostream>// Basic stream I/O. 43 #include <cmath>// Basic math functions. 44 #include <complex>// Complex data type and operations. 45 #include <vector>// STL vector class. 46 #include <fstream>// File stream I/O. 47 #include <algorithm>// Iterators. 48 #include <string>// STL string class. 49 #include <stdexcept>// Exceptions. 50 #include <limits>// Numeric limits, e.g. floating point. 51 #include <iomanip>// Floating point formating. 52 53 usingnamespace std ; // Let us omit the std:: prefix for every STL class. 54 55 #include "FFT.h"// Our own FFT stuff. 56 57 string legalNotice
58 (
59 "\n" 60 "FFT Version 1.2 - An FFT utility in C++.\n" 61 "Copyright (C) 2005-2024 by Sean Erik O'Connor. All Rights Reserved.\n" 62 "\n" 63 "FFT comes with ABSOLUTELY NO WARRANTY; for details see the\n" 64 "GNU General Public License. This is free software, and you are welcome\n" 65 "to redistribute it under certain conditions; see the GNU General Public License\n" 66 "for details.\n\n" 67 ) ;
68 69 70 /*============================================================================= 71 | 72 | NAME 73 | 74 | main 75 | 76 | DESCRIPTION 77 | 78 | This program tests the FFT implementation by performing forward and 79 | inverse transforms on the data and testing error conditions. 80 | 81 | INPUT 82 | 83 | File containing 84 | 85 | - Number of complex points in the transform. 86 | - The data points (complex numbers). 87 | 88 | OUTPUT 89 | 90 | - DFT of the input. 91 | - inverse DFT of the first DFT. 92 | 93 | EXAMPLE 94 | 95 | From the command line, call 96 | 97 | $ FastFourierTransform fftIn.txt fftOut.txt 98 | 99 | where the input file fftIn.txt contains100 |101 | 4102 | (1, 1)103 | (2, 2)104 | (3, 3)105 | (4, 4)106 | 107 | and the output file fftOut.txt contains108 |109 | num_points = 4110 | Point num. 0= (1,1)111 | Point num. 1= (2,2)112 | Point num. 2= (3,3)113 | Point num. 3= (4,4)114 | The transform of the input data is . . . 115 | 116 | Point num. 0 = (5,5)117 | Point num. 1 = (-2,-1.11022e-016)118 | Point num. 2 = (-1,-1)119 | Point num. 3 = (0,-2)120 | The inverse transform of the transform above is . . .121 | 122 | Point num. 0 = (1,1)123 | Point num. 1 = (2,2)124 | Point num. 2 = (3,3)125 | Point num. 3 = (4,4)126 | Point num. 0 = (1,1)127 | Point num. 1 = (2,2)128 | Point num. 2 = (3,3)129 | Point num. 3 = (4,4)130 |131 | Calling FFT with too many points.132 | FFT has too many points for its sin/cos table133 |134 +============================================================================*/135 136 int main( int argc, char * argv[] )
137 {
138 cout << legalNotice ;
139 140 int num_points = 0 ; // Number of points in the transform.141 int i = 0 ;
142 143 // Open file streams for input and output getting file names144 // from the command line.145 if (argc != 3)
146 {
147 cerr << "Error: Missing the names of the input and output "148 "files or too many files on the command line." ;
149 return1 ;
150 }
151 152 char * inputFileName = argv[ 1 ] ;
153 char * outputFileName = argv[ 2 ] ;
154 155 fstream fin ( inputFileName, ios::in ) ;
156 fstream fout( outputFileName, ios::out ) ;
157 158 // Set a really high precision for our printouts.159 // Of course some binary numbers will have repeating decimals, but we are checking to see if we get the same rounding results160 // different computer architectures.161 fout << scientific << setprecision( 70 ) << setw( 80 );
162 163 if (!(fin.is_open() && fout.is_open()))
164 {
165 cerr << "Error: Problem opening the input or output files. "166 "Check the file names and permissions." ;
167 return1 ;
168 }
169 170 try171 {
172 //////////////////////////////////////////////////////////////////////////173 /// Single precision FFT174 //////////////////////////////////////////////////////////////////////////175 176 // Create our DFT object.177 FastFourierTransform<float> dft_single ;
178 179 // Clear the fft object.180 dft_single.clear() ;
181 182 // Read lines until end of file.183 while (!fin.eof())
184 {
185 // Read the number of points in the transform.186 fin >> num_points ;
187 188 // Check to see that the number of points in the FFT is valid.189 if (num_points <= 0)
190 throw FastFourierTransformException( "The number of points is out of range.\n" ) ;
191 192 // Read in the complex data points from file into the vector.193 for (i = 0 ; i < num_points ; ++i)
194 {
195 complex<float> c ;
196 fin >> c ;
197 198 dft_single.push_back( c ) ;
199 }
200 }
201 202 fout << "* * * * * * Single precision FFT" << endl ;
203 204 // Print out the points just read in.205 fout << "\n\n\nNumber of points = " << num_points << endl ;
206 207 for (i = 0 ; i <= num_points - 1 ; ++i)
208 fout << "Point num. " << i << "= " << dft_single[ i ] << endl ;
209 210 // Try a forward FFT.211 dft_single.fft( true ) ;
212 213 // Print the transform results.214 fout << "\n\nThe discrete Fourier transform of the input data is . . . \n\n" ;
215 216 for (i = 0 ; i <= num_points - 1 ; ++i)
217 fout << "Point num. " << i << " = " << dft_single[ i ] << endl ;
218 219 // Try an inverse transform to get the input data back, with roundoff error.220 dft_single.fft( false ) ;
221 222 // Print the inverse transform.223 fout << "\n\nThe inverse transform of the transform above is . . .\n\n" ;
224 225 for (i = 0 ; i <= num_points - 1 ; ++i)
226 fout << "Point num. " << i << " = " << dft_single[ i ] << endl ;
227 228 //////////////////////////////////////////////////////////////////////////229 /// Double precision FFT230 //////////////////////////////////////////////////////////////////////////231 232 // Rewind the input file to the begining.233 fin.clear() ; // Forget we hit the end of file.234 fin.seekg( 0, ios::beg ) ; // Move to the start of the file.235 236 // Create our DFT object.237 FastFourierTransform<double> dft_double ;
238 239 // Clear the fft object.240 dft_double.clear() ;
241 242 // Read lines until end of file.243 while (!fin.eof())
244 {
245 // Read the number of points in the transform.246 fin >> num_points ;
247 248 // Check to see that the number of points in the FFT is valid.249 if (num_points <= 0)
250 throw FastFourierTransformException( "The number of points is out of range.\n" ) ;
251 252 // Read in the complex data points from file into the vector.253 for (i = 0 ; i < num_points ; ++i)
254 {
255 complex<double> c ;
256 fin >> c ;
257 258 dft_double.push_back( c ) ;
259 }
260 }
261 262 fout << "* * * * * * Double precision FFT" << endl ;
263 264 // Print out the points just read in.265 fout << "\n\n\nNumber of points = " << num_points << endl ;
266 267 for (i = 0 ; i <= num_points - 1 ; ++i)
268 fout << "Point num. " << i << "= " << dft_double[ i ] << endl ;
269 270 // Try a forward FFT.271 dft_double.fft( true ) ;
272 273 // Print the transform results.274 fout << "\n\nThe discrete Fourier transform of the input data is . . . \n\n" ;
275 276 for (i = 0 ; i <= num_points - 1 ; ++i)
277 fout << "Point num. " << i << " = " << dft_double[ i ] << endl ;
278 279 // Try an inverse transform to get the input data back, with roundoff error.280 dft_double.fft( false ) ;
281 282 // Print the inverse transform.283 fout << "\n\nThe inverse transform of the transform above is . . .\n\n" ;
284 285 for (i = 0 ; i <= num_points - 1 ; ++i)
286 fout << "Point num. " << i << " = " << dft_double[ i ] << endl ;
287 288 // =======================================289 // Test the copy constructor.290 // =======================================291 292 fout << "\n\nTest the copy constructor" << endl ;
293 294 FastFourierTransform copy_of_dft_single( dft_single ) ;
295 for (i = 0 ; i <= num_points - 1 ; ++i)
296 fout << "Point num. " << i << " = " << copy_of_dft_single[ i ] << endl ;
297 298 // =======================================299 // Floating point implementation details.300 // =======================================301 302 // Set a really high precision for our printouts.303 fout << scientific << setprecision( 50 ) << setw( 60 ) << endl ;
304 fout << "= floating point epsilons =" << endl ;
305 fout << "single precision floating point epsilon = " << numeric_limits<float>::epsilon() << endl ;
306 fout << "single precision test: (1.0f + 0.999999f * epsilon) - 1.0f = " << ((1.0f + 0.5f * numeric_limits<float>::epsilon()) - 1.0f) << endl ;
307 fout << "single precision test: (1.0f + epsilon - 1.0f = " << ((1.0f + numeric_limits<float>::epsilon()) - 1.0f) << endl << endl ;
308 309 fout << "double precision floating point epsilon = " << numeric_limits<double>::epsilon() << endl ;
310 fout << "double precision test: (1.0f + 0.999999 * epsilon) - 1.0f = " << ((1.0f + 0.5f * numeric_limits<double>::epsilon()) - 1.0f) << endl ;
311 fout << "double precision test: (1.0f + epsilon - 1.0f = " << ((1.0f + numeric_limits<double>::epsilon()) - 1.0f) << endl;
312 313 // =======================================314 // Test error conditions. Only the first will execute due to315 // the throw, so comment out and recompile to test them all.316 // =======================================317 318 // Number of points is less than 2.319 fout << "\n\nCalling FFT with 1 point.\n" ;
320 321 // Make sure we flush all the data to file, when we're debugging. Then close the files.322 fout << flush ;
323 fin.close() ;
324 fout.close() ;
325 326 FastFourierTransform<double> dft_single2 ;
327 dft_single2.push_back( complex<double>(1.0, 2.0) ) ;
328 dft_single2.fft( true ) ;
329 }
330 331 // Catch exceptions and print their text explanations.332 // Clean up resources.333 catch( FastFourierTransformException & e )
334 {
335 fout << "FFT exception: " << e.what() << endl ;
336 337 fin.close() ;
338 fout.close() ;
339 340 return1 ;
341 }
342 catch( runtime_error & e )
343 {
344 fout << "Unexpected runtime exception: Please email the author." << e.what() << endl ;
345 346 fin.close() ;
347 fout.close() ;
348 349 return1 ;
350 }
351 catch( ... )
352 {
353 fout << "Unexpected exception: Please email the author." << endl ;
354 355 fin.close() ;
356 fout.close() ;
357 358 return1 ;
359 }
360 361 362 return0 ;
363 364 } /* ====================== end of function main ============================ */
