快速傅里葉變換模塊(fft)


什么是傅里葉變換?

法國科學家傅里葉提出,任何一條周期曲線,無論多么跳躍或不規則,都能表示成一組光滑正弦曲線疊加之和。

傅里葉變換的目的是可將時域(即時間域)上的信號轉變為頻域(即頻率域)上的信號,隨着域的不同,對同一個事物的了解角度也就隨之改變,因此在時域中某些不好處理的地方,在頻域就可以較為簡單的處理。這就可以大量減少處理信號存儲量。

例如:彈鋼琴

假設有一時間域函數:y = f(x),根據傅里葉的理論它可以被分解為一系列正弦函數的疊加,他們的振幅A,頻率ω或初相位φ不同:

所以傅里葉變換可以把一個比較復雜的函數轉換為多個簡單函數的疊加,看問題的角度也從時間域轉到了頻率域,有些的問題處理起來就會比較簡單。

傅里葉變換相關函數

導入快速傅里葉變換所需模塊

import numpy.fft as nf

通過采樣數與采樣周期求得傅里葉變換分解所得曲線的頻率序列

freqs = np.fft.fftfreq(采樣數量, 采樣周期)

通過原函數值的序列j經過快速傅里葉變換得到一個復數數組,復數的模代表的是振幅,復數的輻角代表初相位

np.fft.fft(原函數值序列) -> 目標函數值序列(復數)

通過一個復數數組(復數的模代表的是振幅,復數的輻角代表初相位)經過逆向傅里葉變換得到合成的函數值數組

np.fft.ifft(目標函數值序列(復數))->原函數值序列

案例:針對合成波做快速傅里葉變換,得到一組復數序列;再針對該復數序列做逆向傅里葉變換得到新的合成波並繪制。

ffts = nf.fft(sigs6)
sigs7 = nf.ifft(ffts).real
mp.plot(times, sigs7, label=r'$\omega$='+str(round(1 / (2 * np.pi),3)), alpha=0.5, linewidth=6)

 

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

# 得到分解波的頻率序列
freqs = nf.fftfreq(times.size, times[1] - times[0])
# 復數的模為信號的振幅(能量大小)
ffts = nf.fft(sigs6)
pows = np.abs(ffts)

mp.subplot(122)
mp.title('Frequency Domain', fontsize=16)
mp.xlabel('Frequency', fontsize=12)
mp.ylabel('Power', fontsize=12)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')
mp.plot(freqs[freqs >= 0], pows[freqs >= 0], c='orangered', label='Frequency Spectrum')
mp.legend()
mp.tight_layout()
mp.show()

 

 

"""
傅里葉變換拆方波
"""
import numpy as np
import matplotlib.pyplot as mp
import numpy.fft as nf

x = np.linspace(-2*np.pi, 2*np.pi, 1000)
y = np.zeros(x.size)
for i in range(1, 1000):
    y += 4*np.pi/(2*i-1) * np.sin((2*i-1)*x)

mp.figure('FFT',facecolor='lightgray')
mp.subplot(121)
mp.title('Time Domain',fontsize=18)
mp.grid(linestyle=":")
mp.plot(x, y, label=r'$y$')
#針對方波y,做fft
comp_ary = nf.fft(y)
y2 = nf.ifft(comp_ary)
mp.plot(x,y2,color = 'orangered',linewidth=5,alpha=0.5,label='$y$')
#繪制頻域圖像(頻率/能量)
mp.subplot(122)
freqs = nf.fftfreq(y.size,x[1]-x[0])
pows = np.abs(comp_ary)#負數的模
mp.title('Frequency Domain',fontsize=18)
mp.grid(linestyle=":")
mp.plot(freqs[freqs>0],pows[freqs>0],color= 'orangered',label='frequency')

mp.legend()
mp.show()

 

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

含噪信號是高能信號與低能噪聲疊加的信號,可以通過傅里葉變換的頻域濾波實現降噪。

通過FFT使含噪信號轉換為含噪頻譜,去除低能噪聲,留下高能頻譜后再通過IFFT留下高能信號。

案例:基於傅里葉變換的頻域濾波為音頻文件去除噪聲。

讀取音頻文件,獲取音頻文件基本信息:采樣個數,采樣周期,與每個采樣的聲音信號值。繪制音頻時域的:時間/位移圖像。

import numpy as np
import numpy.fft as nf
import scipy.io.wavfile as wf
import matplotlib.pyplot as mp

sample_rate, noised_sigs = wf.read('../data/noised.wav')
noised_sigs = noised_sigs / 2 ** 15
times = np.arange(len(noised_sigs)) / sample_rate
mp.figure('Filter', facecolor='lightgray')
mp.subplot(221)
mp.title('Time Domain', fontsize=16)
mp.ylabel('Signal', fontsize=12)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')
mp.plot(times[:178], noised_sigs[:178],c='orangered', label='Noised')
mp.legend()
mp.show()

 

基於傅里葉變換,獲取音頻頻域信息,繪制音頻頻域的:頻率/能量圖像。

freqs = nf.fftfreq(times.size, 1 / sample_rate)
noised_ffts = nf.fft(noised_sigs)
noised_pows = np.abs(noised_ffts)
mp.subplot(222)
mp.title('Frequency Domain', fontsize=16)
mp.ylabel('Power', fontsize=12)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')
mp.semilogy(freqs[freqs >= 0],noised_pows[freqs >= 0], c='limegreen',label='Noised')
mp.legend()

 

將低能噪聲去除后繪制音頻頻域的:頻率/能量圖像。

fund_freq = freqs[noised_pows.argmax()]
noised_indices = np.where(freqs != fund_freq)
filter_ffts = noised_ffts.copy()
filter_ffts[noised_indices] = 0
filter_pows = np.abs(filter_ffts)

mp.subplot(224)
mp.xlabel('Frequency', fontsize=12)
mp.ylabel('Power', fontsize=12)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')
mp.plot(freqs[freqs >= 0], filter_pows[freqs >= 0],c='dodgerblue', label='Filter')
mp.legend() 

 

基於逆向傅里葉變換,生成新的音頻信號,繪制音頻時域的:時間/位移圖像。

filter_sigs = nf.ifft(filter_ffts).real
mp.subplot(223)
mp.xlabel('Time', fontsize=12)
mp.ylabel('Signal', fontsize=12)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')
mp.plot(times[:178], filter_sigs[:178],c='hotpink', label='Filter')
mp.legend()

 

重新生成音頻文件。

wf.write('../../data/filter.wav',sample_rate,(filter_sigs * 2 ** 15).astype(np.int16))

 

 

"""
頻域濾波
"""
import numpy as np
import numpy.fft as nf
import matplotlib.pyplot as mp
import scipy.io.wavfile as wf

# 采樣率,采樣位移
sample_rate, sigs = wf.read('noised.wav')
print(sample_rate)  # 44100
print(sigs.shape)  # (220500,)
sigs = sigs/(2**15)
times = np.arange(sigs.size) / sample_rate
print(times)
"""
[0.00000000e+00 2.26757370e-05 4.53514739e-05 ... 4.99993197e+00
 4.99995465e+00 4.99997732e+00]
"""
# 繪圖
mp.figure('Filter', facecolor='lightgray')
mp.subplot(221)
mp.title('Time Domain', fontsize=16)
mp.ylabel('Signal', fontsize=14)
mp.grid(linestyle=':')
mp.plot(times[:178], sigs[:178], color='dodgerblue', label='Noised Signal')
mp.legend()

# fft變換后繪制頻域圖(頻率/能量)
comp_ary = nf.fft(sigs)
pows = np.abs(comp_ary)
freqs = nf.fftfreq(sigs.size, times[1] - times[0])
mp.subplot(222)
mp.title('Frequency Domain', fontsize=16)
mp.ylabel('pow', fontsize=14)
mp.grid(linestyle=":")
mp.semilogy(freqs[freqs > 0], pows[freqs > 0], color='orangered', label='Freq')

mp.legend()

# 在頻譜中去除噪聲
maxpow_freq = freqs[pows.argmax()]
noised_inds = np.where(freqs != maxpow_freq)
# 噪聲的復數全抹為0
comp_ary[noised_inds] = 0
# 繪制降噪之后的頻譜圖像
pows = np.abs(comp_ary)
mp.subplot(224)
mp.ylabel('pow', fontsize=14)
mp.grid(linestyle=':')
mp.plot(freqs[freqs > 0], pows[freqs > 0],
        color='orangered', label='Freq')
mp.legend()


# 逆向傅里葉變換,輸出時域函數圖像
filter_sigs = nf.ifft(comp_ary)
mp.subplot(223)
mp.ylabel('Signal', fontsize=14)
mp.grid(linestyle=":")
mp.plot(times[:178], filter_sigs[:178], color='dodgerblue', label='Filter Signal')
mp.legend()
# 保存文件
wf.write('filter.wav', sample_rate, (filter_sigs * 2 ** 15).astype(np.int16))

mp.tight_layout()
mp.show()

 

 

 

 

 

 

 

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM