Recently I implemented FOURIER() formula for LibreOffice Calc that computes Discrete Fourier Transform [DFT] of a real/complex data sequence. Computation is done using a couple of Fast Fourier Transform algorithms (all implemented from scratch). I’d like to thank Collabora Productivity for a fully funded hack week and lots of encouragement that enabled me to work on this feature.
The syntax of the formula is :
FOURIER(Array, GroupedByColumns, Inverse, Polar, MinimumMagnitude)
First and second arguments : First argument is the data array range. There are some restrictions on its shape. If the array is grouped by columns(rows), you need to indicate that by setting the second argument GroupedByColumns = TRUE(FALSE). In case of grouped by columns(rows) the “Array” can contain 1 or 2 columns(rows), where the first column(row) should contain the real part of the input sequence and second column(row) if present should contain the imaginary part of the input sequence. If there is only 1 column(row), the input sequence is assumed to be purely real. There is no restriction on the length of the input sequence. For example the length need not be an exact power of 2 like MS Excel’s Fourier analysis tool requires. The DFT computed by FOURIER formula is exact for the given length of the input sequence meaning it does not pad zeroes to the end of the input sequence to make the length an even power of 2.
The third argument “Inverse” is a boolean flag to indicate whether an inverse DFT needs to be computed. This argument is optional and the default value is FALSE.
The fourth argument “Polar” is a boolean flag to indicate whether the final output needs to be in polar coordinates. This argument is optional and the default value is FALSE.
The fifth argument “MinimumMagnitude” is a numeric argument and is used only if Polar=TRUE. All output components(frequency components if Inverse=FALSE) with magnitude less than MinimumMagnitude will be suppressed by a zero magnitude-phase entry. This is very useful when looking at the magnitude-phase spectrum of a signal (especially when its spectrum is sparse) because there is always some very tiny amount of rounding error when doing FFT algorithms and results in incorrect non-zero phase for some non-existent frequencies. By providing a suitable value to this parameter, these non-existent frequency components can be suppressed. The default value of this parameter is 0.0, so no suppression is done by default.
The result of this formula consists of two columns – first column contains the real parts (or the magnitudes if Polar=TRUE) and second column contains the imaginary parts (or the phases if Polar=TRUE).
Below is a screenshot of a sample usage of FOURIER() formula.
The data x[n] column was generated by the formula “=10*COS(2*PI()*COS(2*PI()*A2:A257/128))” (a typical frequency modulation example). Here a Phase-Magnitude spectrum is computed using FOURIER formula
=FOURIER(B2:B257,1,0,1, 1E-10)
Note that the Polar argument is set to 1, and MinimumMagnitude is set to 1E-10, to suppress the non-existent frequency components due to rounding errors. The plots here just visualize the data. The top plot shows the time-domain waveform and the next two plots represent the magnitude and phase spectra.
Fourier Analysis Tool
A Fourier analysis tool is also added to Statistics menu. All features of FOURIER formula are available graphically in this tool. Here is a screenshot of Fourier Analysis tool in action.
Implementation details of FOURIER formula
Basically two different Fast Fourier Transform (FFT) algorithms are implemented.
- Radix-2 Decimation in Time FFT : This is used when the length of input sequence is an even power of 2.
- Bluestein’s algorithm : This algorithm is used when the length of data sequence is not an even power of 2.
Both algorithms have asymptotic time complexity of O(N lgN), but Bluestein’s algorithm is much slower in practice. However many optimizations are in place to make the computation faster, especially when the input sequence is purely real.
Below are the patches that shows the evolution of the implementation of FOURIER formula over time.
- tdf#74664 : Adds FOURIER() formula
- tdf#74664 : Compute the phase correctly using atan2
- tdf#74664 : optimize the computation of twiddle factors.
- tdf#74664 : Unit test for FOURIER formula
- FOURIER : use Bluestein’s algorithm for N != 2^k
- More test cases for FOURIER formula
- tdf#74664: FOURIER: add 5th optional parameter MinimumMagnitude
Performance
Below are some perf measurements for various input cases. These numbers are only for the DFT computation part and does not include the time for writing data to the spreadsheet model. We don’t use multi-threading yet for FOURIER, so these are numbers for just one cpu-core.
Case#1 : Data length is a power of 2
| Real Input length | Time (ms) |
| 32768 | 2 |
| 65536 | 4 |
| 262144 | 21 |
| 1048576 | 185 |
| Complex Input length | Time (ms) |
| 32768 | 2 |
| 65536 | 5 |
| 262144 | 30 |
| 1048576 | 400 |
Case#2 : Data length is even but not a power of 2
| Real Input length | Time (ms) |
| 32766 | 8 |
| 65534 | 20 |
| 262142 | 105 |
| 1048574 | 1440 |
| Complex Input length | Time (ms) |
| 32766 | 15 |
| 65534 | 36 |
| 262142 | 313 |
| 1048574 | 3332 |
Case#3 : Data length is odd but not a power of 2
| Real Input length | Time (ms) |
| 32767 | 16 |
| 65535 | 39 |
| 262143 | 313 |
| 1048575 | 2729 |
| Complex Input length | Time (ms) |
| 32767 | 16 |
| 65535 | 38 |
| 262143 | 309 |
| 1048575 | 2703 |
Stacktrace 