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.6 - 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
53using namespace std ; // Let us omit the std:: prefix for every STL class.
54using namespace std::complex_literals ; // Let's us do complex c = 1 +2i
55
56#include "FFT.hpp" // Our own FFT stuff.
57
58// Pick just one.
59#define SIMPLE_WAY_TO_READ_NUMBERS_FROM_FILE
60//#define COMPLICATED_WAY_TO_READ_NUMBERS_FROM_FILE
61//#define HARDCODED_DATA
62
63string legalNotice
64(
65 "\n"
66 "FFT Version 1.2 - An FFT utility in C++.\n"
67 "Copyright (C) 2005-2024 by Sean Erik O'Connor. All Rights Reserved.\n"
68 "\n"
69 "FFT comes with ABSOLUTELY NO WARRANTY; for details see the\n"
70 "GNU General Public License. This is free software, and you are welcome\n"
71 "to redistribute it under certain conditions; see the GNU General Public License\n"
72 "for details.\n\n"
73) ;
74
75
76/*=============================================================================
77|
78| NAME
79|
80| main
81|
82| DESCRIPTION
83|
84| This program tests the FFT implementation by performing forward and
85| inverse transforms on the data and testing error conditions.
86|
87| INPUT
88|
89| File containing
90|
91| - Number of complex points in the transform.
92| - The data points (complex numbers).
93|
94| OUTPUT
95|
96| - DFT of the input.
97| - inverse DFT of the first DFT.
98|
99| EXAMPLE
100|
101| From the command line, call
102|
103| $ FastFourierTransform fftIn.txt fftOut.txt
104|
105| where the input file fftIn.txt contains
106|
107| 4
108| (1, 1)
109| (2, 2)
110| (3, 3)
111| (4, 4)
112|
113| and the output file fftOut.txt contains
114|
115| num_points = 4
116| Point num. 0= (1,1)
117| Point num. 1= (2,2)
118| Point num. 2= (3,3)
119| Point num. 3= (4,4)
120| The transform of the input data is . . .
121|
122| Point num. 0 = (5,5)
123| Point num. 1 = (-2,-1.11022e-016)
124| Point num. 2 = (-1,-1)
125| Point num. 3 = (0,-2)
126| The inverse transform of the transform above is . . .
127|
128| Point num. 0 = (1,1)
129| Point num. 1 = (2,2)
130| Point num. 2 = (3,3)
131| Point num. 3 = (4,4)
132| Point num. 0 = (1,1)
133| Point num. 1 = (2,2)
134| Point num. 2 = (3,3)
135| Point num. 3 = (4,4)
136|
137| Calling FFT with too many points.
138| FFT has too many points for its sin/cos table
139|
140+============================================================================*/
141
142int main( int argc, char * argv[] )
143{
144 cout << legalNotice ;
145
146 int num_points = 0 ; // Number of points in the transform.
147 int i = 0 ;
148
149 // Open file streams for input and output getting file names
150 // from the command line.
151 if (argc != 3)
152 {
153 cerr << "Error: Missing the names of the input and output "
154 "files or too many files on the command line." ;
155 return 1 ;
156 }
157
158 char * inputFileName = argv[ 1 ] ;
159 char * outputFileName = argv[ 2 ] ;
160
161 fstream fin ( inputFileName, ios::in ) ;
162 fstream fout( outputFileName, ios::out ) ;
163
164 // Set a really high precision for our printouts.
165 // Of course some binary numbers will have repeating decimals, but we are checking to see if we get the same rounding results
166 // different computer architectures.
167 fout << scientific << setprecision( 100 ) << setw( 100 );
168
169 if (!(fin.is_open() && fout.is_open()))
170 {
171 cerr << "Error: Problem opening the input or output files. "
172 "Check the file names and permissions." ;
173 return 1 ;
174 }
175
176 try
177 {
178 //////////////////////////////////////////////////////////////////////////
179 /// Single precision FFT
180 //////////////////////////////////////////////////////////////////////////
181
182 // Create our DFT object.
183 FastFourierTransform<float> dft_single ;
184
185 // Clear the fft object.
186 dft_single.clear() ;
187
188 // Read lines until end of file.
189 while (!fin.eof())
190 {
191 // Read the number of points in the transform.
192 fin >> num_points ;
193
194 // Check to see that the number of points in the FFT is valid.
195 if (num_points <= 0)
196 throw FastFourierTransformException( "The number of points is out of range.\n" ) ;
197
198 // Read in the complex data points from file into the vector.
199 for (i = 0 ; i < num_points ; ++i)
200 {
201 complex<float> c ;
202 fin >> c ;
203
204 dft_single.push_back( c ) ;
205 }
206 }
207
208 fout << "* * * * * * Single precision FFT" << endl ;
209
210 // Print out the points just read in.
211 fout << "\n\n\nNumber of points = " << num_points << endl ;
212
213 for (i = 0 ; i <= num_points - 1 ; ++i)
214 fout << "Point num. " << i << " = " << dft_single[ i ] << endl ;
215
216 // Try a forward FFT.
217 dft_single.fft( true ) ;
218
219 // Print the transform results.
220 fout << "\n\nThe discrete Fourier transform of the input data is . . . \n\n" ;
221
222 for (i = 0 ; i <= num_points - 1 ; ++i)
223 fout << "Point num. " << i << " = " << dft_single[ i ] << endl ;
224
225 // Try an inverse transform to get the input data back, with roundoff error.
226 dft_single.fft( false ) ;
227
228 // Print the inverse transform.
229 fout << "\n\nThe inverse transform of the transform above is . . .\n\n" ;
230
231 for (i = 0 ; i <= num_points - 1 ; ++i)
232 fout << "Point num. " << i << " = " << dft_single[ i ] << endl ;
233
234 //////////////////////////////////////////////////////////////////////////
235 /// Double precision FFT
236 //////////////////////////////////////////////////////////////////////////
237
238 fout << endl << endl << endl << "* * * * * * Double precision FFT" << endl ;
239
240 // Create our DFT object.
241 FastFourierTransform<double> dft_double ;
242
243 // Clear the fft object.
244 dft_double.clear() ;
245
246#if defined( SIMPLE_WAY_TO_READ_NUMBERS_FROM_FILE ) || defined( COMPLICATED_WAY_TO_READ_NUMBERS_FROM_FILE )
247 // Rewind the input file to the begining.
248 fin.clear() ; // Forget we hit the end of file.
249 fin.seekg( 0, ios::beg ) ; // Move to the start of the file.
250#endif // defined( SIMPLE_WAY_TO_READ_NUMBERS_FROM_FILE ) || defined( COMPLICATED_WAY_TO_READ_NUMBERS_FROM_FILE )
251
252#ifdef COMPLICATED_WAY_TO_READ_NUMBERS_FROM_FILE
253 // Store the lines of file as a vector of strings.
254 vector<string> buffer ;
255 buffer.clear() ;
256
257 // Add each line of the file to the vector.
258 string line_of_file ;
259 while (getline( fin, line_of_file ))
260 buffer.push_back( line_of_file ) ;
261
262 // Now the whole file is in memory.
263 // Parse the first line which is the number of points.
264 istringstream is ;
265 is.clear() ;
266 is.str( buffer[ 0 ] ) ;
267 is >> num_points ;
268
269 // Check to see that the number of points in the FFT is valid.
270 if (num_points <= 0)
271 throw FastFourierTransformException( "The number of points is out of range.\n" ) ;
272
273 // Now parse all the other lines which are the data points.
274 int line_num = 1 ;
275 while (line_num < buffer.size())
276 {
277 // Get the next line in memory and advance.
278 string line_of_file = buffer[ line_num++ ] ;
279
280 // Parse the complex number.
281 istringstream is ;
282 is.clear() ;
283 is.str( line_of_file ) ;
284 complex<double> c ;
285 is >> c ;
286
287 dft_double.push_back( c ) ;
288 }
289#endif // COMPLICATED_WAY_TO_READ_NUMBERS_FROM_FILE
290
291#ifdef SIMPLE_WAY_TO_READ_NUMBERS_FROM_FILE
292 if (!fin.eof())
293 {
294 // Read the number of points in the transform.
295 fin >> num_points ;
296
297 // Check to see that the number of points in the FFT is valid.
298 if (num_points <= 0)
299 {
300 ostringstream os ;
301 os << "The number of points = " << num_points << " is out of range.\n" ;
302 throw FastFourierTransformException( os.str() ) ;
303 }
304 }
305
306 // Read in the complex data points from file into the vector.
307 for (int index = 1 ; index <= num_points ; ++index)
308 {
309 complex<double> c = 0.0 + 0.0i ; // Floating literals default to double.
310
311 if (!fin.eof())
312 fin >> c ;
313 else
314 {
315 ostringstream os ;
316 os << "The number of data points was too short.\n" ;
317 throw FastFourierTransformException( os.str() ) ;
318 }
319
320 dft_double.push_back( c ) ;
321 }
322#endif // SIMPLE_WAY_TO_READ_NUMBERS_FROM_FILE
323
324#ifdef HARDCODED_DATA
325 num_points = 16 ;
326 dft_double.push_back( 1.0 + 1.0i ) ;
327 dft_double.push_back( 2.0 + 2.0i ) ;
328 dft_double.push_back( 3.0 + 3.0i ) ;
329 dft_double.push_back( 4.0 + 4.0i ) ;
330 dft_double.push_back( 5.0 + 5.0i ) ;
331 dft_double.push_back( 6.0 + 6.0i ) ;
332 dft_double.push_back( 7.0 + 7.0i ) ;
333 dft_double.push_back( 8.0 + 8.0i ) ;
334 dft_double.push_back( 9.0 + 9.0i ) ;
335 dft_double.push_back( 10.0 + 10.0i ) ;
336 dft_double.push_back( 11.0 + 11.0i ) ;
337 dft_double.push_back( 12.0 + 12.0i ) ;
338 dft_double.push_back( 13.0 + 13.0i ) ;
339 dft_double.push_back( 14.0 + 14.0i ) ;
340 dft_double.push_back( 15.0 + 15.0i ) ;
341 dft_double.push_back( 16.0 + 16.0i ) ;
342#endif // HARDCODED_DATA
343
344 // Print out the points just read in.
345 fout << "\n\n\nNumber of points = " << num_points << endl ;
346
347 for (i = 0 ; i <= num_points - 1 ; ++i)
348 fout << "Point num. " << i << "= " << dft_double[ i ] << endl ;
349
350 // Try a forward FFT.
351 dft_double.fft( true ) ;
352
353 // Print the transform results.
354 fout << "\n\nThe discrete Fourier transform of the input data is . . . \n\n" ;
355
356 for (i = 0 ; i <= num_points - 1 ; ++i)
357 fout << "Point num. " << i << " = " << dft_double[ i ] << endl ;
358
359 // Try an inverse transform to get the input data back, with roundoff error.
360 dft_double.fft( false ) ;
361
362 // Print the inverse transform.
363 fout << "\n\nThe inverse transform of the transform above is . . .\n\n" ;
364
365 for (i = 0 ; i <= num_points - 1 ; ++i)
366 fout << "Point num. " << i << " = " << dft_double[ i ] << endl ;
367
368 // =======================================
369 // Test the copy constructor.
370 // =======================================
371
372 fout << "\n\nTest the copy constructor" << endl ;
373
374 FastFourierTransform copy_of_dft_single( dft_single ) ;
375 for (i = 0 ; i <= num_points - 1 ; ++i)
376 fout << "Point num. " << i << " = " << copy_of_dft_single[ i ] << endl ;
377
378 // =======================================
379 // Floating point implementation details.
380 // =======================================
381
382 // Set a really high precision for our printouts.
383 fout << scientific << setprecision( 50 ) << setw( 60 ) << endl ;
384 fout << "= floating point epsilons =" << endl ;
385 fout << "single precision floating point epsilon = " << numeric_limits<float>::epsilon() << endl ;
386 fout << "single precision test: (1.0f + 0.999999f * epsilon) - 1.0f = " << ((1.0f + 0.5f * numeric_limits<float>::epsilon()) - 1.0f) << endl ;
387 fout << "single precision test: (1.0f + epsilon - 1.0f = " << ((1.0f + numeric_limits<float>::epsilon()) - 1.0f) << endl << endl ;
388
389 fout << "double precision floating point epsilon = " << numeric_limits<double>::epsilon() << endl ;
390 fout << "double precision test: (1.0f + 0.999999 * epsilon) - 1.0f = " << ((1.0f + 0.5f * numeric_limits<double>::epsilon()) - 1.0f) << endl ;
391 fout << "double precision test: (1.0f + epsilon - 1.0f = " << ((1.0f + numeric_limits<double>::epsilon()) - 1.0f) << endl;
392
393 // =======================================
394 // Test error conditions. Only the first will execute due to
395 // the throw, so comment out and recompile to test them all.
396 // =======================================
397
398 // Number of points is less than 2.
399 fout << "\n\nCalling FFT with 1 point.\n" ;
400
401 // Make sure we flush all the data to file, when we're debugging. Then close the files.
402 fout << flush ;
403 fin.close() ;
404 fout.close() ;
405
406 FastFourierTransform<double> dft_single2 ;
407 dft_single2.push_back( complex<double>(1.0, 2.0) ) ;
408 dft_single2.fft( true ) ;
409 }
410
411 // Catch exceptions and print their text explanations.
412 // Clean up resources.
413 catch( FastFourierTransformException & e )
414 {
415 fout << "FFT exception: " << e.what() << endl ;
416
417 fin.close() ;
418 fout.close() ;
419
420 return 1 ;
421 }
422 catch( runtime_error & e )
423 {
424 fout << "Unexpected runtime exception: Please email the author." << e.what() << endl ;
425
426 fin.close() ;
427 fout.close() ;
428
429 return 1 ;
430 }
431 catch( ... )
432 {
433 fout << "Unexpected exception: Please email the author." << endl ;
434
435 fin.close() ;
436 fout.close() ;
437
438 return 1 ;
439 }
440
441
442 return 0 ;
443
444} /* ====================== end of function main ============================ */
