LTspice：「順調にレベルが上がっているじゃないか」
画蔵：「ニューイシューを手に入れたから」
EXCEL：「横文字でごまかしているような」
scilab：「・・・（働き過ぎで言葉も出ない）」
前回の解説と計算はどうも怪しげさが満載で、何をやりたかったかというとウィキペディアに載っている画像を再現させて見たかったのです。雰囲気はそんな感じでと思うのですが、元画像の数値が読み取れないので上手くいっているのかはよく分かりません。
ですが、あまりここの再現にこだわるよりも使ったときに何が起こるという方が重要だと思うのでどんどん先に進むことにします。
まずはお試しと云うことで、4Hz と 13Hz が同レベルで混ざっている信号を分析してみます。
スクリプトは以下のようです。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
fs1=4 // 低い信号周波数
fs2=13 // 観測したい信号の周波数
amp1=1 // 低い信号周波数の振幅
amp2=1 // 観測したい信号の振幅
Kaiser=window('kr',N,%pi*3) // Kaiser窓 αは 3
S_sin=amp1*sin(2*%pi*t*fs1)+amp2*sin(2*%pi*t*fs2)
WS_sin=Kaiser.*S_sin
f_sin=fft(S_sin) // 関数 fft を呼び出す。
M_sin=abs(f_sin)*2/N // 振幅を計算
wf_sin=fft(WS_sin)
WM_sin=abs(wf_sin)*2/N
fig_axis=[0,1,-2,2]
scf(10);clf;plot(t,S_sin,t,WS_sin)
title('4Hz + 13Hz signal','fontsize',4)
mtlb_axis(fig_axis)
legend('Rectangular window','Kaiser window')
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=20 // ここで表示幅を変える。
ft_axis=[0,fmax,-100,10]
scf(20);plot([0:fmax],20.*log10(M_sin(1:fmax+1)),[0:fmax],20.*log10(WM_sin(1:fmax+1)))
title(' 4Hz + 13Hz Spectrum','fontsize',4)
mtlb_axis(ft_axis)
legend('Rectangular window','Kaiser window')
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
分析する波形は以下のようになります。青線が元の混ざった波形で、緑線が窓関数を掛けた波形です。時刻 0 の後ろと時刻 1 秒 の手前で振幅がほぼ 0 になっているのが特徴といえます。
FFT した結果は以下のようになります。
今回は観測時間の中で周期性が完結する信号の組み合わせだったので、矩形窓でもきちっと 4 Hz、13Hz が抽出できていてその他の成分は全く見当たりません。それに対してカイザー窓を掛けたものは 4 Hz、13 Hz のレベルが下がっている上にその周辺の周波数にも信号が現れています。さらに非常にレベルが低いとはいえ、19 Hz、20 Hz のところが盛り上がっていますね。その先はどうなっているのでしょうか。ということで、スクリプト中の fmax を 50 にして表示周波数幅を広げてみます。
19 Hz から 29 Hz の間に盛り上がりがあるように見えますね。これがウィキペディアなどで説明されているサイドローブと思って良さそうです。
なんだこれじゃあ、窓関数使う価値ないじゃん！と思いたくなります。
この信号の範囲では実はその通りで、観測時間で完結する信号しかなければ「矩形窓が理論上最も周波数分解能が良い」とウィキペディアに書いてあるとおりで、その結果を確認したと云うことになります。
LTspice：「他の窓関数は？」
画蔵：「特性を見る限り程度の差ということのようだ」
で、少しは窓関数の肩を持たなくてはいけない、というわけで今回はその前振りまでしておきます。
スクリプトの amp2 を 1 から 0.01 にして 13 Hz がいるにはいるが非常に小さくて、でも何とか抽出したいという想定で計算してみます。
結果はこれです。
波形には 13 Hz の痕跡はほとんど見えません。ですが FFT の結果には確かに現れています。
窓関数が掛かっていてもその付近の周波数も付いてきてしまうという事情はそのままですが、見えています。
次回はこれらの周波数が微妙にズレ出すとどうなるかを見てみようと思います。
画蔵：「ニューイシューを手に入れたから」
EXCEL：「横文字でごまかしているような」
scilab：「・・・（働き過ぎで言葉も出ない）」
前回の解説と計算はどうも怪しげさが満載で、何をやりたかったかというとウィキペディアに載っている画像を再現させて見たかったのです。雰囲気はそんな感じでと思うのですが、元画像の数値が読み取れないので上手くいっているのかはよく分かりません。
ですが、あまりここの再現にこだわるよりも使ったときに何が起こるという方が重要だと思うのでどんどん先に進むことにします。
まずはお試しと云うことで、4Hz と 13Hz が同レベルで混ざっている信号を分析してみます。
スクリプトは以下のようです。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
fs1=4 // 低い信号周波数
fs2=13 // 観測したい信号の周波数
amp1=1 // 低い信号周波数の振幅
amp2=1 // 観測したい信号の振幅
Kaiser=window('kr',N,%pi*3) // Kaiser窓 αは 3
S_sin=amp1*sin(2*%pi*t*fs1)+amp2*sin(2*%pi*t*fs2)
WS_sin=Kaiser.*S_sin
f_sin=fft(S_sin) // 関数 fft を呼び出す。
M_sin=abs(f_sin)*2/N // 振幅を計算
wf_sin=fft(WS_sin)
WM_sin=abs(wf_sin)*2/N
fig_axis=[0,1,-2,2]
scf(10);clf;plot(t,S_sin,t,WS_sin)
title('4Hz + 13Hz signal','fontsize',4)
mtlb_axis(fig_axis)
legend('Rectangular window','Kaiser window')
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=20 // ここで表示幅を変える。
ft_axis=[0,fmax,-100,10]
scf(20);plot([0:fmax],20.*log10(M_sin(1:fmax+1)),[0:fmax],20.*log10(WM_sin(1:fmax+1)))
title(' 4Hz + 13Hz Spectrum','fontsize',4)
mtlb_axis(ft_axis)
legend('Rectangular window','Kaiser window')
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
分析する波形は以下のようになります。青線が元の混ざった波形で、緑線が窓関数を掛けた波形です。時刻 0 の後ろと時刻 1 秒 の手前で振幅がほぼ 0 になっているのが特徴といえます。
FFT した結果は以下のようになります。
今回は観測時間の中で周期性が完結する信号の組み合わせだったので、矩形窓でもきちっと 4 Hz、13Hz が抽出できていてその他の成分は全く見当たりません。それに対してカイザー窓を掛けたものは 4 Hz、13 Hz のレベルが下がっている上にその周辺の周波数にも信号が現れています。さらに非常にレベルが低いとはいえ、19 Hz、20 Hz のところが盛り上がっていますね。その先はどうなっているのでしょうか。ということで、スクリプト中の fmax を 50 にして表示周波数幅を広げてみます。
19 Hz から 29 Hz の間に盛り上がりがあるように見えますね。これがウィキペディアなどで説明されているサイドローブと思って良さそうです。
なんだこれじゃあ、窓関数使う価値ないじゃん！と思いたくなります。
この信号の範囲では実はその通りで、観測時間で完結する信号しかなければ「矩形窓が理論上最も周波数分解能が良い」とウィキペディアに書いてあるとおりで、その結果を確認したと云うことになります。
LTspice：「他の窓関数は？」
画蔵：「特性を見る限り程度の差ということのようだ」
で、少しは窓関数の肩を持たなくてはいけない、というわけで今回はその前振りまでしておきます。
スクリプトの amp2 を 1 から 0.01 にして 13 Hz がいるにはいるが非常に小さくて、でも何とか抽出したいという想定で計算してみます。
結果はこれです。
波形には 13 Hz の痕跡はほとんど見えません。ですが FFT の結果には確かに現れています。
窓関数が掛かっていてもその付近の周波数も付いてきてしまうという事情はそのままですが、見えています。
次回はこれらの周波数が微妙にズレ出すとどうなるかを見てみようと思います。