从零开始：用神经网络实现离散傅里叶变换（DFT）

在我们日常生活中，从天体观测到MP3播放器上的频谱显示，都离不开傅里叶变换。本教程将深入讲解如何利用神经网络实现离散傅里叶变换（DFT），并通过训练让神经网络自主学会傅里叶变换。你会发现，DFT与神经网络之间存在着令人惊叹的关联。

1. DFT与神经网络：一个直观的类比

通俗地说，离散傅里叶变换（DFT）能把复杂波形分解成不同频率的成分。例如声音，如果用声波记录仪来显示，你会发现绝大多数声音都是极其复杂甚至杂乱无章的。而通过傅里叶变换，就能将这些杂乱声波转化为正弦波，也就是我们常见的音乐频谱图。

但在实际计算中，这个过程相当复杂。如果把声波视为一个连续函数，它可以唯一地表示为一堆三角函数的叠加。在叠加过程中，每个三角函数的加权系数各不相同，有的要加高，有的要压低，甚至有些不加。傅里叶变换的任务就是找出这些三角函数以及它们各自的权重。

巧合的是，这种“寻找”过程与神经网络的工作机制非常相似。神经网络的本质是逼近一个函数，那么是否可以通过训练神经网络的方式来解算傅里叶变换？这个设想确实可行，而且最近已经有人在网上分享了他们的训练过程和结果。

2. DFT=神经网络：数学结构完全一致

那么，如何训练神经网络呢？这位网友给出的思路是：首先将离散傅里叶变换（DFT）看作一个人工神经网络——这是一个单层网络，没有偏置项、没有激活函数，并且权重有特定值。其输出节点的数量等于傅里叶变换计算后频率的数量。

具体方法如下：

这是一个DFT表达式：

k表示每N个样本的循环次数；
N表示信号的长度；
表示信号在样本n处的值。

一个信号可以表示为所有正弦信号的叠加。

是一个复数值，它给出了信号x中频率为k的正弦信号的信息；从我们可以计算出正弦波的振幅和相位。

转换成矩阵形式就变成了这样：

这里给出了特定值k的傅里叶值。通常情况下，我们要计算全频谱——即k从[0,1,…N-1]的值，这可以用一个矩阵来表示（k按列递增，n按行递增）：

简化后得到：

看到这里应该很熟悉了，因为它正是一个没有偏置项和激活函数的神经网络层。

指数矩阵中包含权值，我们可以称之为复合傅里叶权值（Complex Fourier weights）。一般情况下我们不知道神经网络的权重，但在这里却可以预先确定。

2.1 不用复数：将权重分解为实部和虚部

通常我们不会在神经网络中使用复数，为了适应这一情况，需要将矩阵大小翻倍，使得左边部分包含实部，右边部分包含虚部。

将代入DFT，可以得到：

然后用实部（cos形式）来表示矩阵的左半部分，用虚部（sin形式）来表示矩阵的右半部分：

简化后可得：

将称为傅里叶权重；

需要注意的是，和实际上包含相同的信息，但不使用复数，所以它的长度是的两倍。

换句话说，我们可以用或来表示振幅和相位，但通常我们会使用。

小提示：在构建神经网络时，我们通常用实数张量表示复数权重，将实部和虚部分开存放。这样可以直接使用标准的深度学习框架（如TensorFlow、PyTorch）进行训练。

3. 用傅里叶权重计算傅里叶变换（验证阶段）

现在就可以将傅里叶层加入网络中了。利用傅里叶权重计算傅里叶变换，并用快速傅里叶变换（FFT）检验其正确性。

import matplotlib.pyplot as plt

y_real = y[:, :signal_length]
y_imag = y[:, signal_length:]

tvals = np.arange(signal_length).reshape([-1, 1])
freqs = np.arange(signal_length).reshape([1, -1])
arg_vals = 2 * np.pi * tvals * freqs / signal_length
sinusoids = (y_real * np.cos(arg_vals) - y_imag * np.sin(arg_vals)) / signal_length

reconstructed_signal = np.sum(sinusoids, axis=1)
print('rmse:', np.sqrt(np.mean((x - reconstructed_signal)**2)))

plt.subplot(2, 1, 1)
plt.plot(x[0,:])
plt.title('Original signal')
plt.subplot(2, 1, 2)
plt.plot(reconstructed_signal)
plt.title('Signal reconstructed from sinusoids after DFT')
plt.tight_layout()
plt.show()

rmse: 2.3243522568191728e-15

得到的这个微小误差值表明，计算结果完全符合预期。

另一种方法是重构信号：

import matplotlib.pyplot as plt

y_real = y[:, :signal_length]
y_imag = y[:, signal_length:]

tvals = np.arange(signal_length).reshape([-1, 1])
freqs = np.arange(signal_length).reshape([1, -1])
arg_vals = 2 * np.pi * tvals * freqs / signal_length
sinusoids = (y_real * np.cos(arg_vals) - y_imag * np.sin(arg_vals)) / signal_length

reconstructed_signal = np.sum(sinusoids, axis=1)
print('rmse:', np.sqrt(np.mean((x - reconstructed_signal)**2)))

plt.subplot(2, 1, 1)
plt.plot(x[0,:])
plt.title('Original signal')
plt.subplot(2, 1, 2)
plt.plot(reconstructed_signal)
plt.title('Signal reconstructed from sinusoids after DFT')
plt.tight_layout()
plt.show()

rmse: 2.3243522568191728e-15

最终可以看到，DFT后从正弦信号重建的信号与原始信号完美重合。

4. 通过梯度下降学习傅里叶变换

现在进入真正的学习环节——让神经网络自主训练，不再需要事先计算权重值。

4.1 方法一：用FFT作为监督信号

首先，利用快速傅里叶变换（FFT）作为标签来训练神经网络学习DFT：

import tensorflow as tf

signal_length = 32

# Initialise weight vector to train:
W_learned = tf.Variable(np.random.random([signal_length, 2 * signal_length]) - 0.5)

# Expected weights, for comparison:
W_expected = create_fourier_weights(signal_length)

losses = []
rmses = []

for i in range(1000):
    # Generate a random signal each iteration:
    x = np.random.random([1, signal_length]) - 0.5

    # Compute the expected result using the FFT:
    fft = np.fft.fft(x)
    y_true = np.hstack([fft.real, fft.imag])

    with tf.GradientTa pe() as tape:
        y_pred = tf.matmul(x, W_learned)
        loss = tf.reduce_sum(tf.square(y_pred - y_true))

    # Train weights, via gradient descent:
    W_gradient = tape.gradient(loss, W_learned)
    W_learned = tf.Variable(W_learned - 0.1 * W_gradient)

    losses.append(loss)
    rmses.append(np.sqrt(np.mean((W_learned - W_expected)**2)))

# Final loss value 1.6738563548424711e-09
# Final weights' rmse value 3.1525832404710523e-06

上述结果证实，神经网络确实能够通过学习掌握离散傅里叶变换。

4.2 方法二：通过重建输入信号来学习DFT（自编码器方式）

除了使用FFT作为监督信号，还可以让网络通过重建输入信号来学习DFT，这种方式类似于自编码器（autoencoder）。

自编码器（autoencoder, AE）是一类在半监督学习和非监督学习中常用的人工神经网络（Artificial Neural Networks, ANNs），其目标是将输入信息作为学习目标，对输入信息进行表征学习（representation learning）。

W_learned = tf.Variable(np.random.random([signal_length, 2 * signal_length]) - 0.5)

tvals = np.arange(signal_length).reshape([-1, 1])
freqs = np.arange(signal_length).reshape([1, -1])
arg_vals = 2 * np.pi * tvals * freqs / signal_length
cos_vals = tf.cos(arg_vals) / signal_length
sin_vals = tf.sin(arg_vals) / signal_length

losses = []
rmses = []

for i in range(10000):
    x = np.random.random([1, signal_length]) - 0.5

    with tf.GradientTa pe() as tape:
        y_pred = tf.matmul(x, W_learned)
        y_real = y_pred[:, 0:signal_length]
        y_imag = y_pred[:, signal_length:]

        sinusoids = y_real * cos_vals - y_imag * sin_vals
        reconstructed_signal = tf.reduce_sum(sinusoids, axis=1)

        loss = tf.reduce_sum(tf.square(x - reconstructed_signal))

    W_gradient = tape.gradient(loss, W_learned)
    W_learned = tf.Variable(W_learned - 0.5 * W_gradient)

    losses.append(loss)
    rmses.append(np.sqrt(np.mean((W_learned - W_expected)**2)))

# Final loss value 4.161919455121241e-22
# Final weights' rmse value 0.20243339269590094

作者用这一模型进行了大量测试，最终得到的权重虽然不像方法一那样接近标准的傅里叶权值，但重建的信号却与原始信号完全一致。

4.3 换成输入振幅和相位试试看

W_learned = tf.Variable(np.random.random([signal_length, 2 * signal_length]) - 0.5)

losses = []
rmses = []

for i in range(10000):
    x = np.random.random([1, signal_length]) - .5

    with tf.GradientTa pe() as tape:
        y_pred = tf.matmul(x, W_learned)
        y_real = y_pred[:, 0:signal_length]
        y_imag = y_pred[:, signal_length:]

        amplitudes = tf.sqrt(y_real**2 + y_imag**2) / signal_length
        phases = tf.atan2(y_imag, y_real)

        sinusoids = amplitudes * tf.cos(arg_vals + phases)
        reconstructed_signal = tf.reduce_sum(sinusoids, axis=1)

        loss = tf.reduce_sum(tf.square(x - reconstructed_signal))

    W_gradient = tape.gradient(loss, W_learned)
    W_learned = tf.Variable(W_learned - 0.5 * W_gradient)

    losses.append(loss)
    rmses.append(np.sqrt(np.mean((W_learned - W_expected)**2)))

# Final loss value 2.2379359316633115e-21
# Final weights' rmse value 0.2080118219691059

可以看出，重建信号再次保持一致；但与之前类似，输入振幅和相位所得到的权值并不完全等于傅里叶权值，不过非常接近。

由此可以得出结论：虽然最终得到的权重并非绝对精确，但网络依然能够找到局部最优解。这样，神经网络就成功学会了傅里叶变换！

5. 常见问题与说明

Q1：为什么DFT可以看作是一个神经网络层？
A：DFT的数学形式本质上是一个线性变换（矩阵乘法），而一个不带偏置和激活函数的全连接层也同样是矩阵乘法。DFT的权重矩阵由复指数构成，相当于一个预先确定的权重矩阵。因此，完全可以将DFT视为一个单层神经网络。

Q2：训练出来的权重和标准傅里叶权重有什么区别？
A：在使用FFT作为监督信号时（方法一），训练出的权重非常接近标准傅里叶权重（rmse约为3e-6）。而在使用自编码器方式时（方法二和方法三），训练出的权重虽然在数值上不严格等于标准权重，但重建信号的效果完全一致，说明网络找到了局部最优解。本质上，自编码器学习的是能够完美重建信号的任意线性变换，并不要求权重必须严格等于傅里叶权重，只要变换可逆即可。

小提示：如果你在实际项目中将傅里叶层嵌入神经网络，建议用预计算的傅里叶权重初始化（即方法一中的W_expected），然后允许微调，这样既能保留频率分析的先验知识，又能获得适应任务的灵活性。

6. 有待探讨的问题

值得一提的是，这个方法目前仍存在一些疑问：

• 首先，它并没有解释计算出的权值与真正的傅里叶权值之间的差异有多大；
• 其次，也没有说明将傅里叶层放入模型中能带来哪些具体益处。

编辑：黄飞