13. Fast Fourier transform (FFT)

Note – The example worksheets are available for download. See 1. How to obtain and use example worksheets.

The discrete Fourier transform (DFT) computes the Fourier transform of the discrete data numerically. The inverse transform of DFT is called as the inverse discrete Fourier transform (IDFT). The Fast Fourier Transform (FFT) is the algorithm to calculate DFT/IDFT with a small number of arithmetic operations.

Examples of FFT computations by XLPack are shown below.

Example 1

The following signal is the mixture of the 400 Hz sine wave and the 1400 Hz sine wave.

We will compute DFT of the data of the above wave of 1/10000 second interval (refer to the reference manual on the data format of the transformed data), and then confirm the consistency of the transform by comparing the IDT of the DFT with original data.

WRfft1f worksheet function is used to compute DFT. The following is the example of DFT of 50 data points.

When DFT values are computed, IDT of those values are computed by WRfft1b worksheet function.

It has been confirmed as below that the original data are correctly restored.

Example 2 (Power spectrum)

One of the typical applications of FFT is the power spectrum analysis.

Using the DFT result of N data f_j, the power at frequency f_k is computed as follows. Where f_k = (k/N) f_s , f_s is the sampling frequency.

  
  P0 = P(f0) = N |F0|2 = N ak2
  Pk = P(fk) = N (|Fk|2 + |FN-k|2) = (N/2)(ak2 + bk2)  (k=1 to N/2-1)
  PN/2 = P(fN/2) = N |FN/2|2 = N aN/22  (if N is the even number)

 
The total energy is computed as the sum of the power over all time.

  
  E = Σ Pk  (Σ for k=0 to N/2)

 
Let’s compute the power spectrum of the DFT values of the above example.

The following VBA function is used to compute the power.

Function PsdR(N As Long, VR As Variant) As Variant
    Dim R() As Double, P() As Double, I As Long
    ReDim R(N - 1), P(N / 2, 0)
    For I = 1 To N
        R(I - 1) = VR(I)
    Next
    P(0, 0) = N * R(0) ^ 2
    For I = 1 To N - 3 Step 2
        P((I + 1) / 2, 0) = (N / 2) * (R(I) ^ 2 + R(I + 1) ^ 2)
    Next
    P(N / 2, 0) = N * R(N - 1) ^ 2
    PsdR = P
End Function

The input data is set to this function as below.

A plot of power versus frequency is as follows (this is called the periodogram).

We see that the signal consists of 400 Hz and 1400 Hz sine waves.