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

時間  2020-05-08
標籤 快速 傅里葉 變換 python 代碼 實現 欄目 Python 简体版
原文   原文鏈接

0 前言

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

　　咱們複習一下三角函數的標準式：python

$$y=A\cos (\omega x+\theta )+k$$git

　　$A$表明振幅，函數週期是$\frac{2\pi}{w}$，頻率是週期的倒數$\frac{w}{2\pi}$，$\theta $是函數初相位，$k$在信號處理中稱爲直流份量。這個信號在頻域就是一條豎線。github

　　咱們再來假設有一個比較複雜的時域函數$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)$$數組

　　該信號在頻域有三條豎線組成，而豎線圖咱們把它稱爲頻譜圖，你們能夠經過下面的動畫了解app

　　如圖可知，經過時域到頻域的變換，咱們獲得了一個從側面看的頻譜，可是這個頻譜並無包含時域中所有的信息。由於頻譜只表明每個對應的正弦波的振幅是多少，而沒有提到相位。基礎的正弦波$Asin(wt+\theta )$中，振幅，頻率，相位缺一不可，不一樣相位決定了波的位置，因此對於頻域分析，僅僅有頻譜（振幅譜）是不夠的，咱們還須要一個相位譜。ide

　　我依稀記得高中學正弦函數的是時候，$\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()
python代碼實現 

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(\frac{Y各項}{採樣長度})$

初相位：$atan2(\frac{Y的虛部}{Y的實部})$轉角度製表示，rad2deg(atan2(imag(Y),real(Y)))

 

基於傅里葉變換的頻域濾波

從某條曲線中除去一些特定的頻率成份，這在工程上稱爲「濾波」。

含噪信號是高能信號與低能噪聲疊加的信號，能夠經過傅里葉變換的頻域濾波實現降噪。

經過FFT使含噪信號轉換爲含噪頻譜，去除低能噪聲，留下高能頻譜後再經過IFFT留下高能信號。

案例：基於傅里葉變換的頻域濾波爲音頻文件去除噪聲(noiseed.wav數據集地址)。

　　一、讀取音頻文件，獲取音頻文件基本信息：採樣個數，採樣週期，與每一個採樣的聲音信號值。繪製音頻時域的：時間/位移圖像

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

 

　　二、基於傅里葉變換，獲取音頻頻域信息，繪製音頻頻域的：頻率/能量圖像

# 傅里葉變換後，繪製頻域圖像
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()

　　三、將低頻噪聲去除後繪製音頻頻域的：頻率/能量圖像

# 尋找能量最大的頻率值
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()

　　四、基於逆向傅里葉變換，生成新的音頻信號，繪製音頻時域的：時間/位移圖像

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

　　五、從新生成音頻文件

# 生成音頻文件
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矩陣，shape =（1 + $\frac{n_{fft} }{2}$，t）

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)
View Code

 

frequency bin

在讀paper的時候老是會遇到 frequency bin （頻率窗口）這個詞，

　　raw data 通過FFT後獲得的頻譜圖中，頻率軸的頻率間隔或分辨率，一般取決採樣速率和採樣點數據。
　　功率譜中的頻率點或線的數量爲$\frac{N_{recode}}{2}$  ，$N_{recode}$  是信號在時域的採樣點數。
　　功率譜的第一個頻點始終爲直流（頻率=0），最後一個頻點爲$\frac{f_{sample}}{2}-\frac{f_{sample}}{N_{recode}}$  。 頻點採用相等的間隔$\frac{f_{sample}}{N_{recode}}$  ，一般用頻率窗口或FFT窗口表示。窗口（Bin）能夠由數據轉換器的採樣週期計算：
$$bin=採樣率/採樣點數=\frac{f_{sample}}{N_{recode}}=\frac{1}{N_{recode}*\Delta t_{sample}}$$
例如：咱們能夠做用82MHz的採樣頻率，取得8192個數據記錄，頻率間隔爲10kHz。

參考

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

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

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

原文出處：https://www.cnblogs.com/LXP-Never/p/11558302.html

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程