MATLABで学ぶ信号処理　 窓関数処理後の振幅とパワーについて

数式や画像などが正しく表記できない場合は、下記リンクを参照いただけると幸いです。

https://mech-eng-uni.com/matlabで学ぶ信号処理　-窓関数処理後の振幅とパワー/

---

前回の記事でfftと窓関数の関係についてざっくりと説明しました。本記事では、窓処理がfft結果の振幅やパワースペクトルに及ぼす影響について説明します。また、振幅とパワースペクトルの関係について詳しく説明します。

https://mech-eng-uni.com/matlabで学ぶ信号処理　-窓関数について/

はじめ

上記記事の繰り返しになりましたが、窓関数（ハニング窓）を使うと、信号は下図のようになります。

この窓関数により、リーケージ誤差を低減できるメリットがあります。
しかし、窓関数処理により、本来の信号よりも振幅が小さくなってしまいます。
そのため、窓処理の補正をせずに、そのままfftでの振幅やパワースペクトルを求めると、結果がおかしくなってしまいます。
本記事では、窓処理がfft結果の振幅やパワースペクトルに及ぼす影響について説明します。

窓関数について

一例として、ハニング窓の関数を以下に示します。

振幅

では、上記のハニング窓をすることで、どれだけ振幅が小さくなったのかを考えましょう。
ハニング窓関数の平均値は下式で求まります。とだいたいの教科書や資料には書いてあります。

正直私は数式ばかりで説明するよりも、頭の中のイメージを可視化しながら説明する方がわかりやすいと思うので、あまり数式に頼らずに説明したいと思います。
上式を解かなくても、普通にハニング窓の形状を考えてみてください。
「ハニング窓で振幅が小さくなった分はどのくらい？」と言われたら、ハニング窓処理をしない場合の正方形の面積が1なので、ハニング窓の面積を求めればいいことになります。
この面積は右図のように考えると、0.5というのが容易にわかります。

つまり、ハニング窓処理により振幅が0.5小さくなり、これを補正する必要があります。

パワー

パワーは振幅の二乗なので、窓処理による影響は窓関数の二乗になります。

ハニング窓の二乗は下図になります。
（すいません。これは積分して求めるしかないっぽいです。）
だいたいの窓は正弦波で表現でき、正弦波^2の積分になるので、振幅での補正値の二乗とはならないということを覚えておいてください。

fft後の補正について

fftについては下記でまとめているので、一度ご確認いただけると幸いです。
https://mech-eng-uni.com/matlabで学ぶ信号処理　-高速フーリエ変換（fast-fourier-transformfft）/

振幅

では、実際にMATLABでプログラムを組んで、補正の効果を確認しましょう。100Hzの正弦波に対して、①窓処理なし、②窓処理あり補正なし、③窓処理あり補正ありを比較します。下図がfft結果です。①窓処理なしでは100Hzの振幅が1になっています。一方、②窓処理あり補正なしでは100Hzの振幅が0.5となっています。ただし、補正後の③窓処理あり補正ありでは100Hzの振幅が1になっています。

次に、99Hzおよび101Hzの振幅を確認します。本来であればこの周波数における振幅は0でありますが、誤差程度の振幅を持ちます。①窓処理なしでは振幅が0.02と非常に小さい値です。一方、③窓処理あり補正ありでは振幅が0.5程度であることがわかります。

ここからわかることは、窓処理をすることで最大振幅が正確に把握できるようになる、ということです。また、窓処理により最大振幅近傍の周波数でも振幅値を持つようになることを理解しておく必要があります。

clear all;close all;
smplingFreq=5000;
time=linspace(0,1,smplingFreq);
y=sin(100*2*pi*time); %100Hzの正弦波
hannWin = max(0.5*ones(1,smplingFreq) - 0.5*cos(2.0*pi*time),eps);
% fft
Y1=fft(y,smplingFreq)/length(y)*2;
Y2=fft(y.*hannWin,smplingFreq)/length(y)*2;
% 作図
figure
plot(0:smplingFreq/2,abs(Y1(1:end/2+1)),'ko-')
hold on
plot(0:smplingFreq/2,abs(Y2(1:end/2+1)),'bo-')
plot(0:smplingFreq/2,abs(Y2(1:end/2+1))*2,'ro-')
xlim([90 110])
xlabel('周波数 Hz')
ylabel('振幅')
legend('窓処理なし','窓処理あり-補正なし','窓処理あり-補正あり')

パワースペクトル

振幅と同様に、次はパワースペクトルで補正の効果を確認しましょう。100Hzの正弦波に対して、①窓処理なし、②窓処理あり補正なし、③窓処理あり補正ありを比較します。下図がパワースペクトルの結果です。①窓処理なしでは100Hzのパワーが0.5になっています。一方、補正後の③窓処理あり補正ありでは100Hzの振幅が0.333になっています。これはなぜでしょうか？

これは、パワースペクトルではエネルギーが一致するように窓関数をかけているためです。言い換えればO.A.値などの積算のパワーが一致するようになっています。
下図がパワーの総和をグラフ化したものです。（グラフのx軸の1が①窓処理なし....に対応しています。）

clear all;close all;
smplingFreq=5000;
time=linspace(0,1,smplingFreq);
y=sin(100*2*pi*time); %100Hzの正弦波
hannWin = max(0.5*ones(1,smplingFreq) - 0.5*cos(2.0*pi*time),eps);
% fft
Y1=fft(y,smplingFreq)/length(y);
Y2=fft(y.*hannWin,smplingFreq)/length(y);
% 作図
figure
plot(0:smplingFreq/2,abs(Y1(1:end/2+1)).^2*2,'ko-')
hold on
plot(0:smplingFreq/2,abs(Y2(1:end/2+1)).^2*2,'bo-')
plot(0:smplingFreq/2,abs(Y2(1:end/2+1)).^2*2*8/3,'ro-')
xlim([90 110])
legend('窓処理なし','窓処理あり-補正なし','窓処理あり-補正あり')
xlabel('周波数 Hz')
ylabel('パワー')

figure
bar([sum(abs(Y1(1:end/2+1)).^2*2) sum(abs(Y2(1:end/2+1)).^2*2) sum(abs(Y2(1:end/2+1)).^2*2*8/3)])

※注意！振幅スペクトル^2≠パワースペクトル

これまでのグラフから振幅1の正弦波のfft結果は振幅1となり、パワーは0.5となることがわかりました。教科書などではパワースペクトルは振幅の二乗という記載がされています。教科書の数式展開には間違いはありません。そう考えると、パワーは1ではないのか？と感じる方が多いと思います。

しかし、パワーは0.5が正解です。
この問題は振動・騒音を研究している方でも誤解している人がたくさんいるのでここで説明しておこうと思います。

例えば、教科書などでは下記のように記載されています。
"時間信号をs(t)として、この信号を周波数次元に変換したものをS(ω)と表します。ここで、s(t)のパワースペクトルをP(ω)とすると、P(ω)=S(ω)^2となります。"
ここで問題なのが、周波数変換したときの負の周波数成分の補正をいつするのか！ということです。
下記の記事で述べたようにfftでスペクトルを求めるときは、負の周波数成分（虚像）を折り返して共役で足し合わせる必要があります。

https://mech-eng-uni.com/matlabで学ぶ信号処理　-高速フーリエ変換（fast-fourier-transformfft）/

この計算をいつ行うのか！ということが非常に問題です。
本記事の「振幅」を求めるときはプログラム上でfftをする際に補正しています。（2倍して補正）
では、パワースペクトルを求める場合はどうでしょうか？
このときは下記の2パターンが考えられます。
A1：fft結果を2倍して、「振幅」にして、「振幅」を二乗してパワースペクトルを求める
A2：fft結果を二乗して、その結果を2倍して、パワースペクトルを求める

これはA2が正解です。
なぜなら、単位時間当たりの平均エネルギーであるパワースペクトル密度には下記関係があるからです。

上記式から、時間信号を二乗してものを積分していることがわかります。
つまり、周波数次元でも二乗した後に、負の周波数成分（虚像）を折り返して共役で足し合わせる必要があります。
上記は非常に重要です。
この理解ができていないと、振幅を分析したあとにパワーやPSDを分析するときなどに、必ず間違えます。
また騒音を分析するときも、音圧[Pa]で分析するときと、音圧[dB]で分析するときで間違えが生じます。
ここらへんはまた詳しく説明するかもしれません。

今回はこのへんでGood luck

私の個人ブログもよろしくお願いします。
https://mech-eng-uni.com/