Summary
The #Recommended method’s combination of the DFT and FFT is in fftdftr.zip.
A more general DFT code that works with complex data is described in DFT code.doc .htm. The codes described in this document assume that the input data is real.
The codes described in this section are designed for experimental data that starts with d[0]=0 and ends at another value. This data is reflected (#Data Extension) inside the Fourier transform codes when requested. The data transformed then has the property that it can be zero padded to any size. A truncated Fourier transform is a least squares fit to this data. The #Recommended method is to find a Df appropriate for the fast Fourier transform, use the slow transform with the experimental data and data extension with a limited number of Fourier coefficients NF. Then back transform this data using the FFT to find the data and optionally its derivative with a different Dt.
The DFTR code (dftras.zip) in the forward direction LF=.TRUE. reads d[i] values 0 to NT-1 from the file namet. These are used to calculate complex transform values D[m] from 0 to NF. If LEXT is also true, a simple expansion of the data is used to make the data set appropriate for zero padding. There is almost no internal storage. The code can be competitive with the Fast Fourier Transform owing to its small size and direct use of the external files. It is especially useful when NT and NF are both much smaller than N.
The FFTR code differs from the above in that it has an internal complex array of size N=2M. This array allows the use of the Fast Fourier Transform algorithm which changes the total number of operations from NT´NF to N´log2(N).
Discrete Fourier Transform
The Discrete Fourier Transform is defined in Discrete Fourier Transform.doc .htm (3.7) and (3.8)
Equation (1.2) defines D[m] for all m. The back sum is
Equation (1.3) defines d[i] for all i. It is periodic with period N so that d[i+nN] = d[i]. If
With this value of N, (1.2) and (1.3) are periodic with period N so that for n any integer n D[m+nN]
For real data the sum in (1.3) is shown in RealityCondition.doc .htm to be equivalent to
Odd values of N/2 are truncated. That is 9/2 = 4, so that for N=9, the above sum is from -4 to 4. Because N-m has its maximum at N/2, the largest frequencies in (1.5) are at the ends of the summation range. There is no imaginary part to D[0].
Figure 1 Typical data set sin2(t) 15 points from from 0 to 1.
The value of N given in (1.4) of the evenly spaced data over N data points requires
(1.6)
Figure 2 Forward transform, followed by back transform with the same data.
Data Extension
Choose a Dt such that
N’ = 1/(DtDf) > N
Then evaluate (1.3). The data interpolates with a few bumps.
Figure 3 Df = 1/(N´Dt), Dt à Dt /10 black = data, red = interpolation.
N’ frequency values were summed evaluating each element of (1.3). N’-N of these were zero values making for an immediate transform. The periodicity of the results are determined by T=1/Df which is unchanged in this sum. When the new Dt is a multiple of the old, all of the points corresponding to the old Dt are on the line representing the original data. The oscillations are due to the sudden discontinuity in the data where the end of the range periodically continues to the beginning of the range. This can be corrected.
If the data is extended by extending the data points from NT-1 to 1 (d[NT] = d[NT-2], …,d[2NT-3]=d[1], the result is shown in figure 1.
Figure 4 Original data - black dots - extended data red line.
This data starts at zero and ends at zero. It is appropriate for ZeroPadding.doc .htm.[1] If the data is assumed zero beyond 2NT-3, equation (1.2) becomes
In the second sum let i’ = 2N-2-i
Code evaluating (2.2) using the sine and cosine directly is in dft0.zip
Change the order of the sum and factor out the terms in the exponential that are not varying
(2.3)
Details about the fastest way to find the exponentials are in Advancing the sin.doc .htm.
The second set of code in DFT code.doc .htm was used to advance the complex exponentials directly once the beginning frequency (forward) or time (backward) is found. This reduces the calculation of these to four multiplications and two additions. The resulting code is in dftras.zip.
Reflection
Reflection can be made for real data with a first point equal to zero. The NT used in this section becomes NText if the data has been extended as discussed in #Data Extension.
Reflect the data such that
(3.1)
Then (1.2) with d[0]=0 as in the #Assumed_data_form so that
It might seem that the number of coefficients has been halved by this reflection owing to the absence of real(D[m]). The Time parameter, however, has been doubled making the remainder of the D’s spaced half the distance apart in frequency space.
Smoothing
Suppose that we replace the data by
Where
The second line in (4.2) is fundamental. It differs from the first when the trap rule integral taking d and h to integrals is not well approximated by the sum. The sums always are over a single period of the data or of the extended data. If d or h goes to zero with all derivatives 0 at t = Time/2, the two lines in (4.2) will be the same. Details about this are discussed in
DiscreteConvolution.doc .htm and Nyquist.doc .htm
A Gaussian H(f) is optimal for many transforms gabor\Gabor Transforms.doc. Details about the analytic transform are in TransformOfGaussian.doc .htm. Taking f0 = 0
(4.3)
Defining
(4.4)
Equation (4.5) is the second approximation of a delta function given in ..\definitions\Deltafun.doc .htm. For small values of wt (4.1) makes dsm(t) = d(t) and H(f) = 1. The frequencies extend from –Nf/Time to Nf/Time. Thus the end values of H are
(4.6)
Figure 6 2000 data points – extended, reflected, zero padded and smoothed. wt = 0.1.
Note that the sin2 used as a sample function does not naturally reflect. The code in expreflsm.zip does this using the dft only.
In the frequency domain, the smoothing is
(4.7)
As in the time domain, there are theoretical reasons to believe that a Gaussian H is ideal
A Gaussian H(f) is optimal for many transforms gabor\Gabor Transforms.doc. Details about the analytic transform are in TransformOfGaussian.doc .htm. Taking f0 = 0
(4.8)
(4.9)
Or
(4.10)
At the end of the transform
(4.11)
The FFT
The FFT factors the Fourier sum so that operations combine to speed up the formation of the transform.
(5.1)
Powers of 2 are not the only possible factorization, but they are convenient. If there are fewer than NP2 values, the rest need to be explicitly added as zeroes. The transform requires NP2 log2(NP2) times a small factor involving table look up, 6 multiplications and two additions. It is fast.
The DFT used above uses NT points extended and reflected to be the same as 4´NT points. The data is normally smoothed and truncated in frequency space so that the total number of operations is a small factor times NT´NF.
DFT : NT ´ NF
FFT : NP2 ´ 4 ´ log2(NP2)
The second term is about 40 for NT ~ 1000. Thus for 1000 data points, represented by 160 coefficients, the FFT is not faster than the DFT. By 10000 data points and 1600 coefficients, there is no question that the FFT is 4 times faster. This continues until at some millions of data points the memory overflows leading to a factor of 1000 slowdown.
Start in the time domain with a given value of Dt and 4´NT data points. Make NP2 ³ 4´NT a power of 2. Then make a forward transform of zero padded data with
This transform gives N = NF complex frequency values that are largely redundant since D(-m) = D*(m). If the data has been properly extended the set of D’s converge fast enough that it can be truncated at NF < N/2.
To transform to a different Dt hat it is necessary to find a different N so that
(5.3)
The number of points Nhat can be twice the original so that the new data points are halfway between the old etc, but arbitrary spacing is best done with the DFT.
Recommended method
The recommended method for fixed Dt and Dt hat is to first find a reasonably large Nhat = 2M such that Nhat ³ 2NT-2 (use Lext = true). Then use (5.2) to find
(5.4)
The DFT code in the section #Data Extension should next be used with the original Dt, NT, Df and NF < Nhat/2 to find D[m]. Then this can be used with
The subroutine FFTRALL has parameters (LF,LEXT,LDER,NP2,NT,NF,TIME,NAMET,NAMEF)
These are F,F,T,
If the number of data points is ~ 1000, the extended data going back to zero requires 2000 data points. A set of 100 complex frequency values represents one fitting parameter for every 10 data points. Thus the DFT code requires on the order of 6´2000´200=2.4´106 operations for the transform that finds 100 complex D[m] values. .
The FFT with a such that 2000 à 2048 requires 6 ´ 2048 ´ 11 operations. The 6 is the 4 required to advance the sine, plus one for rearranging the array plus one for finding the elements in the rearranged array. The total is 0.135´106 operations for each transform. This is a saving of about a factor of 20 over repeated use of the DFT code.
Code for using the FFT including the calculation of the derivative is in FFTRALL.zip
Derivative
The form of d shown in (1.3) is deceptive for taking derivatives due to the fact that terms from m = N/2 to N do not oscillate with their apparent frequency. For example m=N-1 is through periodic repeats exactly equivalent to m = -1. This does not happen for the d(t) form given by (1.5)
(5.5)
This is the back transform of
(5.6)
Figure 7 Sin2 and its derivative using 2000 Poisson distributed data points, data extension and data reflection. (method described above). Data reflection makes for the obvious error at the reflected value near 10.5.
The code described above, including a Pascal version, is in fftdftr.zip. The function sin2(t) is generated with Dt. The resulting points by design begin at 0 and extend to the peak of the data hump in Figure 3. The value of Df needed for an appropriately extended set of points N2P is calculated. The DFT code with LEXT = true is used to find the transform of the full hump in figure 3 at a set of data points appropriate for an FFT. The code FFTRAL calculates the back transform of N2P points which is plotted above. Then a second call is made with LDER = true to calculate the derivative shown.
Code using the DFT is in dddt.zip, the derivative is an option in #FFTRALL above.
Integral
The function d(ti) is given by (1.5). this becomes d(t), with the appropriate extensions and reflections that allow the sums to be zero padded to any region.
If this form is used in its present form, the first term in (5.7) will integrate to t introduces a discontinuity in the periodic function representing the integral. Use #Reflection to remove the D[0] term so that d is represented by the continuous version of (3.3).
(5.8)
Integrate term by term
(5.9)
This is a new series
(5.10)
(5.11)
This has a term at zero and will require work to integrate further, but we already knew that integration would have difficulties.
Figure 3 Ihtegral of sin2(pt/2 ) numerical - black dots – 150 points with random component 0.01*Gauss*sqrt(f), extended, reflected and smoothed with 0.4.