Q&A
窓関数に付いては、
fft 窓 不連続
でググって見つかったwebページを見てください
質問
現在、PythonのscipyライブラリのFFTを使用し、計測機器から得られたデータを解析しています。
その中で、下記のサイト様を参考にFFT解析を行っています。
手順として、①オーバーラップ処理、②ハニング窓、③FFT、④③の平均化という流れで解析を行っているのですが、一つ疑問として、なぜ①、②の工程が必要なのでしょうか。
データをすべてscipy.fft.fft()のライブラリに代入してしまえば、似たような波形が出てくると思っているのですが、、、
[質問①] 何か理由があれば、ご教示いただきたいです。
[質問②] またscipyとnumpyのFFTの処理は何か違う点等あるのでしょうか。
参考URL：https://watlab-blog.com/2019/04/21/python-fft/
実現したいこと
また現実したいこととしては、FFT変換後、任意の周波数帯だけ指定し、
IFFT変換したいと考えています。
[質問③] もし、上記のみで(③の処理のみ)問題のない場合、以下のような手順でFFT、IFFT変換を行えばよかったでしょうか。またこの場合、dataとdata_ifftは瓜二つのデータ群になるのでしょうか。
Python
1data = np.read_csv(..., encoding="shift-jis") 2 3data_fft = numpy.fft.fft(data) 4data_ifft = numpy.fft.ifft(data_fft)
参考ページ内の｢PythonでFFTをする前にオーバーラップ処理をしよう！｣( https://watlab-blog.com/2019/04/17/python-overlap/ )を読んで理解出来なかったところを質問に追記ください
>データをすべてscipy.fft.fft()のライブラリに代入してしまえば、似たような波形が出てくると思っているのですが、、、
と記載の通り、違いがなかったので何が違うのか、なぜこのような手順を踏まなければいけいないのかが分かっていません。。オーバーラップ処理や、それにより必要となってくるハニング窓関数の処理などは分かったのですが、これらの処理をしたものとscipy.fft.fft()のライブラリに入れたものとの違いが判らなかったので、そちらについてご教授をお願いいただければと思っております。よろしくお願いいたします。
> 違いがなかった
下記を実行してみてください
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft
x = np.arange(0, 2*np.pi*8.5, np.pi/180)
y = np.sin(x)
h_window = np.hamming(len(y))
yw = y * h_window
plt.plot(x, y, label='no window')
plt.plot(x, yw, label='window')
plt.legend()
plt.show()
yf = np.abs(fft(y)) * 2 / len(y)
ywf = np.abs(fft(yw)) * 2 / np.sum(h_window)
plt.plot(yf, label='no window')
plt.plot(ywf, label='window')
plt.xlim(0, 16)
plt.grid()
plt.legend()
plt.show()
元データが8.5周期なので、fft後の結果は横軸8, 9以外は小さくなってほしいですよね
窓関数無しだと、横軸8, 9から離れたところも、値がまあまああります
窓関数有りだと、横軸8, 9と、そのすぐ近くを除くと、値は0に近いです
jbpb0様、ご回答ありがとうございます。
窓関数を使用しないと、必要な周波数帯以外の部分もFFT処理で計測されてしまうことが分かりました。
理解しました。ありがとうございます。
また回答いただきましたプログラムで1点追加で質問させていただきます。
yf = np.abs(fft(y)) * 2 / len(y)
ywf = np.abs(fft(yw)) * 2 / np.sum(h_window)
yfは、他の質問で積分誤差を減らすために(2/len(y))を割るとの記事を見たのですが、
ywfはなぜ2 / np.sum(h_window)で割るのでしょうか。
自分の中では、ハミング窓による減衰を相殺する目的という認識であってますでしょうか。
何度も申し訳ございませんが、ご回答のほどよろしくお願いいたします。
追記いたします。私の質問より、参考URLのサイトでは、
ywf = np.abs(fft(yw)) * 2 / np.sum(h_window)の”2 / np.sum(h_window)”部分が
”1 / np.mean(h_window)”となっていたのですが、
これはどちらの認識があっていますでしょうか。
> ハミング窓による減衰を相殺する目的という認識であってますでしょうか。
はい
> ”2 / np.sum(h_window)”部分が”1 / np.mean(h_window)”となっていた
フーリエ変換の正規化は、フーリエ変換→フーリエ逆変換 で元に戻ればいいのであって、フーリエ変換だけの正規化には、どれが正解というのはありません
あくまでも、フーリエ変換の正規化と、フーリエ逆変換の正規化がセットで意味があるのです
私のコメントのコードで、
x = np.arange(0, 2*np.pi*8.5, np.pi/180)
↓ 変更
x = np.arange(0, 2*np.pi*8, np.pi/180)
だけ変えて、他はそのままで実行してみてください
元のデータの周期が8なので横軸8の値を見ると、窓関数有り無しどちらも1.0ですよね
これは、元のデータの振幅±1.0に合わせて表示されるように正規化したからです
次に、
yf = np.abs(fft(y)) * 2 / len(y)
ywf = np.abs(fft(yw)) * 2 / np.sum(h_window)
↓ 変更
yf = np.abs(fft(y))
ywf = np.abs(fft(yw)) / np.mean(h_window)
も変更して、実行してみてください
横軸8の値は、1.0ではないけど、窓関数有り無しでやはり一致してますよね
なので、これはこれで正しいのです
繰り返しますが、フーリエ変換→フーリエ逆変換 で元に戻ればいいのであって、フーリエ変換だけの正規化には、どれが正解というのはありません
なので、質問者さんの目的に合うように正規化してください
なるほど。FFTのライブラリで出力されるのが比率で計算結果が出力されるのですね。
自分の目的にあった処理をするようにいたします。
jbpb0様の回答で疑問を解消することができました。ありがとうございました。
回答2件
0
自己解決
jbpb0様のコメントで疑問を解消しました。
投稿2023/04/17 04:46
総合スコア7
0
なぜ必要なのか疑問なら、一度それらの処理を省いてFFTを掛けるなりしてみよう
どうなりますか？
投稿2023/04/13 06:32
総合スコア88173
あなたの回答
tips
プレビュー
