文章來源：https://www.cnblogs.com/LXP-Never/p/11558302.html

---

快速傅里葉變換及python代碼實現

目錄

一、前言
??傅里葉變換相關函數
??基于傅里葉變換的頻域濾波
??離散傅里葉變換（DFT）
二、短時傅里葉變換stft
三、frequency bin
參考

一、前言

　　我想認真寫好快速傅里葉變換（Fast Fourier Transform，FFT），所以這篇文章會由淺到細，由窄到寬的講解，但是傅里葉變換對于尋常人并不是很容易理解的，所以對于基礎不牢的人我會通過前言普及一下相關知識。

　　我們復習一下三角函數的標準式：

y=Acos(ωx+θ)+ky=Acos?(ωx+θ)+k

　　AA代表振幅，函數周期是2πw2πw，頻率是周期的倒數w2πw2π，θθ是函數初相位，kk在信號處理中稱為直流分量。這個信號在頻域就是一條豎線。

　　我們再來假設有一個比較復雜的時域函數y=f(t)y=f(t)，根據傅里葉的理論，任何一個周期函數可以被分解為一系列振幅A，頻率ωω或初相位θθ正弦函數的疊加

?

y=A1sin(ω1t+θ1)+A2sin(ω2t+θ2)+A3sin(ω3t+θ3)y=A1sin(ω1t+θ1)+A2sin(ω2t+θ2)+A3sin(ω3t+θ3)

　　該信號在頻域有三條豎線組成，而豎線圖我們把它稱為頻譜圖，大家可以通過下面的動畫了解

　　如圖可知，通過時域到頻域的變換，我們得到了一個從側面看的頻譜，但是這個頻譜并沒有包含時域中全部的信息。因為頻譜只代表每個正弦波對應頻率的振幅是多少，而沒有提到相位。基礎的正弦波Asin(wt+θ)Asin(wt+θ)中，振幅，頻率，相位缺一不可，不同相位決定了波的位置，所以對于頻域分析，僅僅有頻譜（振幅譜）是不夠的，我們還需要一個相位譜。

　　我依稀記得高中學正弦函數的是時候，θθ的多少決定了正弦波向右移動多少。當然那個時候橫坐標是相位角度，而時域信號的橫坐標是時間，因此我們只需要將時間轉換為相位角度就得到了初相位。相位差則是時間差在一個周期中所占的比例

?

θ=2πtTθ=2πtT

　　所以傅里葉變換可以把一個比較復雜的函數轉換為多個簡單函數的疊加，將時域（即時間域）上的信號轉變為頻域（即頻率域）上的信號，看問題的角度也從時間域轉到了頻率域，因此在時域中某些不好處理的地方，在頻域就可以較為簡單的處理，這就可以大量減少處理信號計算量。信號經過傅里葉變換后，可以得到頻域的幅度譜以及相位譜，信號的幅度譜和相位譜是信號傅里葉變換后頻譜的兩個屬性。

傅里葉用途

• 時域復雜的函數，在頻域就是幾條豎線
• 求解微分方程，傅里葉變換則可以讓微分和積分在頻域中變為乘法和除法

傅里葉變換相關函數

　　假設我們的輸入信號的函數是

S=0.2+0.7?cos(2π?50t+20180π)+0.2?cos(2π?100t+70180π)S=0.2+0.7?cos?(2π?50t+20180π)+0.2?cos?(2π?100t+70180π)

可以發現直流分量是0.2，以及兩個余弦函數的疊加，余弦函數的幅值分別為0.7和0.2，頻率分別為50和100，初相位分別為20度和70度。

freqs = np.fft.fftfreq(采樣數量, 采樣周期)--通過采樣數與采樣周期得到時域序列經過傅里葉變換后的頻率序列

np.fft.fft(原序列)--原函數值的序列經過快速傅里葉變換得到一個復數數組，復數的模代表的是振幅，復數的輻角代表初相位

np.fft.ifft(復數序列)--復數數組 經過逆向傅里葉變換得到合成的函數值數組

案例：針對合成波做快速傅里葉變換，得到分解波數組的頻率、振幅、初相位數組，并繪制頻域圖像。

python代碼實現

import matplotlib.pyplot as plt import numpy as np import numpy.fft as fftplt.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()python代碼實現

MATLAB實現

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-')MATLAB實現

直流分量：就是傅里葉變換的第一個值，一般來說是非復數

幅值：2?abs(Y各項采樣長度)2?abs(Y各項采樣長度)

初相位：atan2(Y的虛部Y的實部)atan2(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_rateplt.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 +?nfft2nfft2，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 + 1results = []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的大小"""# 功率譜等于每一點的幅度平方/NFFTreturn 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后得到的頻譜圖（頻域率）中，頻率軸的頻率間隔或分辨率，通常取決采樣率和采樣點。

　　

例子1：我們可以作用82MHz的采樣頻率，取得8200個數據記錄，frequency bin=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（頻率箱）。

參考

知乎Heinrich博主文章——傅里葉分析之掐死教程

【音頻處理】離散傅里葉變換

【音頻處理】短時傅里葉變換

總結

以上是生活随笔為你收集整理的快速傅里叶变换及python代码实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。