你還弄不懂的傅里葉變換，神經(jīng)網(wǎng)絡(luò)只用了30多行代碼就學(xué)會(huì)了

明敏 發(fā)自 凹非寺
量子位 報(bào)道 | 公眾號(hào) QbitAI

在我們的生活中，大到天體觀測、小到MP3播放器上的頻譜，沒有傅里葉變換都無法實(shí)現(xiàn)。

通俗來講，離散傅里葉變換（DFT）就是把一串復(fù)雜波形中分成不同頻率成分。

比如聲音，如果用聲波記錄儀顯示聲音的話，其實(shí)生活中絕大部分聲音都是非常復(fù)雜、甚至雜亂無章的。

而通過傅里葉變換，就能把這些雜亂的聲波轉(zhuǎn)化為正弦波，也就是我們平常看到的音樂頻譜圖的樣子。

不過在實(shí)際計(jì)算中，這個(gè)過程其實(shí)非常復(fù)雜。

如果把聲波視作一個(gè)連續(xù)函數(shù)，它可以唯一表示為一堆三角函數(shù)相疊加。不過在疊加過程中，每個(gè)三角函數(shù)的加權(quán)系數(shù)不同，有的要加高一些、有的要壓低一些，有的甚至不加。

傅里葉變換要找到這些三角函數(shù)以及它們各自的權(quán)重。

這不就巧了，這種找啊找的過程，像極了神經(jīng)網(wǎng)絡(luò)。

神經(jīng)網(wǎng)絡(luò)的本質(zhì)其實(shí)就是逼近一個(gè)函數(shù)。

那豈不是可以用訓(xùn)練神經(jīng)網(wǎng)絡(luò)的方式來搞定傅里葉變換？

這還真的可行，并且最近有人在網(wǎng)上發(fā)布了自己訓(xùn)練的過程和結(jié)果。

DFT=神經(jīng)網(wǎng)絡(luò)

該怎么訓(xùn)練神經(jīng)網(wǎng)絡(luò)呢？這位網(wǎng)友給出的思路是這樣的：

首先要把離散傅里葉變換（DFT）看作是一個(gè)人工神經(jīng)網(wǎng)絡(luò)，這是一個(gè)單層網(wǎng)絡(luò)，沒有bias、沒有激活函數(shù)，并且對(duì)于權(quán)重有特定的值。它輸出節(jié)點(diǎn)的數(shù)量等于傅里葉變換計(jì)算后頻率的數(shù)量。

具體方法如下：

這是一個(gè)DFT：

• k表示每N個(gè)樣本的循環(huán)次數(shù)；
• N表示信號(hào)的長度；
• 表示信號(hào)在樣本n處的值。

一個(gè)信號(hào)可以表示為所有正弦信號(hào)的和。

yk是一個(gè)復(fù)值，它給出了信號(hào)x中頻率為k的正弦信號(hào)的信息；從yk我們可以計(jì)算正弦的振幅和相位。

換成矩陣式，它就變成了這樣：

這里給出了特定值k的傅里葉值。

不過通常情況下，我們要計(jì)算全頻譜，即k從[0,1,…N-1]的值，這可以用一個(gè)矩陣來表示（k按列遞增，n按行遞增）：

簡化后得到：

看到這里應(yīng)該還很熟悉，因?yàn)樗且粋€(gè)沒有bias和激活函數(shù)的神經(jīng)網(wǎng)絡(luò)層。

指數(shù)矩陣包含權(quán)值，可以稱之為復(fù)合傅里葉權(quán)值（Complex Fourier weights），通常情況下我們并不知道神經(jīng)網(wǎng)絡(luò)的權(quán)重，不過在這里可以。

• 不用復(fù)數(shù)

通常我們也不會(huì)在神經(jīng)網(wǎng)絡(luò)中使用復(fù)數(shù)，為了適應(yīng)這種情況，就需要把矩陣的大小翻倍，使其左邊部分包含實(shí)數(shù)，右邊部分包含虛數(shù)。

將

帶入DFT，可以得到：

然后用實(shí)部（cos形式）來表示矩陣的左半部分，用虛部（sin形式）來表示矩陣的右半部分：

簡化后可以得到：

將

稱為傅里葉權(quán)重；

需要注意的是，y^和y實(shí)際上包含相同的信息，但是y^

不使用復(fù)數(shù)，所以它的長度是y的兩倍。

換句話說，我們可以用

或

表示振幅和相位，但是我們通常會(huì)使用

現(xiàn)在，就可以將傅里葉層加到網(wǎng)絡(luò)中了。

用傅里葉權(quán)重計(jì)算傅里葉變換

現(xiàn)在就可以用神經(jīng)網(wǎng)絡(luò)來實(shí)現(xiàn)

，并用快速傅里葉變換（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

得到的這個(gè)微小誤差值可以證明，計(jì)算的結(jié)果是我們想要的。

• 另一種方法是重構(gòu)信號(hào)：

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后從正弦信號(hào)重建的信號(hào)和原始信號(hào)能夠很好地重合。

通過梯度下降學(xué)習(xí)傅里葉變換

現(xiàn)在就到了讓神經(jīng)網(wǎng)絡(luò)真正來學(xué)習(xí)的部分，這一步就不需要向之前那樣預(yù)先計(jì)算權(quán)重值了。

首先，要用FFT來訓(xùn)練神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)離散傅里葉變換：

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.GradientTape() 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

得出結(jié)果如上，這證實(shí)了神經(jīng)網(wǎng)絡(luò)確實(shí)能夠?qū)W習(xí)離散傅里葉變換。

訓(xùn)練網(wǎng)絡(luò)學(xué)習(xí)DFT

除了用快速傅里葉變化的方法，還可以通過網(wǎng)絡(luò)來重建輸入信號(hào)來學(xué)習(xí)DFT。(類似于autoencoders自編碼器)。

自編碼器（autoencoder, AE）是一類在半監(jiān)督學(xué)習(xí)和非監(jiān)督學(xué)習(xí)中使用的人工神經(jīng)網(wǎng)絡(luò)（Artificial Neural Networks, ANNs），其功能是通過將輸入信息作為學(xué)習(xí)目標(biāo)，對(duì)輸入信息進(jìn)行表征學(xué)習(xí)（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.GradientTape() 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

作者用這一模型進(jìn)行了很多測試，最后得到的權(quán)重不像上面的例子中那樣接近傅里葉權(quán)值，但是可以看到重建的信號(hào)是一致的。

換成輸入振幅和相位試試看呢。

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.GradientTape() 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

可以看到，重建信號(hào)再次一致；

不過，和此前一樣，輸入振幅和相位最終得到的權(quán)值也不完全等同于傅里葉權(quán)值(但非常接近)。

由此可以得出結(jié)論，雖然最后得到的權(quán)重還不是最準(zhǔn)確的，但是也能夠獲得局部的最優(yōu)解。

這樣一來，神經(jīng)網(wǎng)絡(luò)就學(xué)會(huì)了傅里葉變換！

• 值得一提的是，這個(gè)方法目前還有疑問存在：

首先，它沒有解釋計(jì)算出的權(quán)值和真正的傅里葉權(quán)值相差多少；

而且，也沒有說明將傅里葉層放到模型中能帶來哪些益處。

原文鏈接：
https://sidsite.com/posts/fourier-nets/