numpy5傅立葉變換示例
阿新 • • 發佈：2018-12-14
目錄
介紹
傅立葉這一張如果之後會用到，我會在進行擴充套件，這裡先講解最基本的用法。
定義
如我講解有問題：希望各位指出。
這裡有一個介紹傅立葉變換的部落格，講得很好，如果不知道傅立葉變換是什麼，可以在這裡看。https://www.cnblogs.com/h2zZhou/p/8405717.html
關鍵字解釋：
傅立葉變換和反傅立葉變換是將聲音在時域和頻域之間的轉換。
時域：以時間作為參照來觀察動態世界的方法我們稱其為時域分析。
頻域：以震動頻率為參照。
原文中以音樂和樂譜的關係來講解時域和頻域，其實是很容易讓人誤解的，因為樂譜畢竟也是按照時間先後順序進行排列的。
但是頻域不同，頻域的概念更像你在同一時刻把樂譜上所有的聲音都彈出來，而彈出來之後，聲音會隨著時間的變化自動變成你需要的音樂。因為兩者之間的音符不同，頻域的音符是連續的，是從頭貫穿到尾的，但是時域的音符則是間斷的，是隻有那一時刻才有的。
很難理解，所以說文中這麼形容頻域：世界是永恆不變的
在頻域的世界裡，音樂在播放的第一時刻就已經決定了後面的所有的形容。
快速傅立葉變換模組（fft）
傅立葉定理：任何一個周期函式總可以被分解為若干不同振幅、頻率和相位的正弦函式的疊加。
f(x) = 4/1pi sin(1x) + 4/3pi sin(3x) + 4/5pi sin(5x) +4/7pi sin(7x) + …
非周期函式可被視為週期無窮大的周期函式，傅立葉級數變成傅立葉積分。在工程上，從這些連續的頻率中再選取離散的頻率取樣，得到離散化的頻率序列，這就是離散傅立葉形式。
時間域的離散 -> 頻率域的離散 x1->y1 f1 -> A1 fai1 A1sin(f12pi x+fai1) x2->y2 f2 -> A2 fai2 A2sin(f22pi x+fai2) ... ... xn->yn fn -> An fain Ansin(fn2pi x+fai2) \_________FFT變換_________^ +)------------------------ y = f(x) ^_______________IFFT反變換_________________/ import np.fft as nf nf.fftfreq(時間域離散樣本的個數，時間域離散樣本的間隔（即取樣週期）) 取樣週期=1/取樣頻率 ->頻率域的離散樣本系列(f1,f2,...,fn)，單位Hz nf.fft([y1,y2,...,yn])->[(A1,fai1), (A2,fai2), ...] \ / a+bj 模正比A1 幅角就是fai1 nf.ifft([(A1,fai1), (A2,fai2), ...])->[y1,y2,...,yn]
import numpy as np
import numpy.fft as nf
import matplotlib.pyplot as mp
times = np.linspace(0, 2 * np.pi, 201)
#畫五根曲線
sigs1 = 4 / (1 * np.pi) * np.sin(1 * times)
sigs2 = 4 / (3 * np.pi) * np.sin(3 * times)
sigs3 = 4 / (5 * np.pi) * np.sin(5 * times)
sigs4 = 4 / (7 * np.pi) * np.sin(7 * times)
sigs5 = 4 
例項：音訊濾波
.wav檔案記錄一段聲音在每一個取樣時刻的位移，即位移關於時間的函式：s=f(t)，以時間域的方式表現聲音。
數學中的函式：y = 3x^2+2x+1
計算機中的函式：多個存在對應關係的陣列
1 2 3 4 5 …
6 17 34 57 86 …
含有噪聲的 FFT 含有噪聲的 頻率濾波 不含噪聲的
時間域訊號 -------> 頻率域譜線 ----------> 頻率域譜線
# -*- coding: utf-8 -*-
from __future__ import unicode_literals
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(
r'C:\Users\Cs\Desktop\資料分析\DS+ML\DS\data\noised.wav')
print(sample_rate,noised_sigs)
#儲存的時候放大了2**15，（不然有些浮點數太小會直接被忽略？）所以這裡先還原數值
noised_sigs = noised_sigs / 2 ** 15
#記錄的次數，每個之間跨度為1/比率
times = np.arange(len(noised_sigs)) / sample_rate
print("times:",times)
print("rate:",sample_rate)
#頻率序列
freqs = nf.fftfreq(times.size, 1 / sample_rate)
#進行傅立葉變換
noised_ffts = nf.fft(noised_sigs)
#獲得複數的模，即振幅絕對值
noised_pows = np.abs(noised_ffts)
#armax()返回最大值的下標，將能量值較小的噪聲過濾掉
fund_freq = np.abs(freqs[noised_pows.argmax()])
print("fund_freq:",fund_freq)
#生成篩子陣列，這裡是噪聲的下標
noised_indices = np.where(np.abs(freqs) != fund_freq)
# 過濾噪聲
filter_ffts = noised_ffts.copy()
filter_ffts[noised_indices] = 0
#過濾噪聲之後的聲音的模
filter_pows = np.abs(filter_ffts)
#逆傅立葉變換，得到的是類似： \
# 0.02333286-3.29975402e-17j格式的資料，這裡過濾掉虛數部分，為啥要過濾掉虛數？？
filter_sigs = nf.ifft(filter_ffts).real
# 儲存過濾後的資料
wf.write(r'C:\Users\Cs\Desktop\資料分析\DS+ML\DS\data\filter2.wav', sample_rate,
(filter_sigs * 2 ** 15).astype(np.int16))
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.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()
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()
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()
mp.tight_layout()
mp.show()
