フーリエ逆変換(窓関数の効果) | python3Xのブログ

python3Xのブログ

ここでは40代、50代の方が日々の生活で役に立つ情報や私の趣味であるプログラム、Excelや科学に関する内容で投稿する予定です。

テーマ:

今回は

Sampling rate: 4000

Sampling length: 512

Frequency: 44Hz, 60Hz(Amplitude: x 2), 80Hz(Amplitude: x 3)

Noise: np.random.randn(L)

カットオフ周波数:300

カットオフ振幅強度:0.2

Window無しとhamming window、hanning window、bartlett window、blackman window

の窓関数を用いた場合の比較を行いました

ノイズ除去に関しては、窓関数の効果は認められませんでした。

 

●周波数によるカット

●振幅の大きさによるカット

#coding UTF-8
from
#coding UTF-8
from matplotlib import pyplot as plt
import numpy as np
fs = 4000
L = 512  
# 440Hzのsin波を作る
sin44 = np.sin(2. * np.pi * np.arange(L) * 44. / fs)
# 600Hzのsin波を作る
sin60= 2 * np.sin(2. * np.pi * np.arange(L) * 60. / fs)
# 800Hzのsin波を作る
sin80 = 3 * np.sin(2. * np.pi * np.arange(L) * 80. / fs)
# 3つを足し合わせる
sinall = sin44 + sin60 + sin80 + 0.3 * np.random.randn(L)
corect = sin44 + sin60 + sin80

# 窓関数:ハミング窓(0→最大→0に近付くイメージ)
hammingWindow = np.hamming(L)
hanningWindow = np.hanning(L)
bartlettWindow = np.bartlett(L)
blackmanWindow = np.blackman(L)
# フーリエ変換
spect_nW = np.fft.fft(sinall) # 窓関数なし
hammingWin = np.fft.fft(sinall * hammingWindow)
hanningWin = np.fft.fft(sinall * hanningWindow)
bartlettWin = np.fft.fft(sinall * bartlettWindow)
blackmanWin = np.fft.fft(sinall * blackmanWindow)
spect_nw = abs(spect_nW)
hamming = abs(hammingWin)
hanning = abs(hanningWin)
bartlett = abs(bartlettWin)
blackman = abs(blackmanWin)
freq = np.linspace(0, fs, L)
f1 = np.copy(2*spect_nW )
f2 = np.copy(2*hammingWin)
f3 = np.copy(2*hanningWin)
f4 = np.copy(2*bartlettWin)
f5 = np.copy(2*blackmanWin)
fc = 300
f1[freq > fc] = 0
f2[freq > fc] = 0
f3[freq > fc] = 0
f4[freq > fc] = 0
f5[freq > fc] = 0
f1_abs = np.abs(f1) / 2
f2_abs = np.abs(f1) / 2
f3_abs = np.abs(f1) / 2
f4_abs = np.abs(f1) / 2
f5_abs = np.abs(f1) /2
# フーリエ逆変換
inv_f1 = np.fft.ifft(f1)
inv_f2 = np.fft.ifft(f2)
inv_f2 /= hammingWindow
inv_f3 = np.fft.ifft(f3)
inv_f3/= hanningWindow
inv_f4 = np.fft.ifft(f4)
inv_f4 /= bartlettWindow
inv_f5 = np.fft.ifft(f5)
inv_f5 /= blackmanWindow
f1_real = inv_f1.real
f2_real = inv_f2.real
f3_real = inv_f3.real
f4_real = inv_f4.real
f5_real = inv_f5.real
f11 = np.copy(spect_nW )
f12 = np.copy(hammingWin)
f13 = np.copy(hanningWin)
f14 = np.copy(bartlettWin)
f15 = np.copy(blackmanWin)
ac = 0.2
f11[f1_abs < ac] = 0
f12[f1_abs < ac] = 0
f13[f1_abs < ac] = 0
f14[f1_abs < ac] = 0
f15[f1_abs < ac] = 0
inv_f11 = np.fft.ifft(2*f11)
inv_f12 = np.fft.ifft(2*f12)
inv_f12 /= hammingWindow
inv_f13 = np.fft.ifft(2*f13)
inv_f13/= hanningWindow
inv_f14 = np.fft.ifft(2*f14)
inv_f14 /= bartlettWindow
inv_f15 = np.fft.ifft(2*f15)
inv_f15 /= blackmanWindow

# 図を表示
fig = plt.figure(figsize=(20, 70))
plt.rcParams["font.size"] = 20
fig.add_subplot(511)
plt.plot(sinall, label="ノイズあり")
plt.plot(corect, label="ノイズなし")
plt.xlim([0, L])
plt.ylim([-6, 6])
plt.xlabel("$time$")
plt.ylabel("$signal$")
plt.grid(True)
plt.legend(loc='upper right')
plt.title("1. 入力信号")
fig.add_subplot(512)
plt.plot(inv_f1, label="窓関数なし")
plt.plot(inv_f2, label="ハミング窓(周波数でカット)")
plt.xlim([0, L])
plt.ylim([-6, 6])
plt.xlabel("$time$")
plt.ylabel("$signal$")
plt.grid(True)
plt.legend(loc='upper right')
plt.title("2. 比較:ハミング窓(周波数でカット)")
fig.add_subplot(513)
plt.plot(inv_f1, label="窓関数なし")
plt.plot(inv_f3, label="ハン窓")
plt.xlim([0, L])
plt.ylim([-6, 6])
plt.xlabel("$time$")
plt.ylabel("$signal$")
plt.grid(True)
plt.legend(loc='upper right')
plt.title("3. 比較:ハン窓(周波数でカット)")
fig.add_subplot(514)
plt.plot(inv_f1, label="窓関数なし")
plt.plot(inv_f4, label="バートレット窓")
plt.xlim([0, L])
plt.ylim([-6, 6])
plt.xlabel("$time$")
plt.ylabel("$signal$")
plt.grid(True)
plt.legend(loc='upper right')
plt.title("4. 比較:バートレット窓(周波数でカット)")
fig.add_subplot(515)
plt.plot(inv_f1, label="窓関数なし")
plt.plot(inv_f5, label="ブラックマン窓")
plt.xlim([0, L])
plt.ylim([-6, 6])
plt.xlabel("$time$")
plt.ylabel("$signal$")
plt.grid(True)
plt.legend(loc='upper right')
plt.title("5. 比較:ブラックマン窓(周波数でカット)")
plt.subplots_adjust(hspace=0.5) # グラフ間の隙間調整
fig = plt.figure(figsize=(20, 70))
plt.rcParams["font.size"] = 20
fig.add_subplot(511)
plt.plot(sinall, label="ノイズあり")
plt.plot(corect, label="ノイズなし")
plt.xlim([0, L])
plt.ylim([-6, 6])
plt.xlabel("$time$")
plt.ylabel("$signal$")
plt.grid(True)
plt.legend(loc='upper right')
plt.title("1. 入力信号")
fig.add_subplot(512)
plt.plot(inv_f11, label="窓関数なし")
plt.plot(inv_f12, label="ハミング窓(周波数でカット)")
plt.xlim([0, L])
plt.ylim([-6, 6])
plt.xlabel("$time$")
plt.ylabel("$signal$")
plt.grid(True)
plt.legend(loc='upper right')
plt.title("2. 比較:ハミング窓(振幅の大きさでカット)")
fig.add_subplot(513)
plt.plot(inv_f11, label="窓関数なし")
plt.plot(inv_f13, label="ハン窓")
plt.xlim([0, L])
plt.ylim([-6, 6])
plt.xlabel("$time$")
plt.ylabel("$signal$")
plt.grid(True)
plt.legend(loc='upper right')
plt.title("3. 比較:ハン窓(振幅の大きさでカット)")
fig.add_subplot(514)
plt.plot(inv_f11, label="窓関数なし")
plt.plot(inv_f14, label="バートレット窓")
plt.xlim([0, L])
plt.ylim([-6, 6])
plt.xlabel("$time$")
plt.ylabel("$signal$")
plt.grid(True)
plt.legend(loc='upper right')
plt.title("4. 比較:バートレット窓(振幅の大きさでカット)")
fig.add_subplot(515)
plt.plot(inv_f11, label="窓関数なし")
plt.plot(inv_f15, label="ブラックマン窓")
plt.xlim([0, L])
plt.ylim([-6, 6])
plt.xlabel("$time$")
plt.ylabel("$signal$")
plt.grid(True)
plt.legend(loc='upper right')
plt.title("5. 比較:ブラックマン窓(振幅の大きさでカット))")
plt.subplots_adjust(hspace=0.5) # グラフ間の隙間調整
plt.show()