DFT は時間方向に離散なら周波数方向にも離散的です。
そのため信号の観測期間が有限なら、測定できる周波数分解能も限られてしまうので、今回はここをもう少し詳しく見てみます。
前に出てきた記事をもう一度見て欲しいのですが、
FFT 使い手 Level 5 にありがちなこと（時間制限があるとは？）
http://blogs.yahoo.co.jp/susanoo2001_hero/12333034.html?type=folderlist
ここにも再掲しますが、たとえば 1 秒間観測した画面の中に 2.4Hz の信号があったとして、こんな風に見えたとします。
これを FFT に掛けると云うことは実はこのような連続波であることを前提に分析していると云うことなんです。
これでは本来観測したい波形とは違ってしまっていますね。
実際に 2Hz と 3Hz の間で中途半端な周波数を入れるとどんなスペクトラムが観測されるか試してみます。
スクリプトはこれになります。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
// phi=int(0/100*N) // 位相をずらす量をパーセントで記入。
i=0
j=1
// 基本周波数を入れる
fs=2
for i=0:0.25:1
fs=fs+i // 周波数を端数にする。
S_sin(j,:)=sin(2*%pi*t*fs)
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 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 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 の結果はこうです。
きちっと整数になっている、2Hz と 3Hz はその周波数だけがきちっとスペクトラムが立っていますが、2.2Hz, 2.5Hz, 2.75Hz は広がりを持ってしまっています。
ちなみに 10 - 11 Hz の場合はこうなります。最初の方の fs=2 を fs=10 にするだけです。
だいたい周波数軸方向にずれたようになりますね。
DC 成分だけについて考察しておくと、観測期間の中の平均値（積分値）がいくつになっているかということです。周波数が整数ならばゼロですが、それ以外の場合はなんらかの値が残ります。この値は周波数が低くて端数が 0.5 の時に大きくなると云うのは直感的に分かると思います。
グラフだけ見ると、「2.5Hz、10.5Hz は 2 と 3、10 と 11 が同じぐらいの値になっているようだし、その他のは 2 または 3、10 または 11 に丸めるでいいでねぇの？」と思いたくなりますし、そこだけ見ればそれでいいのですがそういうわけにもいかないということを次回やってみたいと思います。
そのため信号の観測期間が有限なら、測定できる周波数分解能も限られてしまうので、今回はここをもう少し詳しく見てみます。
前に出てきた記事をもう一度見て欲しいのですが、
FFT 使い手 Level 5 にありがちなこと（時間制限があるとは？）
http://blogs.yahoo.co.jp/susanoo2001_hero/12333034.html?type=folderlist
ここにも再掲しますが、たとえば 1 秒間観測した画面の中に 2.4Hz の信号があったとして、こんな風に見えたとします。
これを FFT に掛けると云うことは実はこのような連続波であることを前提に分析していると云うことなんです。
これでは本来観測したい波形とは違ってしまっていますね。
実際に 2Hz と 3Hz の間で中途半端な周波数を入れるとどんなスペクトラムが観測されるか試してみます。
スクリプトはこれになります。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
// phi=int(0/100*N) // 位相をずらす量をパーセントで記入。
i=0
j=1
// 基本周波数を入れる
fs=2
for i=0:0.25:1
fs=fs+i // 周波数を端数にする。
S_sin(j,:)=sin(2*%pi*t*fs)
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 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 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 の結果はこうです。
きちっと整数になっている、2Hz と 3Hz はその周波数だけがきちっとスペクトラムが立っていますが、2.2Hz, 2.5Hz, 2.75Hz は広がりを持ってしまっています。
ちなみに 10 - 11 Hz の場合はこうなります。最初の方の fs=2 を fs=10 にするだけです。
だいたい周波数軸方向にずれたようになりますね。
DC 成分だけについて考察しておくと、観測期間の中の平均値（積分値）がいくつになっているかということです。周波数が整数ならばゼロですが、それ以外の場合はなんらかの値が残ります。この値は周波数が低くて端数が 0.5 の時に大きくなると云うのは直感的に分かると思います。
グラフだけ見ると、「2.5Hz、10.5Hz は 2 と 3、10 と 11 が同じぐらいの値になっているようだし、その他のは 2 または 3、10 または 11 に丸めるでいいでねぇの？」と思いたくなりますし、そこだけ見ればそれでいいのですがそういうわけにもいかないということを次回やってみたいと思います。