快速傅里叶变换python_快速傅里叶变换及python代码实现
一、前言
我想認(rèn)真寫好快速傅里葉變換（Fast Fourier Transform，FFT），所以這篇文章會由淺到細(xì)，由窄到寬的講解，但是傅里葉變換對于尋常人并不是很容易理解的，所以對于基礎(chǔ)不牢的人我會通過前言普及一下相關(guān)知識。
我們復(fù)習(xí)一下三角函數(shù)的標(biāo)準(zhǔn)式：
$$y=A\cos (\omega x+\theta )+k$$
$A$代表振幅，函數(shù)周期是$\frac{2\pi}{w}$，頻率是周期的倒數(shù)$\frac{w}{2\pi}$，$\theta $是函數(shù)初相位，$k$在信號處理中稱為直流分量。這個信號在頻域就是一條豎線。
我們再來假設(shè)有一個比較復(fù)雜的時域函數(shù)$y=f(t)$，根據(jù)傅里葉的理論，任何一個周期函數(shù)可以被分解為一系列振幅A，頻率$\omega$或初相位$\theta $正弦函數(shù)的疊加
$$y = A_1sin(\omega_1t+\theta_1) + A_2sin(\omega_2t+\theta_2) + A_3sin(\omega_3t+\theta_3)$$
該信號在頻域有三條豎線組成，而豎線圖我們把它稱為頻譜圖，大家可以通過下面的動畫了解
如圖可知，通過時域到頻域的變換，我們得到了一個從側(cè)面看的頻譜，但是這個頻譜并沒有包含時域中全部的信息。因為頻譜只代表每個正弦波對應(yīng)頻率的振幅是多少，而沒有提到相位。基礎(chǔ)的正弦波$Asin(wt+\theta )$中，振幅，頻率，相位缺一不可，不同相位決定了波的位置，所以對于頻域分析，僅僅有頻譜（振幅譜）是不夠的，我們還需要一個相位譜。
我依稀記得高中學(xué)正弦函數(shù)的是時候，$\theta $的多少決定了正弦波向右移動多少。當(dāng)然那個時候橫坐標(biāo)是相位角度，而時域信號的橫坐標(biāo)是時間，因此我們只需要將時間轉(zhuǎn)換為相位角度就得到了初相位。相位差則是時間差在一個周期中所占的比例
$$\theta=2\pi \frac{t}{T}$$
所以傅里葉變換可以把一個比較復(fù)雜的函數(shù)轉(zhuǎn)換為多個簡單函數(shù)的疊加，將時域（即時間域）上的信號轉(zhuǎn)變?yōu)轭l域（即頻率域）上的信號，看問題的角度也從時間域轉(zhuǎn)到了頻率域，因此在時域中某些不好處理的地方，在頻域就可以較為簡單的處理，這就可以大量減少處理信號計算量。信號經(jīng)過傅里葉變換后，可以得到頻域的幅度譜以及相位譜，信號的幅度譜和相位譜是信號傅里葉變換后頻譜的兩個屬性。
傅里葉用途
時域復(fù)雜的函數(shù)，在頻域就是幾條豎線
求解微分方程，傅里葉變換則可以讓微分和積分在頻域中變?yōu)槌朔ê统?/p>
傅里葉變換相關(guān)函數(shù)
假設(shè)我們的輸入信號的函數(shù)是
$$S=0.2+0.7*\cos (2\pi*50t+\frac{20}{180}\pi)+0.2*\cos (2\pi*100t+\frac{70}{180}\pi)$$
可以發(fā)現(xiàn)直流分量是0.2，以及兩個余弦函數(shù)的疊加，余弦函數(shù)的幅值分別為0.7和0.2，頻率分別為50和100，初相位分別為20度和70度。
freqs = np.fft.fftfreq(采樣數(shù)量, 采樣周期) 通過采樣數(shù)與采樣周期得到時域序列經(jīng)過傅里葉變換后的頻率序列
np.fft.fft(原序列) 原函數(shù)值的序列經(jīng)過快速傅里葉變換得到一個復(fù)數(shù)數(shù)組，復(fù)數(shù)的模代表的是振幅，復(fù)數(shù)的輻角代表初相位
np.fft.ifft(復(fù)數(shù)序列) 復(fù)數(shù)數(shù)組 經(jīng)過逆向傅里葉變換得到合成的函數(shù)值數(shù)組
案例：針對合成波做快速傅里葉變換，得到分解波數(shù)組的頻率、振幅、初相位數(shù)組，并繪制頻域圖像。
importmatplotlib.pyplot as pltimportnumpy as npimportnumpy.fft as fft
plt.rcParams['font.sans-serif']=['SimHei'] #用來正常顯示中文標(biāo)簽
plt.rcParams['axes.unicode_minus']=False #用來正常顯示符號
Fs= 1000; #采樣頻率
T = 1/Fs; #采樣周期
L = 1000; #信號長度
t = [i*T for i inrange(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])#復(fù)數(shù)的模為信號的振幅（能量大小）
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代碼實現(xiàn)
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實現(xiàn)
直流分量：就是傅里葉變換的第一個值，一般來說是非復(fù)數(shù)
幅值：$2*abs(\frac{Y各項}{采樣長度})$
初相位：$atan2(\frac{Y的虛部}{Y的實部})$轉(zhuǎn)角度制表示，rad2deg(atan2(imag(Y),real(Y)))
基于傅里葉變換的頻域濾波
從某條曲線中除去一些特定的頻率成份，這在工程上稱為“濾波”。
含噪信號是高能信號與低能噪聲疊加的信號，可以通過傅里葉變換的頻域濾波實現(xiàn)降噪。
通過FFT使含噪信號轉(zhuǎn)換為含噪頻譜，去除低能噪聲，留下高能頻譜后再通過IFFT留下高能信號。
案例：基于傅里葉變換的頻域濾波為音頻文件去除噪聲(noiseed.wav數(shù)據(jù)集地址)。
1、讀取音頻文件，獲取音頻文件基本信息：采樣個數(shù)，采樣周期，與每個采樣的聲音信號值。繪制音頻時域的：時間/位移圖像
importnumpy as npimportnumpy.fft as nfimportscipy.io.wavfile as wfimportmatplotlib.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=':')#指數(shù)增長坐標(biāo)畫圖
plt.semilogy(freqs[freqs > 0], pows[freqs > 0], c='limegreen', label='Noised')
plt.legend()
3、將低頻噪聲去除后繪制音頻頻域的：頻率/能量圖像
#尋找能量最大的頻率值
fund_freq =freqs[pows.argmax()]#where函數(shù)尋找那些需要抹掉的復(fù)數(shù)的索引
noised_indices = np.where(freqs !=fund_freq)#復(fù)制一個復(fù)數(shù)數(shù)組的副本，避免污染原始數(shù)據(jù)
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）對有限長時域離散信號的頻譜進(jìn)行等間隔采樣，頻域函數(shù)被離散化了，便于信號的計算機(jī)處理。DFT的運(yùn)算量太大，FFT是離散傅里葉變換的快速算法。
說白了FFT和DFT它倆就是一個東東，只不過復(fù)雜度不同，
有時候我們能夠看到N點傅里葉變換，我個人理解是這個N點是信號前面N個連續(xù)的數(shù)值，即N點FFT意思就是截取前面N個信號進(jìn)行FFT，這樣就要求我們的前N個采樣點必須包含當(dāng)前信號的一個周期，不然提取的余弦波參數(shù)與正確的疊加波的參數(shù)相差很大。
如果在N點FFT的時候，如果這N個采樣點不包含一個周期呢？或者說我們的信號壓根不是一個周期函數(shù)咋辦？或者有一段是噪音數(shù)據(jù)呢？如果用FFT計算，就會對整體結(jié)果影響很大，然后就有人想通過局部來逼近整體，跟微積分的思想很像，將信號分成一小段一小段，然后對每一小段做FFT，就跟分段函數(shù)似的，無數(shù)個分段函數(shù)能逼近任意的曲線((⊙o⊙)…應(yīng)該沒錯吧)，這樣每一段都不會互相影響到了。
二、短時傅里葉變換stft
在短時傅里葉變換過程中，窗的長度決定頻譜圖的時間分辨率和頻率分辨率，窗長越長，截取的信號越長，信號越長，傅里葉變換后頻率分辨率越高，時間分辨率越差；相反，窗長越短，截取的信號就越短，頻率分辨率越差，時間分辨率越好，也就是說短時傅里葉變換中，時間分辨率和頻率分辨率之間不能兼得，應(yīng)該根據(jù)具體需求進(jìn)行取舍。
計算短時傅里葉變換，需要指定的有：
每個窗口的長度：nsc
每相鄰兩個窗口的重疊率：nov
每個窗口的FFT采樣點數(shù)：nff
可以計算的有：
海明窗：w=hamming(nsc, 'periodic')
信號被分成了多少片
短時傅里葉變換：
python庫librosa實現(xiàn)
librosa.stft(y, n_fft=2048, hop_length=None, win_length=None,
window='hann', center=True, pad_mode='reflect')
短時傅立葉變換（STFT），返回一個復(fù)數(shù)矩陣使得D(f,t)
復(fù)數(shù)的實部：np.abs(D(f,t))頻率的振幅
復(fù)數(shù)的虛部：np.angle(D(f,t))頻率的相位
參數(shù)：
y：音頻時間序列
n_fft：FFT窗口大小，n_fft=hop_length+overlapping
hop_length：幀移，如果未指定，則默認(rèn)win_length / 4。
win_length：每一幀音頻都由window（）加窗。窗長win_length，然后用零填充以匹配N_FFT。默認(rèn)win_length=n_fft。
window：字符串，元組，數(shù)字，函數(shù) shape =（n_fft, )
窗口（字符串，元組或數(shù)字）；
窗函數(shù)，例如scipy.signal.hanning
長度為n_fft的向量或數(shù)組
center：bool
如果為True，則填充信號y，以使幀 D [:, t]以y [t * hop_length]為中心。
如果為False，則D [:, t]從y [t * hop_length]開始
dtype：D的復(fù)數(shù)值類型。默認(rèn)值為64-bit complex復(fù)數(shù)
pad_mode：如果center = True，則在信號的邊緣使用填充模式。默認(rèn)情況下，STFT使用reflection padding。
返回：
STFT矩陣D，shape =（1 + $\frac{n_{fft} }{2}$，t）
librosa.magphase(D, power=1)
將復(fù)值頻譜D分離成其幅度（S）和相位（P）
參數(shù)：
D：復(fù)值譜圖，np.ndarray [shape =（d，t），dtype = complex]
power：幅度譜圖的指數(shù)，例如，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)
#將復(fù)值頻譜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實現(xiàn)
一句話實現(xiàn)分幀、加窗、STFT變換
#[batch_size, signal_length]. batch_size and signal_length 可能會不知道
signals =tf.placeholder(tf.float32, [None, None])#`stfts` 短時傅里葉變換：就是對信號中每一幀信號進(jìn)行傅里葉變換#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點數(shù)，在有些地方也表示為NFFT。
注意：FFT的長度必須大于或者等于win的長度或者幀長。以獲得更高的頻域分辨率
FFT后的分辨率（頻率的間隔）為fs/NFFT。當(dāng)NFFT>wlen時就是在數(shù)據(jù)補(bǔ)零后做FFT，做的FFT得到的頻譜等于以wlen長數(shù)據(jù)FFT的頻譜中內(nèi)插。
numpy庫實現(xiàn)
上面的一行代碼相當(dāng)于下面一大段代碼
defwav_to_frame(wave_data, win_len, win_shift):"""進(jìn)行分幀操作
:param wave_data: 原始的數(shù)據(jù)
:param win_len: 滑動窗長
:param win_shift: 滑動間隔
:return: 分幀之后的結(jié)果，輸出一個幀矩陣"""num_frames= (len(wave_data) - win_len) // win_shift + 1results=[]for i inrange(num_frames):
results.append(wave_data[i*win_shift:i*win_shift +win_len])returnnp.array(results)defspectrum_power(frames, NFFT):"""計算每一幀傅立葉變換以后的功率譜
參數(shù)說明：
frames:audio2frame函數(shù)計算出來的幀矩陣
NFFT:FFT的大小"""
#功率譜等于每一點的幅度平方/NFFT
return 1.0/NFFT *np.square(spectrum_magnitude(frames, NFFT))defspectrum_magnitude(frames, NFFT):"""計算每一幀經(jīng)過FFT變幻以后的頻譜的幅度，若frames的大小為N*L,則返回矩陣的大小為N*NFFT
參數(shù)：
frames:即audio2frame函數(shù)中的返回值矩陣，幀矩陣
NFFT:FFT變換的數(shù)組大小,如果幀長度小于NFFT，則幀的其余部分用0填充鋪滿"""complex_spectrum= np.fft.rfft(frames, NFFT) #對frames進(jìn)行FFT變換
#返回頻譜的幅度值
return np.absolute(complex_spectrum)
View Code
三、frequency bin
在讀paper的時候總是會遇到 frequency bin （頻率窗口）這個詞，
frequency bin 是指：raw data 經(jīng)過FFT后得到的頻譜圖（頻域率）中，頻率軸的頻率間隔或分辨率，通常取決采樣率和采樣點。
$$frequency \quad bin=\frac{采樣率}{采樣點數(shù)}=\frac{f_{sample}}{N_{recode}}$$
$N_{recode}$ 是信號在時域的采樣點數(shù)，頻譜中的頻率點或線的數(shù)量為$\frac{N_{recode}}{2}$ （奈奎斯特采樣定理）
頻譜的第一個頻點始終為直流（頻率=0），最后一個頻點為$\frac{f_{sample}}{2}-\frac{f_{sample}}{N_{recode}}$ 。頻點采用相等的間隔，這間隔通常用frequency bin（頻率窗口）或FFT bin表示。
例子1：我們可以作用82MHz的采樣頻率，取得8200個數(shù)據(jù)記錄，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（頻率箱）。
參考
總結(jié)
以上是生活随笔為你收集整理的快速傅里叶变换python_快速傅里叶变换及python代码实现的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 为什么新疆人吃的羊肉膻味那么重？
- 下一篇： 7月11号，大连小雨