FFT を行う際どんな窓関数を使うかは、ダイナミックレンジ（S / N）と周波数分解能のトレードオフを考慮しながら決める、と前に書きましたが、今回は周波数分解能を確認してみます、
いきなりスクリプトです。10Hz と 11Hz の組み合わせ、20Hz と 22Hz の組み合わせ、30Hz と 33Hz の組み合わせをそれぞれ確認してみます。振幅は全部 1 です。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
phi=0/N // 正弦波の位相：ここを変化させると分離の様子が変わります。
fs11=10 // 10Hz
fs12=11 // 11Hz
fs21=20 // 20Hz
fs22=22 // 22Hz
fs31=30 // 30Hz
fs32=33 // 33Hz
amp11=1 //
amp12=1 //
amp21=1 //
amp22=1 //
amp31=1 //
amp32=1 //
Kaiser=window('kr',N,%pi*2) // Kaiser窓
S_sin1=amp11*sin(2*%pi*t*fs11)+amp12*sin(2*%pi*(t+phi)*fs12)
S_sin2=amp21*sin(2*%pi*t*fs21)+amp22*sin(2*%pi*(t+phi)*fs22)
S_sin3=amp31*sin(2*%pi*t*fs31)+amp32*sin(2*%pi*(t+phi)*fs32)
WS_sin1=Kaiser.*S_sin1
WS_sin2=Kaiser.*S_sin2
WS_sin3=Kaiser.*S_sin3
f_sin1=fft(S_sin1) // 関数 fft を呼び出す。
M_sin1=abs(f_sin1)*2/N // 振幅を計算
wf_sin1=fft(WS_sin1)
WM_sin1=abs(wf_sin1)*2/N
f_sin2=fft(S_sin2) // 関数 fft を呼び出す。
M_sin2=abs(f_sin2)*2/N // 振幅を計算
wf_sin2=fft(WS_sin2)
WM_sin2=abs(wf_sin2)*2/N
f_sin3=fft(S_sin3) // 関数 fft を呼び出す。
M_sin3=abs(f_sin3)*2/N // 振幅を計算
wf_sin3=fft(WS_sin3)
WM_sin3=abs(wf_sin3)*2/N
fig_axis=[0,1,-2,2]
scf(10);clf;
plot(t,WS_sin1,t,WS_sin2,t,WS_sin3)
title('Kaiser window waveform α=2','fontsize',4)
mtlb_axis(fig_axis)
legend('10Hz + 11Hz','20Hz + 22Hz','30Hz + 33Hz')
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=50 // ここで表示幅を変える。
ft_axis=[0,fmax,-100,10]
scf(20);clf;
subplot(2,1,1)
plot([0:fmax],20.*log10(M_sin1(1:fmax+1)),[0:fmax],20.*log10(M_sin2(1:fmax+1)),[0:fmax],20.*log10(M_sin3(1:fmax+1)))
title('Rectangular window separation','fontsize',4)
legend('10Hz + 11Hz','20Hz + 22Hz','30Hz + 33Hz')
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
subplot(2,1,2)
plot([0:fmax],20.*log10(WM_sin1(1:fmax+1)),[0:fmax],20.*log10(WM_sin2(1:fmax+1)),[0:fmax],20.*log10(WM_sin3(1:fmax+1)))
title('Kaiser window separation α=2','fontsize',4)
mtlb_axis(ft_axis)
legend('10Hz + 11Hz','20Hz + 22Hz','30Hz + 33Hz')
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
カイザー窓を掛けた場合の波形はこうなりました。10Hz などはビートが出ているように見えますね。スクリプト中位相（phi）をいじるとこの様子も変わります。カイザー窓の場合は FFT の結果も変わります。
FFT の結果です。矩形窓とカイザー窓を並べてみました。
さてどうでしょうか。矩形窓の 10Hz と 11Hz は普通のプロットにしたのでくっついて見えますが、要はいずれの周波数もちゃんとレベルが検知できている、ということです。
カイザー窓の方も 10Hz と 11Hz はいるにはいますが、問題は 9Hz や 12Hz も同じぐらいのレベルがいるということですね。さらに 13Hz ぐらいを混ぜるとたくさんの周波数が並び出すような気がします。
で、20Hz と 22Hz も 21Hz まで同じレベルにいるように見えてしまいます。3Hz 差ある 30Hz と 33Hz の組み合わせなら一応分離できてそうだ、というところでしょうか。
Wikipedia などでの解説の通り、矩形窓の方が分解能という意味では優れているといえそうです。
もっとも端数の周波数が現れると、矩形窓は前に出したように端数の前後の周波数のレベルが上がってしまいますから、結局のところ観測時間を少しでも長く取って、サンプル数を増やすのが無理がないということになってしまいます。
普通のスペクトラムアナライザで分析する場合は、あまり気にする必要がないのですが（たくさんデータを取り込むので）、組み込み機器のような場合は時間もサンプル数も制約がありますので、利用する際には気にした方が良いと云うことです。
いきなりスクリプトです。10Hz と 11Hz の組み合わせ、20Hz と 22Hz の組み合わせ、30Hz と 33Hz の組み合わせをそれぞれ確認してみます。振幅は全部 1 です。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
phi=0/N // 正弦波の位相：ここを変化させると分離の様子が変わります。
fs11=10 // 10Hz
fs12=11 // 11Hz
fs21=20 // 20Hz
fs22=22 // 22Hz
fs31=30 // 30Hz
fs32=33 // 33Hz
amp11=1 //
amp12=1 //
amp21=1 //
amp22=1 //
amp31=1 //
amp32=1 //
Kaiser=window('kr',N,%pi*2) // Kaiser窓
S_sin1=amp11*sin(2*%pi*t*fs11)+amp12*sin(2*%pi*(t+phi)*fs12)
S_sin2=amp21*sin(2*%pi*t*fs21)+amp22*sin(2*%pi*(t+phi)*fs22)
S_sin3=amp31*sin(2*%pi*t*fs31)+amp32*sin(2*%pi*(t+phi)*fs32)
WS_sin1=Kaiser.*S_sin1
WS_sin2=Kaiser.*S_sin2
WS_sin3=Kaiser.*S_sin3
f_sin1=fft(S_sin1) // 関数 fft を呼び出す。
M_sin1=abs(f_sin1)*2/N // 振幅を計算
wf_sin1=fft(WS_sin1)
WM_sin1=abs(wf_sin1)*2/N
f_sin2=fft(S_sin2) // 関数 fft を呼び出す。
M_sin2=abs(f_sin2)*2/N // 振幅を計算
wf_sin2=fft(WS_sin2)
WM_sin2=abs(wf_sin2)*2/N
f_sin3=fft(S_sin3) // 関数 fft を呼び出す。
M_sin3=abs(f_sin3)*2/N // 振幅を計算
wf_sin3=fft(WS_sin3)
WM_sin3=abs(wf_sin3)*2/N
fig_axis=[0,1,-2,2]
scf(10);clf;
plot(t,WS_sin1,t,WS_sin2,t,WS_sin3)
title('Kaiser window waveform α=2','fontsize',4)
mtlb_axis(fig_axis)
legend('10Hz + 11Hz','20Hz + 22Hz','30Hz + 33Hz')
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=50 // ここで表示幅を変える。
ft_axis=[0,fmax,-100,10]
scf(20);clf;
subplot(2,1,1)
plot([0:fmax],20.*log10(M_sin1(1:fmax+1)),[0:fmax],20.*log10(M_sin2(1:fmax+1)),[0:fmax],20.*log10(M_sin3(1:fmax+1)))
title('Rectangular window separation','fontsize',4)
legend('10Hz + 11Hz','20Hz + 22Hz','30Hz + 33Hz')
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
subplot(2,1,2)
plot([0:fmax],20.*log10(WM_sin1(1:fmax+1)),[0:fmax],20.*log10(WM_sin2(1:fmax+1)),[0:fmax],20.*log10(WM_sin3(1:fmax+1)))
title('Kaiser window separation α=2','fontsize',4)
mtlb_axis(ft_axis)
legend('10Hz + 11Hz','20Hz + 22Hz','30Hz + 33Hz')
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
カイザー窓を掛けた場合の波形はこうなりました。10Hz などはビートが出ているように見えますね。スクリプト中位相（phi）をいじるとこの様子も変わります。カイザー窓の場合は FFT の結果も変わります。
FFT の結果です。矩形窓とカイザー窓を並べてみました。
さてどうでしょうか。矩形窓の 10Hz と 11Hz は普通のプロットにしたのでくっついて見えますが、要はいずれの周波数もちゃんとレベルが検知できている、ということです。
カイザー窓の方も 10Hz と 11Hz はいるにはいますが、問題は 9Hz や 12Hz も同じぐらいのレベルがいるということですね。さらに 13Hz ぐらいを混ぜるとたくさんの周波数が並び出すような気がします。
で、20Hz と 22Hz も 21Hz まで同じレベルにいるように見えてしまいます。3Hz 差ある 30Hz と 33Hz の組み合わせなら一応分離できてそうだ、というところでしょうか。
Wikipedia などでの解説の通り、矩形窓の方が分解能という意味では優れているといえそうです。
もっとも端数の周波数が現れると、矩形窓は前に出したように端数の前後の周波数のレベルが上がってしまいますから、結局のところ観測時間を少しでも長く取って、サンプル数を増やすのが無理がないということになってしまいます。
普通のスペクトラムアナライザで分析する場合は、あまり気にする必要がないのですが（たくさんデータを取り込むので）、組み込み機器のような場合は時間もサンプル数も制約がありますので、利用する際には気にした方が良いと云うことです。