Python使用scipy.fft進(jìn)行大學(xué)經(jīng)典的傅立葉變換

傅里葉變換是在高數(shù)是一個很重要的知識點(diǎn)，今天將結(jié)合Python代碼實(shí)現(xiàn)傅立葉變換。

傅立葉變換

我們平時是如何去分解一個復(fù)雜的問題呢？一個經(jīng)典的方法就是把這個復(fù)雜的問題分解成為多個簡單的可操作的子問題， 傅立葉變換也是基于這個思想。

傅里葉分析是研究如何將數(shù)學(xué)函數(shù)分解為一系列更簡單的三角函數(shù)的領(lǐng)域。傅里葉變換是該領(lǐng)域的一種工具，用于將函數(shù)分解為其分量頻率。

在本教程中，傅立葉變換是一種工具，可以獲取信號并查看其中每個頻率的功率?？匆豢丛摳盗⑷~變換中的重要術(shù)語：

• 信號：信號是隨時間變化的信息。例如，音頻，視頻和電壓走線都是信號的示例。
• 頻率：頻率是某物重復(fù)的速度。例如，時鐘以1赫茲（Hz）的頻率滴答，或每秒重復(fù)1次。
• 功率：功率表示每個頻率的強(qiáng)度。

下圖是一些正弦波的頻率和功率的直觀演示：

第一個是低頻正弦波，第二個是高頻正弦波，第三個是低頻低功率正弦波，因此低功率正弦波比其它兩個正弦波的峰較小。

時域與頻域

時域與頻域是查看信號的兩種不同方式，即信號的組成頻率或隨時間變化的信息。

在時域中，信號是隨時間（x軸）幅度（y軸）變化的波。您最有可能在時域中查看圖表，例如：

這是一些音頻的圖像，它是一個時域信號。橫軸表示時間，縱軸表示振幅。

在頻域中，信號表示為一系列頻率（x軸），每個頻率都具有關(guān)聯(lián)的功率（y軸）。下圖是經(jīng)過傅立葉變換后的上述音頻信號：

代碼實(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軸以秒為單位表示時間，并且由于每秒鐘的時間都有兩個峰值，因此可以看到正弦波每秒振蕩兩次。

混合音頻

下面將兩個正弦波，混合音頻信號僅包括兩個步驟：

將正弦波加在一起，然后進(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

下一步是歸一化，或縮放信號以適合目標(biāo)格式。由于以后將如何存儲音頻，目標(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)化為音頻，最簡單的方法是使用SciPy的wavfile.write方法將其存儲在WAV文件中。16位整數(shù)是WAV文件的標(biāo)準(zhǔn)數(shù)據(jù)類型，因此需要將信號標(biāo)準(zhǔn)化為16位整數(shù)：

from scipy.io.wavfile import write

# 記住，采樣率=44100赫茲是我們的播放率
write("mysinewave.wav", SAMPLE_RATE, normalized_tone)

這個音頻聽起來音調(diào)很高。

完成此步驟后，就當(dāng)作音頻樣本了。下一步是使用傅立葉變換消除高音調(diào)！

傅立葉變換

現(xiàn)在對生成的音頻上使用FFT了。FFT是一種算法，可實(shí)現(xiàn)傅立葉變換并可以在時域中為信號計(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()

我們可以在正頻率中看到兩個峰值，正頻率峰值位于400 Hz和4000 Hz，與之前生成的音頻的頻率相對應(yīng)。

計(jì)算傅里葉變換

yf = fft(normalized_tone)
xf = fftfreq(N, 1 / SAMPLE_RATE)

上面代碼的功能

• fft() 計(jì)算轉(zhuǎn)換本身。
• fftfreq()計(jì)算的輸出中每個倉中心的頻率fft()。沒有這個，就無法在頻譜上繪制x軸

fft()輸出的頻譜圍繞y軸反射，因此負(fù)半部分是正半部分的鏡像，我們一般只需計(jì)算一半對稱值，即可更快地進(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()

過濾信號

傅里葉變換的一大優(yōu)點(diǎn)是它是可逆的，我們可以利用此優(yōu)勢來過濾音頻并擺脫高音調(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來擺脫它：

yf[target_idx - 1 : target_idx + 2] = 0

plt.plot(xf, np.abs(yf))
plt.show()

由于只有一個高峰，下面應(yīng)用傅立葉逆變換返回時域。

應(yīng)用逆FFT與應(yīng)用FFT相似：

from scipy.fft import irfft

new_sig = irfft(yf)

plt.plot(new_sig[:1000])
plt.show()

由于您正在使用rfft()，因此需要使用irfft()來應(yīng)用反函數(shù)。但是，如果您使用過fft()，則反函數(shù)將是ifft()?，F(xiàn)在，您的繪圖應(yīng)如下所示：

現(xiàn)在有一個以400 Hz振蕩的正弦波，并且您已經(jīng)成功地消除了4000 Hz的噪聲。

對信號進(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)容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！