Efficient Use of FourierTransform and InverseFourierTransform from the DiscreteTransforms Package
As mentioned in the DiscreteTransforms[FourierTransform] help page, the efficiency of computing the input transform using the fast methods (mintime and minstorage) is greatly enhanced if the data length for the transform is highly composite (that is, splits into small prime factors).
A (very) rough approximation of the expense of the FFT algorithms when used with N data points, and having a prime factorization of N=N1⁢N2...Nm with all Ni<Ni+1 is that the expense is proportional to N⋅Nm⋅logNlogNm. The DFT (discrete Fourier transform) algorithm simply has an expense of N2. (For additional information on the algorithm=DFT option, see DiscreteTransforms[FourierTransform].)
The expense of FFT algorithms is highly dependent on the factorization of the data length while the expense of the DFT algorithm is not. If N is prime, the expense of both algorithms is approximately the same.
The best case performance for the FFT algorithms can be observed when the data length is a power of 2. In this case the FFT expense is proportional to N⁢log⁡N, which compares to an expense of N2 for the DFT algorithm. A very significant factor of improvement.
The FourierTransform and InverseFourierTransform commands provide the padding option for padding the input data with zeros to obtain a better data length, but this is not always desirable, as it distorts the high frequency portion of the transform (hence, it is provided as an option).
In addition to time efficiency, storage space can also be a concern when dealing with large data Arrays. The primary consideration here is that the input data should be provided in the form of a single Array (or Vector or Matrix) with datatype=complex8, as this is the datatype used in the implementation of the algorithms, and inplace=true should be specified. If a different datatype is used (for example, two float8 Arrays, or a single complex Array), the input is converted to complex8 datatype, which effectively doubles the required storage for the input data alone.
For a composite data length, the minstorage FFT algorithm requires the least additional storage to compute the Fourier transform. Assuming the work is being done in-place (via inplace=true), in addition to the input data (which is also the output data) it requires additional complex storage of twice the size of the largest prime in the data length.
In contrast (again for a composite data length) the mintime algorithm requires additional complex storage that is equal to the maximum size of any single transform being computed (for a single 1-D transform this is simply the data length) or twice the size of the largest prime in the data length, whichever is larger.
This can become very significant for large data sets.
Note: In general, the run time of the mintime algorithm is always better or equal to the run time of the minstorage algorithm.
For a prime data length, the DFT algorithm requires the least additional storage, needing only 1.5 times the size of the input data length (which is 25% less than the mintime and minstorage algorithm if the data length is prime). But if a transform is of prime length, and sufficiently large so that the memory usage becomes an issue, it would take a very long time to compute the DFT.
As a final note for economy of space, if padding is to be used, but space is a concern, then the padding should be done manually by allocating the data array to be the desired padded size, and populating only the first entries with data values, setting all remaining entries to zero (that is, doing the padding manually).
The following guidelines should be followed for time and memory efficient use of the FourierTransform and InverseFourierTransform commands.
Always use a highly composite number of data points (optimally all factors should be between 2-5).
If the input data length is not highly composite, then manually perform the padding on creation of the input data Array.
Always use datatype=complex8, and use the inplace operation via inplace=true.
If time is a concern, then use algorithm=mintime, but if memory usage is more of a concern, then use algorithm=minstorage.
Digits ≔ 15:
Consider a problem with data length of 1 million. Storage of the complex data alone requires 16 MB. Use of the minstorage algorithm requires virtually no additional data be allocated (the largest prime factor is 5), while in this case use of the (default) mintime algorithm requires allocation of an additional 16 MB.
Note: Use ArrayTools and evalhf to more rapidly construct the input dataset (which is iN+I⁢sin⁡10⁢iN for i=1..N).
N ≔ 1000000:
Z ≔ Vector⁡N,datatype=complex8:
Za ≔ ArrayTools:-ComplexAsFloat⁡Z:
fill := proc(N,Za) local i;
for i to N do
Za[2*i-1] := sqrt(i/N);
Za[2*i] := sin(10*i/N);
And at this point, Maple's memory usage is at ~ 17MB - the size of the data Array.
The transform with 'minstorage':
tt ≔ time⁡:
and Maple's memory usage has stayed roughly constant.
The default 'mintime' transform (after resetting Z),
which for this problem takes noticeably less time, but increases Maple's memory usage to ~ 32 MB (double).
Download Help Document
What kind of issue would you like to report? (Optional)