13. Fast Fourier Transform (FFT)

13.1 Overview

Note – In the definitions of Fourier Transform/Inverse Transform, DFT/IDFT, etc., the normalization factor (such as 1/N) may vary depending on the literature. Please adjust accordingly.

13.1.1 Fourier Transform

When \(f(x)\) is a complex function, the following is called the Fourier Transform:
\[
F(u) = (1/2\pi)\int_{-\infty}^{\infty}f(x)e^{-iux}dx
\] The inverse of this is called the Fourier Inverse Transform:
\[
f(x) = \int_{-\infty}^{\infty}F(u)e^{iux}du
\]

13.1.2 Discrete Fourier Transform (DFT)

The Fourier Transform (its discrete approximation) can be numerically computed using a finite set of input function values (discrete values).

We use the function values \(f(x_j)\), denoted as \(f_j\), at \(N\) equally spaced sample points \(x_j = 2{\pi}j/N (j = 0 \sim N – 1)\).

Assuming \(f(x)\) is a periodic function with period \(2\pi\), and integrating \(F(u)\) over the interval \([0, 2\pi]\) using the composite trapezoidal rule, we obtain:
\[
F(u) = (1/N)\sum_{j=0}^{N-1}f_je^{-iux_j}
\] Setting \(u_k = 2{\pi}k\) and \(F_k = F(u_k)\), we get:
\[
F_k = (1/N)\sum_{j=0}^{N-1}f_je^{-i2{\pi}jk/N} \space (k = 0 \sim N – 1)
\] This is called the Discrete Fourier Transform (DFT), which approximates the Fourier Transform.

Similarly, the Inverse Discrete Fourier Transform (IDFT) is expressed as:
\[
f_j = \sum_{k=0}^{N-1}F_ke^{i2{\pi}kj/N} \space (j = 0 \sim N – 1)
\]

13.1.3 Discrete Fourier Transform for Real Numbers

When the input function values are real numbers, the relation \(F_{N-k} = \overline{F_k}\) holds. (\(\overline{F_k}\) represents the conjugate complex of \(F_k\).) Therefore, calculations only need to be performed for \(k = 0 \sim N/2\). And \(F_0\) and \(F_{N/2}\) are real numbers.

Usually, we substitute \(a_k = 2Re(F_k)\) and \(b_k = -2Im(F_k)\) (where \(Re\) represents the real part and \(Im\) represents the imaginary part), and express the computations using real-number operations as follows:
\[
\begin{align}
& a_0 = (1/N)\sum_{j=0}^{N-1}f_j \\
& a_k = (2/N)\sum_{j=0}^{N-1}f_jcos(2{\pi}jk/N) \space (k = 1 \sim N2) \\
& a_{N/2} = (1/N)\sum_{j=0}^{N-1}(-1)^jf_j \\
& b_k = (2/N)\sum_{j=0}^{N-1}f_jsin(2{\pi}jk/N) \space (k = 1 \sim N2) \\
\end{align}
\] where, \(N2 = N/2 – 1\) (if \(N\) is even), \(N2 = (N – 1)/2\) (if \(N\) is odd).

The IDFT is given as:
\[
\begin{align}
& For even N:
f_j = a_0 + \sum_{k=1}^{N/2}(a_kcos(2{\pi}kj/N) + b_ksin(2{\pi}kj/N)) + (-1)^ka_N/2 \space (j = 0 \sim N – 1) \\
& For odd N:
f_j = a_0 + \sum_{k=1}^{(N-1)/2}(a_kcos(2{\pi}kj/N) + b_ksin(2{\pi}kj/N)) \space (j = 0 \sim N – 1) \\
\end{align}
\]

13.2 Fast Fourier Transform (FFT)

Fast Fourier Transform (FFT) is a numerical computation method for obtaining the DFT and the IDFT with fewer operations.

If the DFT is computed directly according to its definition, \(N\) calculations must be performed within the summation for each of the \(N\) values of \(F_k\), leading to an overall computational complexity of \(N^2\) operations.

However, if \(N\) can be expressed as a product of small prime powers, the FFT algorithm can significantly reduce the computational load. For instance, if \(N = 2^n\) (i.e., a power of 2), the computation can be performed in approximately \(Nlog_2 N\) operations. As an example, when \(N = 1024\), a standard computation would require roughly one million operations, whereas FFT reduces it to about ten thousand, making the process 100 times faster.

13.3 How to Use Fast Fourier Transform (FFT) with XLPack

Using the VBA subroutines Rfft1f, Rfft1b, Rfft1i, or the worksheet functions WRfft1f, WRfft1b, you can compute the one-dimensional real Fourier transform and inverse Fourier transform.

Example Problem

When a 2.5 kHz carrier wave is amplitude-modulated (AM) with a combination of a 100 Hz sine wave and a 200 Hz sine wave, the resulting waveform appears as follows:

Using this data, we will compute the DFT and the IDFT.

13.3.1 Using Worksheet Functions for Computation

To compute the DFT, we prepare 1000 data points sampled at intervals of 1/20000 seconds (i.e., at a 20 kHz sampling rate). For details on how the transformation results are stored, refer to the reference manual. We then compute the IDFT to confirm that the transformation has been correctly performed.

To obtain the DFT, we use the WRfft1f worksheet function.

Once the DFT has been calculated, we apply the WRfft1b worksheet function to perform the IDFT.

As shown below, the original data has been successfully restored.

13.3.2 Power Spectrum

One of applications of FFT is power spectrum analysis, which is equivalent to observing the frequency distribution of an input signal using a spectrum analyzer.

From the DFT output of \(N\) data points \(f_j\), the power \(P_k\) at frequency \(\omega_k\) can be determined as follows, where \(\omega_k = (k/N)\omega_s\) (\(\omega_s\) is the sampling frequency):
\[
\begin{align}
& P_0 = P(f_0) = N |F_0|^2 = N a_k^2 \\
& P_k = P(f_k) = N (|F_k|^2 + |F_{N-k}|^2) = (N/2)(a_k^2 + b_k^2) \space (k = 1 \sim N/2 – 1) \\
& P_N/2 = P(f_N/2) = N |F_N/2|^2 = N a_{N/2}^2 \space (N が偶数の場合) \\
\end{align}
\] The total energy \(E\) can be obtained as the sum of these values:
\[
E = \sum_{k=0}^{N/2} P_k
\] Now, the dataset used in the previous example was AM-modulated waveform data sampled at 20 kHz. Based on the DFT results, we will compute the power spectrum.

Since performing these calculations directly on a worksheet can be complex, the following VBA function is prepared to compute power values.

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

Using this function, we compute the power spectrum.

The DFT results are specified as input.

The power spectrum is obtained. The generated frequency distribution graph appears as follows (this is called a periodogram).

It is confirmed that the AM-modulated wave contains signals around the 2.5 kHz carrier wave at ±100 Hz and ±200 Hz.