前回、観測期間の中で整数にならない周波数を FFT に掛けたときの応答を紹介しました。それだけ見ると「S / N が悪そうだね～」で終わってしまうのですが、他にも見たい周波数があった場合はどうでしょうか。たとえば、家庭での電灯線の周波数は 50Hz / 60Hz です。で、これに含まれているもっと小さな高調波成分の有無、量を調べたいと思ったとします。観測期間を十分に長くするとか、HPF でそれらの周波数を落とすとかが容易に出来れば問題ないのですが（こちらを参考にして下さい）、たとえば 50ms ぐらいしか観測期間が与えられなくて（50Hz / 60Hz のいずれに対しても端数になる）、その中に含まれている 100Hz とか 200Hz の成分を見ておきたい、などという場合どういうことが起きるでしょうか。
今回は、50Hz / 60Hz を 2 ～ 3Hz だと思って、それに振幅の小さい 10Hz が加わったという前提で FFT を掛けてみましょう。
スクリプトです。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
// phi=int(0/100*N) // 位相をずらす量をパーセントで記入。
i=0
j=1
fs=2
fs1=10 // 観測したい信号の周波数
amp1=0.1 // 観測したい信号の振幅
for i=0:0.25:1
fs=2+i // 端数を入れる
S_sin(j,:)=sin(2*%pi*t*fs)+amp1*sin(2*%pi*t*fs1)
f_sin(j,:)=fft(S_sin(j,:)) // 関数 fft を呼び出す。
M_sin(j,:)=abs(f_sin(j,:))*2 ./N // 振幅を計算
disp(max(M_sin(j,:)))
j=j+1
end
fig_axis=[0,1,-2,2]
scf(10);clf;plot(t',S_sin')
title('2.0 - 3.0Hz + 10Hz signal','fontsize',4)
mtlb_axis(fig_axis)
legend('f=2.0','f=2.25','f=2.5','f=2.75','f=3.0')
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=20
ft_axis=[0,fmax,-60,10]
scf(20);plot([0:fmax]',20.*log10(M_sin(1:j-1,1:fmax+1)'))
title('2.0 - 3.0 Hz + 10Hz Spectrum','fontsize',4)
mtlb_axis(ft_axis)
legend('f=2.0','f=2.25','f=2.5','f=2.75','f=3.0')
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
波形はこうです。
FFT の結果はこうです。
まあ、微妙ですが 10Hz が検出できてなくもないというところでしょうか。もちろん、2Hz、3Hz の時はしっかり分離して測定できています。振幅自体は 0.1 なので -20dB ぐらいということでまあまあかなです。
ところが、観測したい信号の振幅を 0.01（-40dB）にするとこうなってしまいます。（amp1=0.01）
ちょっと極端な例かも知れませんが、波形を見ると 2 - 3 Hz だけがしっかり見えていて、10Hz はいるのかな？という信号ですが、観測期間がきちっとしていると 10Hz が測定できて余分な周波数がいる、ということが分かります。しかし観測期間が悪いと存在がはっきりしなくなります。要は中心となる周波数に対して観測期間はきちっと合わせておいた方が、それ以外の成分が見やすいということになります。
ちなみに、10Hz を 10.5Hz にして振幅を 0.1 に戻すとこうなります。
なんかすごいことになってしまいました。
こういうことが起きるのは観測期間の波形がそのまま繰り返されるという前提で FFT が計算するので、多少加工してでも繰り返し性をよく見せる必要があります。これについては次回取り上げます。
今回は、50Hz / 60Hz を 2 ～ 3Hz だと思って、それに振幅の小さい 10Hz が加わったという前提で FFT を掛けてみましょう。
スクリプトです。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
// phi=int(0/100*N) // 位相をずらす量をパーセントで記入。
i=0
j=1
fs=2
fs1=10 // 観測したい信号の周波数
amp1=0.1 // 観測したい信号の振幅
for i=0:0.25:1
fs=2+i // 端数を入れる
S_sin(j,:)=sin(2*%pi*t*fs)+amp1*sin(2*%pi*t*fs1)
f_sin(j,:)=fft(S_sin(j,:)) // 関数 fft を呼び出す。
M_sin(j,:)=abs(f_sin(j,:))*2 ./N // 振幅を計算
disp(max(M_sin(j,:)))
j=j+1
end
fig_axis=[0,1,-2,2]
scf(10);clf;plot(t',S_sin')
title('2.0 - 3.0Hz + 10Hz signal','fontsize',4)
mtlb_axis(fig_axis)
legend('f=2.0','f=2.25','f=2.5','f=2.75','f=3.0')
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=20
ft_axis=[0,fmax,-60,10]
scf(20);plot([0:fmax]',20.*log10(M_sin(1:j-1,1:fmax+1)'))
title('2.0 - 3.0 Hz + 10Hz Spectrum','fontsize',4)
mtlb_axis(ft_axis)
legend('f=2.0','f=2.25','f=2.5','f=2.75','f=3.0')
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
波形はこうです。
FFT の結果はこうです。
まあ、微妙ですが 10Hz が検出できてなくもないというところでしょうか。もちろん、2Hz、3Hz の時はしっかり分離して測定できています。振幅自体は 0.1 なので -20dB ぐらいということでまあまあかなです。
ところが、観測したい信号の振幅を 0.01（-40dB）にするとこうなってしまいます。（amp1=0.01）
ちょっと極端な例かも知れませんが、波形を見ると 2 - 3 Hz だけがしっかり見えていて、10Hz はいるのかな？という信号ですが、観測期間がきちっとしていると 10Hz が測定できて余分な周波数がいる、ということが分かります。しかし観測期間が悪いと存在がはっきりしなくなります。要は中心となる周波数に対して観測期間はきちっと合わせておいた方が、それ以外の成分が見やすいということになります。
ちなみに、10Hz を 10.5Hz にして振幅を 0.1 に戻すとこうなります。
なんかすごいことになってしまいました。
こういうことが起きるのは観測期間の波形がそのまま繰り返されるという前提で FFT が計算するので、多少加工してでも繰り返し性をよく見せる必要があります。これについては次回取り上げます。