妙に端数の周波数にこだわるな、と感じているかも知れませんが、実際に限られた観測時間内で整数である周波数の方が確率は低いでしょう。なので、今回は端数の周波数と位相の関係について着目してみます。
カイザー窓を使うと矩形窓よりも端数の周波数に対して S / N が良くなると云うことを確認しました。では、その端数の周波数が観測時間内に色々な位相で入って来たら（普通はそうなる）FFT の結果はどうなるのか見てみることにします。
スクリプトです。だんだん scilab 自体へのレベルが上がってきたのか、今まで使っていなかった関数も結構ありますので、これなんだ？と思った方はヘルプを見て下さい。ここで質問してくれてもいいです。
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
n=1:5
phi(n)=(n-1)*45 // 正弦波の位相: 0deg, 45deg, 90deg, 135deg, 180deg
deg=string(phi)
fs1=2.5 // 信号周波数、後半では10.5Hz に書き換える。
amp1=1 // 信号周波数の振幅
frequency=strcat([string(fs1),'Hz'])
Kaiser=window('kr',N,%pi*2) // Kaiser窓
S_sin=zeros(5,2048);
for i=1:5
S_sin(i,:)=amp1*sin(2*%pi*(t*fs1+phi(i)/360))
WS_sin(i,:)=Kaiser.*S_sin(i,:)
f_sin(i,:)=fft(S_sin(i,:)) // 関数 fft を呼び出す。
M_sin(i,:)=abs(f_sin(i,:))*2/N // 振幅を計算
wf_sin(i,:)=fft(WS_sin(i,:))
WM_sin(i,:)=abs(wf_sin(i,:))*2/N
end
plot_title=strcat([frequency,' waveform'])
plot_title1=strcat([frequency,' Kaiser window waveform α=2'])
fig_axis=[0,1,-2,2]
scf(10);clf;
subplot(2,1,1)
plot(t',S_sin(n,:)')
title(plot_title,'fontsize',4)
mtlb_axis(fig_axis)
legend(deg(1),deg(2),deg(3),deg(4),deg(5),-1)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
subplot(2,1,2)
plot(t',WS_sin(n,:)')
title(plot_title1,'fontsize',4)
mtlb_axis(fig_axis)
legend(deg(1),deg(2),deg(3),deg(4),deg(5),-1)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=20 // ここで表示幅を変える。
ft_axis=[0,fmax,-100,10]
scf(20);clf;
subplot(2,1,1)
plot([0:fmax]',20.*log10(M_sin(n,1:fmax+1)'))
title('Rectangular window separation','fontsize',4)
legend(deg(1),deg(2),deg(3),deg(4),deg(5),-1)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
subplot(2,1,2)
plot([0:fmax]',20.*log10(WM_sin(n,1:fmax+1)'))
title('Kaiser window separation α=2','fontsize',4)
mtlb_axis(ft_axis)
legend(deg(1),deg(2),deg(3),deg(4),deg(5),-1)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
これは 2.5Hz に対しての結果です。
もう一つ 10.5Hz に対しての結果も載せておきます。
注目したいのは、周波数 0 = DC ですね。90deg の時だけが -∞ ということでレベルがゼロになります。
これは 2.5Hz も 10.5Hz も位相が 90deg 遅れで入ってくると、面積がゼロになるからです。見やすい方 2.5Hz で波形を確認すると赤い線は 90deg の時の波形ですが、信号レベル 1 から始まって 0.5sec のところで 0 になり、1sec（実際はサンプル数が偶数なのでない）のところでは -1 になります。このため DC レベルはただ積分しただけだと 0 ということになります。
ただし、2.5Hz も 10.5Hz の時も高い周波数側でのノイズレベルは他の位相に比べて上がっています。痛し痒しというところでしょうか。
それでも 60dB の S / N があるので良しとしましょうか。
もう一つ注目したいのは、いずれの周波数でもカイザー窓の場合は位相に関係なく端数の両隣の整数のレベルが同じになっていることです。
前回カイザー窓では周波数分解能は、3Hz ぐらいだとしましたが、現実にこのような端数の周波数が入ってきた場合、2Hz と 3Hz が両方あると見るのか、2.5Hz だけがあると見るのかちょっと判断がつきません。
これはクリアできるのでしょうか。色々やってみることにします。
カイザー窓を使うと矩形窓よりも端数の周波数に対して S / N が良くなると云うことを確認しました。では、その端数の周波数が観測時間内に色々な位相で入って来たら（普通はそうなる）FFT の結果はどうなるのか見てみることにします。
スクリプトです。だんだん scilab 自体へのレベルが上がってきたのか、今まで使っていなかった関数も結構ありますので、これなんだ？と思った方はヘルプを見て下さい。ここで質問してくれてもいいです。
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
n=1:5
phi(n)=(n-1)*45 // 正弦波の位相: 0deg, 45deg, 90deg, 135deg, 180deg
deg=string(phi)
fs1=2.5 // 信号周波数、後半では10.5Hz に書き換える。
amp1=1 // 信号周波数の振幅
frequency=strcat([string(fs1),'Hz'])
Kaiser=window('kr',N,%pi*2) // Kaiser窓
S_sin=zeros(5,2048);
for i=1:5
S_sin(i,:)=amp1*sin(2*%pi*(t*fs1+phi(i)/360))
WS_sin(i,:)=Kaiser.*S_sin(i,:)
f_sin(i,:)=fft(S_sin(i,:)) // 関数 fft を呼び出す。
M_sin(i,:)=abs(f_sin(i,:))*2/N // 振幅を計算
wf_sin(i,:)=fft(WS_sin(i,:))
WM_sin(i,:)=abs(wf_sin(i,:))*2/N
end
plot_title=strcat([frequency,' waveform'])
plot_title1=strcat([frequency,' Kaiser window waveform α=2'])
fig_axis=[0,1,-2,2]
scf(10);clf;
subplot(2,1,1)
plot(t',S_sin(n,:)')
title(plot_title,'fontsize',4)
mtlb_axis(fig_axis)
legend(deg(1),deg(2),deg(3),deg(4),deg(5),-1)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
subplot(2,1,2)
plot(t',WS_sin(n,:)')
title(plot_title1,'fontsize',4)
mtlb_axis(fig_axis)
legend(deg(1),deg(2),deg(3),deg(4),deg(5),-1)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=20 // ここで表示幅を変える。
ft_axis=[0,fmax,-100,10]
scf(20);clf;
subplot(2,1,1)
plot([0:fmax]',20.*log10(M_sin(n,1:fmax+1)'))
title('Rectangular window separation','fontsize',4)
legend(deg(1),deg(2),deg(3),deg(4),deg(5),-1)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
subplot(2,1,2)
plot([0:fmax]',20.*log10(WM_sin(n,1:fmax+1)'))
title('Kaiser window separation α=2','fontsize',4)
mtlb_axis(ft_axis)
legend(deg(1),deg(2),deg(3),deg(4),deg(5),-1)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
これは 2.5Hz に対しての結果です。
もう一つ 10.5Hz に対しての結果も載せておきます。
注目したいのは、周波数 0 = DC ですね。90deg の時だけが -∞ ということでレベルがゼロになります。
これは 2.5Hz も 10.5Hz も位相が 90deg 遅れで入ってくると、面積がゼロになるからです。見やすい方 2.5Hz で波形を確認すると赤い線は 90deg の時の波形ですが、信号レベル 1 から始まって 0.5sec のところで 0 になり、1sec（実際はサンプル数が偶数なのでない）のところでは -1 になります。このため DC レベルはただ積分しただけだと 0 ということになります。
ただし、2.5Hz も 10.5Hz の時も高い周波数側でのノイズレベルは他の位相に比べて上がっています。痛し痒しというところでしょうか。
それでも 60dB の S / N があるので良しとしましょうか。
もう一つ注目したいのは、いずれの周波数でもカイザー窓の場合は位相に関係なく端数の両隣の整数のレベルが同じになっていることです。
前回カイザー窓では周波数分解能は、3Hz ぐらいだとしましたが、現実にこのような端数の周波数が入ってきた場合、2Hz と 3Hz が両方あると見るのか、2.5Hz だけがあると見るのかちょっと判断がつきません。
これはクリアできるのでしょうか。色々やってみることにします。