一、前言
我想認真寫好快速傅里葉變換（Fast Fourier Transform，FFT），所以這篇文章會由淺到細，由窄到寬的講解，但是傅里葉變換對於尋常人並不是很容易理解的，所以對於基礎不牢的人我會通過前言普及一下相關知識。
我們復習一下三角函數的標准式：
$$y=A\cos (\omega x+\theta )+k$$
$A$代表振幅，函數周期是$\frac{2\pi}{w}$，頻率是周期的倒數$\frac{w}{2\pi}$，$\theta $是函數初相位，$k$在信號處理中稱為直流分量。這個信號在頻域就是一條豎線。
我們再來假設有一個比較復雜的時域函數$y=f(t)$，根據傅里葉的理論，任何一個周期函數可以被分解為一系列振幅A，頻率$\omega$或初相位$\theta $正弦函數的疊加
$$y = A_1sin(\omega_1t+\theta_1) + A_2sin(\omega_2t+\theta_2) + A_3sin(\omega_3t+\theta_3)$$
該信號在頻域有三條豎線組成，而豎線圖我們把它稱為頻譜圖，大家可以通過下面的動畫了解
如圖可知，通過時域到頻域的變換，我們得到了一個從側面看的頻譜，但是這個頻譜並沒有包含時域中全部的信息。因為頻譜只代表每個正弦波對應頻率的振幅是多少，而沒有提到相位。基礎的正弦波$Asin(wt+\theta )$中，振幅，頻率，相位缺一不可，不同相位決定了波的位置，所以對於頻域分析，僅僅有頻譜（振幅譜）是不夠的，我們還需要一個相位譜。
我依稀記得高中學正弦函數的是時候，$\theta $的多少決定了正弦波向右移動多少。當然那個時候橫坐標是相位角度，而時域信號的橫坐標是時間，因此我們只需要將時間轉換為相位角度就得到了初相位。相位差則是時間差在一個周期中所占的比例
$$\theta=2\pi \frac{t}{T}$$
所以傅里葉變換可以把一個比較復雜的函數轉換為多個簡單函數的疊加，將時域（即時間域）上的信號轉變為頻域（即頻率域）上的信號，看問題的角度也從時間域轉到了頻率域，因此在時域中某些不好處理的地方，在頻域就可以較為簡單的處理，這就可以大量減少處理信號計算量。信號經過傅里葉變換后，可以得到頻域的幅度譜以及相位譜，信號的幅度譜和相位譜是信號傅里葉變換后頻譜的兩個屬性。
傅里葉用途
- 時域復雜的函數，在頻域就是幾條豎線
- 求解微分方程，傅里葉變換則可以讓微分和積分在頻域中變為乘法和除法
傅里葉變換相關函數
假設我們的輸入信號的函數是
$$S=0.2+0.7*\cos (2\pi*50t+\frac{20}{180}\pi)+0.2*\cos (2\pi*100t+\frac{70}{180}\pi)$$
可以發現直流分量是0.2，以及兩個余弦函數的疊加，余弦函數的幅值分別為0.7和0.2，頻率分別為50和100，初相位分別為20度和70度。
freqs = np.fft.fftfreq(采樣數量, 采樣周期) 通過采樣數與采樣周期得到時域序列經過傅里葉變換后的頻率序列
np.fft.fft(原序列) 原函數值的序列經過快速傅里葉變換得到一個復數數組，復數的模代表的是振幅，復數的輻角代表初相位
np.fft.ifft(復數序列) 復數數組 經過逆向傅里葉變換得到合成的函數值數組
案例：針對合成波做快速傅里葉變換，得到分解波數組的頻率、振幅、初相位數組，並繪制頻域圖像。
import matplotlib.pyplot as plt import numpy as np import numpy.fft as fft plt.rcParams['font.sans-serif']=['SimHei'] #用來正常顯示中文標簽 plt.rcParams['axes.unicode_minus']=False #用來正常顯示符號 Fs = 1000; # 采樣頻率 T = 1/Fs; # 采樣周期 L = 1000; # 信號長度 t = [i*T for i in range(L)] t = np.array(t) S = 0.2+0.7*np.cos(2*np.pi*50*t+20/180*np.pi) + 0.2*np.cos(2*np.pi*100*t+70/180*np.pi) ; complex_array = fft.fft(S) print(complex_array.shape) # (1000,) print(complex_array.dtype) # complex128 print(complex_array[1]) # (-2.360174309695419e-14+2.3825789764340993e-13j) ################################# plt.subplot(311) plt.grid(linestyle=':') plt.plot(1000*t[1:51], S[1:51], label='S') # y是1000個相加后的正弦序列 plt.xlabel("t（毫秒）") plt.ylabel("S(t)幅值") plt.title("疊加信號圖") plt.legend() ################################### plt.subplot(312) S_ifft = fft.ifft(complex_array) # S_new是ifft變換后的序列 plt.plot(1000*t[1:51], S_ifft[1:51], label='S_ifft', color='orangered') plt.xlabel("t（毫秒）") plt.ylabel("S_ifft(t)幅值") plt.title("ifft變換圖") plt.grid(linestyle=':') plt.legend() ################################### # 得到分解波的頻率序列 freqs = fft.fftfreq(t.size, t[1] - t[0]) # 復數的模為信號的振幅（能量大小） pows = np.abs(complex_array) plt.subplot(313) plt.title('FFT變換,頻譜圖') plt.xlabel('Frequency 頻率') plt.ylabel('Power 功率') plt.tick_params(labelsize=10) plt.grid(linestyle=':') plt.plot(freqs[freqs > 0], pows[freqs > 0], c='orangered', label='Frequency') plt.legend() plt.tight_layout() plt.show()
clear clc close all Fs = 1000; % Sampling frequency T = 1/Fs; % Sampling period L = 1000; % Length of signal t = (0:L-1)*T; % Time vector S = 0.2-0.7*cos(2*pi*50*t+20/180*pi) + 0.2*cos(2*pi*100*t+70/180*pi) ; plot(1000*t(1:50),S(1:50)) title('疊加信號圖') xlabel('t (milliseconds)') ylabel('S(t)') figure Y = fft(S); P2 = abs(Y/L); P1 = P2(1:L/2+1); P1(2:end-1) = 2*P1(2:end-1); f = Fs*(0:(L/2))/L; plot(f,P1,'linewidth',2) title('FFT變換') xlabel('頻率(Hz)') ylabel('幅值') figure pred_X=ifft(Y); plot(1000*t(1:50),pred_X(1:50),'r-')
直流分量：就是傅里葉變換的第一個值，一般來說是非復數
幅值：$2*abs(\frac{Y各項}{采樣長度})$
初相位：$atan2(\frac{Y的虛部}{Y的實部})$轉角度制表示，rad2deg(atan2(imag(Y),real(Y)))
基於傅里葉變換的頻域濾波
從某條曲線中除去一些特定的頻率成份，這在工程上稱為“濾波”。
含噪信號是高能信號與低能噪聲疊加的信號，可以通過傅里葉變換的頻域濾波實現降噪。
通過FFT使含噪信號轉換為含噪頻譜，去除低能噪聲，留下高能頻譜后再通過IFFT留下高能信號。
案例：基於傅里葉變換的頻域濾波為音頻文件去除噪聲(noiseed.wav數據集地址)。
1、讀取音頻文件，獲取音頻文件基本信息：采樣個數，采樣周期，與每個采樣的聲音信號值。繪制音頻時域的：時間/位移圖像
import numpy as np import numpy.fft as nf import scipy.io.wavfile as wf import matplotlib.pyplot as plt # 讀取音頻文件 sample_rate, noised_sigs = wf.read('./da_data/noised.wav') print(sample_rate) # sample_rate：采樣率44100 print(noised_sigs.shape) # noised_sigs:存儲音頻中每個采樣點的采樣位移(220500,) times = np.arange(noised_sigs.size) / sample_rate plt.figure('Filter') plt.subplot(221) plt.title('Time Domain', fontsize=16) plt.ylabel('Signal', fontsize=12) plt.tick_params(labelsize=10) plt.grid(linestyle=':') plt.plot(times[:178], noised_sigs[:178], c='orangered', label='Noised') plt.legend()
2、基於傅里葉變換，獲取音頻頻域信息，繪制音頻頻域的：頻率/能量圖像
# 傅里葉變換后，繪制頻域圖像 freqs = nf.fftfreq(times.size, times[1] - times[0]) complex_array = nf.fft(noised_sigs) pows = np.abs(complex_array) plt.subplot(222) plt.title('Frequency Domain', fontsize=16) plt.ylabel('Power', fontsize=12) plt.tick_params(labelsize=10) plt.grid(linestyle=':') # 指數增長坐標畫圖 plt.semilogy(freqs[freqs > 0], pows[freqs > 0], c='limegreen', label='Noised') plt.legend()
3、將低頻噪聲去除后繪制音頻頻域的：頻率/能量圖像
# 尋找能量最大的頻率值 fund_freq = freqs[pows.argmax()] # where函數尋找那些需要抹掉的復數的索引 noised_indices = np.where(freqs != fund_freq) # 復制一個復數數組的副本，避免污染原始數據 filter_complex_array = complex_array.copy() filter_complex_array[noised_indices] = 0 filter_pows = np.abs(filter_complex_array) plt.subplot(224) plt.xlabel('Frequency', fontsize=12) plt.ylabel('Power', fontsize=12) plt.tick_params(labelsize=10) plt.grid(linestyle=':') plt.plot(freqs[freqs >= 0], filter_pows[freqs >= 0], c='dodgerblue', label='Filter') plt.legend()
4、基於逆向傅里葉變換，生成新的音頻信號，繪制音頻時域的：時間/位移圖像
filter_sigs = nf.ifft(filter_complex_array).real plt.subplot(223) plt.xlabel('Time', fontsize=12) plt.ylabel('Signal', fontsize=12) plt.tick_params(labelsize=10) plt.grid(linestyle=':') plt.plot(times[:178], filter_sigs[:178], c='hotpink', label='Filter') plt.legend()
5、重新生成音頻文件
# 生成音頻文件 wf.write('./da_data/filter.wav', sample_rate, filter_sigs) plt.show()
離散傅里葉變換（DFT）
離散傅里葉變換（DFT）對有限長時域離散信號的頻譜進行等間隔采樣，頻域函數被離散化了，便於信號的計算機處理。DFT的運算量太大，FFT是離散傅里葉變換的快速算法。
說白了FFT和DFT它倆就是一個東東，只不過復雜度不同，
有時候我們能夠看到N點傅里葉變換，我個人理解是這個N點是信號前面N個連續的數值，即N點FFT意思就是截取前面N個信號進行FFT，這樣就要求我們的前N個采樣點必須包含當前信號的一個周期，不然提取的余弦波參數與正確的疊加波的參數相差很大。
如果在N點FFT的時候，如果這N個采樣點不包含一個周期呢？或者說我們的信號壓根不是一個周期函數咋辦？或者有一段是噪音數據呢？如果用FFT計算，就會對整體結果影響很大，然后就有人想通過局部來逼近整體，跟微積分的思想很像，將信號分成一小段一小段，然后對每一小段做FFT，就跟分段函數似的，無數個分段函數能逼近任意的曲線((⊙o⊙)…應該沒錯吧)，這樣每一段都不會互相影響到了。
二、短時傅里葉變換stft
在短時傅里葉變換過程中，窗的長度決定頻譜圖的時間分辨率和頻率分辨率，窗長越長，截取的信號越長，信號越長，傅里葉變換后頻率分辨率越高，時間分辨率越差；相反，窗長越短，截取的信號就越短，頻率分辨率越差，時間分辨率越好，也就是說短時傅里葉變換中，時間分辨率和頻率分辨率之間不能兼得，應該根據具體需求進行取舍。
計算短時傅里葉變換，需要指定的有：
- 每個窗口的長度：nsc
- 每相鄰兩個窗口的重疊率：nov
- 每個窗口的FFT采樣點數：nff
可以計算的有：
- 海明窗：w=hamming(nsc, 'periodic')
- 信號被分成了多少片
- 短時傅里葉變換：
python庫librosa實現
librosa.stft(y, n_fft=2048, hop_length=None, win_length=None,
window='hann', center=True, pad_mode='reflect')
短時傅立葉變換（STFT），返回一個復數矩陣使得D(f,t)
- 復數的實部：np.abs(D(f,t))頻率的振幅
- 復數的虛部：np.angle(D(f,t))頻率的相位
參數：
- y：音頻時間序列
- n_fft：FFT窗口大小，n_fft=hop_length+overlapping
- hop_length：幀移，如果未指定，則默認win_length / 4。
- win_length：每一幀音頻都由window（）加窗。窗長win_length，然后用零填充以匹配N_FFT。默認
win_length=n_fft。 - window：字符串，元組，數字，函數 shape =（n_fft, )
- 窗口（字符串，元組或數字）；
- 窗函數，例如
scipy.signal.hanning - 長度為n_fft的向量或數組
- center：bool
- 如果為True，則填充信號y，以使幀 D [:, t]以y [t * hop_length]為中心。
- 如果為False，則D [:, t]從y [t * hop_length]開始
- dtype：D的復數值類型。默認值為64-bit complex復數
- pad_mode：如果center = True，則在信號的邊緣使用填充模式。默認情況下，STFT使用reflection padding。
返回：
- STFT矩陣D，shape =（1 + $\frac{n_{fft} }{2}$，t）
librosa.magphase(D, power=1)
將復值頻譜D分離成其幅度（S）和相位（P）
參數：
- D：復值譜圖，np.ndarray [shape =（d，t），dtype = complex]
- power：幅度譜圖的指數，例如，1代表能量，2代表功率，等等。
返回：
- D_mag ：D的幅值，np.ndarray [shape =（d，t），dtype = real]
- D_phase ：D的相位，np.ndarray [shape =（d，t），dtype = complex]，exp（1.j * phi）其中phi是D的相位
y, _ = librosa.load("p225_002.wav", sr=16000) D = librosa.stft(y, n_fft=2048, hop_length=None, win_length=None, window='hann', center=True, pad_mode='reflect') print(D.shape) # (1025, 127) # 將復值頻譜D分離成其幅度(S)和相位(P)的部件 magnitude, phase = librosa.magphase(D, power=1) # magnitude # 賦值 [shape =（d，t），dtype = real] (1025, 127) # phase.shape # 相位 [shape =（d，t），dtype = complex] (1025, 127) angle = np.angle(phase) # 獲取相角（以弧度為單位）
tensorflow實現
一句話實現分幀、加窗、STFT變換
# [batch_size, signal_length]. batch_size and signal_length 可能會不知道 signals = tf.placeholder(tf.float32, [None, None]) # `stfts` 短時傅里葉變換：就是對信號中每一幀信號進行傅里葉變換 # shape is [batch_size, ?, fft_unique_bins] # 其中 fft_unique_bins = fft_length // 2 + 1 = 513. stfts = tf.contrib.signal.stft(signals, frame_length=1024, frame_step=512, fft_length=1024)
wlen：窗長
frame_length是信號中幀的長度
frame_step是幀移
fft_length：做fft變換的長度，或一種說話：fft變換所取得N點數，在有些地方也表示為NFFT。
注意：FFT的長度必須大於或者等於win的長度或者幀長。以獲得更高的頻域分辨率
FFT后的分辨率（頻率的間隔）為fs/NFFT。當NFFT>wlen時就是在數據補零后做FFT，做的FFT得到的頻譜等於以wlen長數據FFT的頻譜中內插。
numpy庫實現
上面的一行代碼相當於下面一大段代碼
def wav_to_frame(wave_data, win_len, win_shift): """ 進行分幀操作 :param wave_data: 原始的數據 :param win_len: 滑動窗長 :param win_shift: 滑動間隔 :return: 分幀之后的結果，輸出一個幀矩陣 """ num_frames = (len(wave_data) - win_len) // win_shift + 1 results = [] for i in range(num_frames): results.append(wave_data[i*win_shift:i*win_shift + win_len]) return np.array(results) def spectrum_power(frames, NFFT): """ 計算每一幀傅立葉變換以后的功率譜 參數說明： frames:audio2frame函數計算出來的幀矩陣 NFFT:FFT的大小 """ # 功率譜等於每一點的幅度平方/NFFT return 1.0/NFFT * np.square(spectrum_magnitude(frames, NFFT)) def spectrum_magnitude(frames, NFFT): """計算每一幀經過FFT變幻以后的頻譜的幅度，若frames的大小為N*L,則返回矩陣的大小為N*NFFT 參數： frames:即audio2frame函數中的返回值矩陣，幀矩陣 NFFT:FFT變換的數組大小,如果幀長度小於NFFT，則幀的其余部分用0填充鋪滿 """ complex_spectrum = np.fft.rfft(frames, NFFT) # 對frames進行FFT變換 # 返回頻譜的幅度值 return np.absolute(complex_spectrum)
三、frequency bin
在讀paper的時候總是會遇到 frequency bin （頻率窗口）這個詞，
frequency bin 是指：raw data 經過FFT后得到的頻譜圖（頻域率）中，頻率軸的頻率間隔或分辨率，通常取決采樣率和采樣點。
$$frequency \quad bin=\frac{采樣率}{采樣點數}=\frac{f_{sample}}{N_{recode}}$$
$N_{recode}$ 是信號在時域的采樣點數，頻譜中的頻率點或線的數量為$\frac{N_{recode}}{2}$ （奈奎斯特采樣定理）
頻譜的第一個頻點始終為直流（頻率=0），最后一個頻點為$\frac{f_{sample}}{2}-\frac{f_{sample}}{N_{recode}}$ 。頻點采用相等的間隔，這間隔通常用frequency bin（頻率窗口）或FFT bin表示。
例子1：我們可以作用82MHz的采樣頻率，取得8200個數據記錄，frequency bin$=\frac{82000000}{8200}=10000=10kHz$。
例子2：frequency bin是頻率域中采樣點之間的間隔。例如，如果采樣率為100赫茲，FFT為100個點，frequency bin=1，則在[0 100）赫茲之間有100個點。因此，您將整個100赫茲范圍划分為100個間隔，如0-1赫茲、1-2赫茲等。每一個如此小的間隔，比如0-1Hz，都是一個frequency bin（頻率箱）。