傅里葉變換是在高數(shù)是一個(gè)很重要的知識(shí)點(diǎn)，今天將結(jié)合Python代碼實(shí)現(xiàn)傅立葉變換。
傅立葉變換
我們平時(shí)是如何去分解一個(gè)復(fù)雜的問(wèn)題呢？一個(gè)經(jīng)典的方法就是把這個(gè)復(fù)雜的問(wèn)題分解成為多個(gè)簡(jiǎn)單的可操作的子問(wèn)題， 傅立葉變換也是基于這個(gè)思想。
傅里葉分析是研究如何將數(shù)學(xué)函數(shù)分解為一系列更簡(jiǎn)單的三角函數(shù)的領(lǐng)域。傅里葉變換是該領(lǐng)域的一種工具，用于將函數(shù)分解為其分量頻率。
在本教程中，傅立葉變換是一種工具，可以獲取信號(hào)并查看其中每個(gè)頻率的功率。看一看該傅立葉變換中的重要術(shù)語(yǔ)：
- 信號(hào)：信號(hào)是隨時(shí)間變化的信息。例如，音頻，視頻和電壓走線都是信號(hào)的示例。
- 頻率：頻率是某物重復(fù)的速度。例如，時(shí)鐘以1赫茲（Hz）的頻率滴答，或每秒重復(fù)1次。
- 功率：功率表示每個(gè)頻率的強(qiáng)度。
下圖是一些正弦波的頻率和功率的直觀演示：
第一個(gè)是低頻正弦波，第二個(gè)是高頻正弦波，第三個(gè)是低頻低功率正弦波，因此低功率正弦波比其它兩個(gè)正弦波的峰較小。
時(shí)域與頻域
時(shí)域與頻域是查看信號(hào)的兩種不同方式，即信號(hào)的組成頻率或隨時(shí)間變化的信息。
在時(shí)域中，信號(hào)是隨時(shí)間（x軸）幅度（y軸）變化的波。您最有可能在時(shí)域中查看圖表，例如：
這是一些音頻的圖像，它是一個(gè)時(shí)域信號(hào)。橫軸表示時(shí)間，縱軸表示振幅。
在頻域中，信號(hào)表示為一系列頻率（x軸），每個(gè)頻率都具有關(guān)聯(lián)的功率（y軸）。下圖是經(jīng)過(guò)傅立葉變換后的上述音頻信號(hào)：
代碼實(shí)現(xiàn)正弦波
音頻本質(zhì)上是正弦波。
下面是產(chǎn)生正弦波的代碼：
import numpy as np
from matplotlib import pyplot as plt
SAMPLE_RATE = 44100 # 赫茲
DURATION = 5 # 秒
def generate_sine_wave(freq, sample_rate, duration):
x = np.linspace(0, duration, sample_rate * duration, endpoint=False)
frequencies = x * freq
y = np.sin((2 * np.pi) * frequencies)
return x, y
# 產(chǎn)生持續(xù)5秒的2赫茲正弦波
x, y = generate_sine_wave(2, SAMPLE_RATE, DURATION)
plt.plot(x, y)
plt.show()
x軸以秒為單位表示時(shí)間，并且由于每秒鐘的時(shí)間都有兩個(gè)峰值，因此可以看到正弦波每秒振蕩兩次。
混合音頻
下面將兩個(gè)正弦波，混合音頻信號(hào)僅包括兩個(gè)步驟：
將正弦波加在一起，然后進(jìn)行歸一化的操作。
具體實(shí)現(xiàn)的代碼如下。
_, nice_tone = generate_sine_wave(400, SAMPLE_RATE, DURATION)
_, noise_tone = generate_sine_wave(4000, SAMPLE_RATE, DURATION)
noise_tone = noise_tone * 0.3
mixed_tone = nice_tone + noise_tone
下一步是歸一化，或縮放信號(hào)以適合目標(biāo)格式。由于以后將如何存儲(chǔ)音頻，目標(biāo)格式為16位整數(shù)，范圍為-32768到32767：
normalized_tone = np.int16((mixed_tone / mixed_tone.max()) * 32767)
plt.plot(normalized_tone[:1000])
plt.show()
看到的正弦波是生成的400 Hz音調(diào)，將上面的正弦波轉(zhuǎn)化為音頻，最簡(jiǎn)單的方法是使用SciPy的wavfile.write方法將其存儲(chǔ)在WAV文件中。16位整數(shù)是WAV文件的標(biāo)準(zhǔn)數(shù)據(jù)類型，因此需要將信號(hào)標(biāo)準(zhǔn)化為16位整數(shù)：
from scipy.io.wavfile import write
# 記住，采樣率=44100赫茲是我們的播放率
write("mysinewave.wav", SAMPLE_RATE, normalized_tone)
這個(gè)音頻聽(tīng)起來(lái)音調(diào)很高。
完成此步驟后，就當(dāng)作音頻樣本了。下一步是使用傅立葉變換消除高音調(diào)！
傅立葉變換
現(xiàn)在對(duì)生成的音頻上使用FFT了。FFT是一種算法，可實(shí)現(xiàn)傅立葉變換并可以在時(shí)域中為信號(hào)計(jì)算頻譜。
from scipy.fft import fft, fftfreq
# 標(biāo)準(zhǔn)化音調(diào)中的樣本數(shù)
N = SAMPLE_RATE * DURATION
yf = fft(normalized_tone)
xf = fftfreq(N, 1 / SAMPLE_RATE)
plt.plot(xf, np.abs(yf))
plt.show()
我們可以在正頻率中看到兩個(gè)峰值，正頻率峰值位于400 Hz和4000 Hz，與之前生成的音頻的頻率相對(duì)應(yīng)。
計(jì)算傅里葉變換
yf = fft(normalized_tone)
xf = fftfreq(N, 1 / SAMPLE_RATE)
上面代碼的功能
- fft() 計(jì)算轉(zhuǎn)換本身。
- fftfreq()計(jì)算的輸出中每個(gè)倉(cāng)中心的頻率fft()。沒(méi)有這個(gè)，就無(wú)法在頻譜上繪制x軸
fft()輸出的頻譜圍繞y軸反射，因此負(fù)半部分是正半部分的鏡像，我們一般只需計(jì)算一半對(duì)稱值，即可更快地進(jìn)行傅立葉變換。scipy.fft以的形式實(shí)施此速度駭客rfft()。
from scipy.fft import rfft, rfftfreq
# 注意前面多余的“r”
yf = rfft(normalized_tone)
xf = rfftfreq(N, 1 / SAMPLE_RATE)
plt.plot(xf, np.abs(yf))
plt.show()
過(guò)濾信號(hào)
傅里葉變換的一大優(yōu)點(diǎn)是它是可逆的，我們可以利用此優(yōu)勢(shì)來(lái)過(guò)濾音頻并擺脫高音調(diào)頻率。
# 最大頻率為采樣率的一半
points_per_freq = len(xf) / (SAMPLE_RATE / 2)
# 我們的目標(biāo)頻率是4000赫茲 將44100變成4000
target_idx = int(points_per_freq * 4000)
然后，您可以將其設(shè)置yf為0目標(biāo)頻率附近的index來(lái)擺脫它：
yf[target_idx - 1 : target_idx + 2] = 0
plt.plot(xf, np.abs(yf))
plt.show()
由于只有一個(gè)高峰，下面應(yīng)用傅立葉逆變換返回時(shí)域。
應(yīng)用逆FFT與應(yīng)用FFT相似：
from scipy.fft import irfft
new_sig = irfft(yf)
plt.plot(new_sig[:1000])
plt.show()
由于您正在使用rfft()，因此需要使用irfft()來(lái)應(yīng)用反函數(shù)。但是，如果您使用過(guò)fft()，則反函數(shù)將是ifft()。現(xiàn)在，您的繪圖應(yīng)如下所示：
現(xiàn)在有一個(gè)以400 Hz振蕩的正弦波，并且您已經(jīng)成功地消除了4000 Hz的噪聲。
對(duì)信號(hào)進(jìn)行歸一化，然后再將其寫入文件。
norm_new_sig = np.int16(new_sig * (32767 / new_sig.max()))
write("clean.wav", SAMPLE_RATE, norm_new_sig)
到此這篇關(guān)于Python使用scipy.fft進(jìn)行大學(xué)經(jīng)典的傅立葉變換的文章就介紹到這了,更多相關(guān)Python 傅立葉變換內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！
