离散傅里叶变换实现

Introduction

The discrete Fourier transform (DFT) is a basic yet very versatile algorithm for digital signal processing (DSP). This article will walk through the steps to implement the algorithm from scratch. It also provides the final resulting code in multiple programming languages.

The DFT overall is a function that maps a vector of \(n\) complex numbers to another vector of \(n\) complex numbers. Using 0-based indexing, let \(x(t)\) denote the \(t^\text{th}\) element of the input vector and let \(X(k)\) denote the \(k^\text{th}\) element of the output vector. Then the basic DFT is given by the following formula:

\( X(k) = \displaystyle\sum_{t=0}^{n-1} x(t) e^{-2 π i t k / n} \).

The interpretation is that the vector \(x\) represents the signal level at various points in time, and the vector \(X\) represents the signal level at various frequencies. What the formula says is that the signal level at frequency \(k\) is equal to the sum of {the signal level at each time \(t\) multiplied by a complex exponential}.

Coding walkthrough

With a basic knowledge of summation notation, complex numbers, and computer programming, it is straightforward to convert the description above into computer code. The purpose of this article is to guide you through the translation process step by step. We will use the Java programming language to illustrate.

Skeleton structure

The first part of the description says that the DFT takes an input vector of \(n\) complex numbers and calculates an output vector of \(n\) complex numbers. Since Java does not have a native complex number type, we will manually emulate a complex number with a pair of real numbers. A vector is a sequence of numbers, which can be represented by an array. Instead of returning the output arrays, we will have them passed in by reference as arguments. Let’s write a skeleton method:

void dft(double[] inreal , double[] inimag,
         double[] outreal, double[] outimag) {
    // Assume all 4 arrays have the same length
    int n = inreal.length;
    // Incomplete: Method body
}

Next, write the outer loop to assign a value to each output element:

void dft(double[] inreal , double[] inimag,
         double[] outreal, double[] outimag) {
    int n = inreal.length;
    for (int k = 0; k < n; k++) {  // For each output element
        outreal[k] = ?;  // Incomplete
        outimag[k] = ?;  // Incomplete
    }
}

Summation

Summation notation, while it might look intimidating, is actually easy to understand. The general form of a finite sum just means this:

\( \displaystyle\sum_{j=a}^{b} f(j) = f(a) + f(a+1) + \cdots + f(b-1) + f(b) \).

See how we substituted the values that \(j\) takes on? In imperative-style code, it looks like this:

double sum = 0;
for (int j = a; j <= b; j++) {
    sum += f(j);
}
// The value of 'sum' now has the desired answer

Complex arithmetic

The addition of complex numbers is easy:

\( (a+bi) + (c+di) = (a+c) + (b+d)i \).

The multiplication of complex numbers is slightly harder, using the distributive law and the identity \(i^2 = -1\):

\( (a+bi)(c+di) = ac + ad\,i + bc\,i - bd = (ac-bd) + (ad+bc)i \).

Euler’s formula tells us that \( e^{xi} = \cos x + i \sin x \), for any real number \(x\). Furthermore, cosine is an even function so \( \cos(-x) = \cos x \), and sine is an odd function so \( \sin(-x) = -(\sin x) \). By substitution:

\( e^{-2 π i t k / n} = e^{(-2 π t k / n) i} = \cos \left( -2 π \frac{t k}{n} \right) + i \sin \left( -2 π \frac{t k}{n} \right) = \cos \left( 2 π \frac{t k}{n} \right) - i \sin \left( 2 π \frac{t k}{n} \right) \).

Let \(\text{Re}(x)\) be the real part of \(x\) and let \(\text{Im}(x)\) be the imaginary part of \(x\). By definition, \(x = \text{Re}(x) + i \, \text{Im}(x)\). Therefore:

\( x(t) \, e^{-2 π i t k / n} = \left[ \text{Re}(x(t)) + i \, \text{Im}(x(t)) \right] \left[ \cos \left( 2 π \frac{t k}{n}\right) - i \sin \left( 2 π \frac{t k}{n} \right) \right] \).

Expand the complex multiplication to give:

\( = \left[ \text{Re}(x(t)) \cos \left( 2 π \frac{t k}{n} \right) + \text{Im}(x(t)) \sin \left( 2 π \frac{t k}{n} \right) \right] + i \left[ -\text{Re}(x(t)) \sin \left( 2 π \frac{t k}{n} \right) + \text{Im}(x(t)) \cos \left( 2 π \frac{t k}{n} \right) \right] \).

So each term in the summation has this code for the real and imaginary parts:

double angle = 2 * Math.PI * t * k / n;
double real =  inreal[t] * Math.cos(angle) + inimag[t] * Math.sin(angle);
double imag = -inreal[t] * Math.sin(angle) + inimag[t] * Math.cos(angle);

Putting it all together

Merge the code for each term of the sum into the code for the summation, and we’re done:

static void dft(double[] inreal , double[] inimag,
                double[] outreal, double[] outimag) {
    int n = inreal.length;
    for (int k = 0; k < n; k++) {  // For each output element
        double sumreal = 0;
        double sumimag = 0;
        for (int t = 0; t < n; t++) {  // For each input element
            double angle = 2 * Math.PI * t * k / n;
            sumreal +=  inreal[t] * Math.cos(angle) + inimag[t] * Math.sin(angle);
            sumimag += -inreal[t] * Math.sin(angle) + inimag[t] * Math.cos(angle);
        }
        outreal[k] = sumreal;
        outimag[k] = sumimag;
    }
}

Source code 

The end product of this DFT tutorial is available in multiple languages for your convenience, and the code is placed in the public domain:

• Java: Dft.java
• C: dft.c
• C++: Dft.cpp
• C#: Dft.cs
• JavaScript: dft.js
• TypeScript: dft.ts
• Python: dft.py
• Rust: dft.rs
• MATLAB: dft.m

Notes:

• Python, C, C++, C#, and MATLAB have built-in support for complex numbers. This feature makes our job easier and the resulting DFT implementation much simpler.

• Each implementation respects the naming convention, formatting style, and code idioms of its own language – they don’t necessarily imitate the tutorial’s Java implementation as closely as possible.

• The flavor of DFT described here does not include normalization/scaling for the sake of simplicity. There are numerous choices in how to scale the DFT – such as scaling only the forward transform by \(1/n\), scaling both the forward and inverse transforms by \(1/\sqrt{n}\), scaling the precomputed trigonometric tables, etc. – with no clear winner.

• The code on this page is a correct but naive DFT algorithm with a slow \(Θ(n^2)\) running time. A much faster algorithm with \(Θ(n \log n)\) run time is what gets used in the real world. See my page Free small FFT in multiple languages for an implementation of such.

More info

• Wikipedia: Discrete Fourier transform
• MathWorld: Discrete Fourier Transform
• The Scientist and Engineer’s Guide to Digital Signal Processing: Chapter 8: The Discrete Fourier Transform