※ スクリプトに抜けがありました。FFTの結果表示を加えました。（2016/01/17）
窓関数として観測時間が十分あって、サンプリング数も十分あれば何も考えずに矩形窓を使えば良いですし、ほとんどの応用ではそれを意識せずに使っていると思います。しかし、組み込み機器などでやすでに採られている限られた時間幅のデータを分析しようとした場合、端数の周波数成分が混ざっているとノイズフロアが上がってしまい、観測したい信号レベルが小さいと見えなくなってしまうことがあるので、窓関数としてカイザー窓を使って分解能を犠牲にしても S / N を良くする方法がある、ということを解説して来ました。
一方ランプ関数などがあった場合、その影響などを取り除くために HPF などを入れて低い周波数成分をカットして 1 秒間の観測時間で精度測定良く出来るのは 10Hz 以上だろう、ということも確認してみました。ですが、そういうものはないから出来るだけ低い周波数まで観測したいということもあるかと思います。そういう場合は DC 成分などは 0 Hz として検出されるだけだから無視しても良さそう、ということにしておきましたが、今回はその DC 成分がカイザー窓を使った場合どのような影響を検出結果に与えるか調べてみます。
スクリプトです。前回のものを加工したので、整理されていません。ご容赦。
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
n=1
phi(n)=(n-1)*45 // 正弦波の位相: 0deg, 45deg, 90deg, 135deg, 180deg
deg=string(phi)
fs1=10 // 信号周波数
amp1=1 // 信号周波数の振幅
dc=1 // 直流成分
frequency=strcat([string(fs1),'Hz'])
Kaiser=window('kr',N,%pi*2) // Kaiser窓
S_sin=zeros(1,N);
i=1
S_sin(i,:)=amp1*sin(2*%pi*(t*fs1+phi(i)/360))+dc
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
plot_title=strcat([frequency,' waveform'])
plot_title1=strcat([frequency,' Kaiser window waveform α=2'])
fig_axis=[0,1,-1,3]
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','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 α=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)
波形はこうなりました。
FFT の結果はこうです。
なんとなく予想されたことですが、DC 成分を含んだ信号をカイザー窓を通してから FFT をかけると、1 Hz、2 Hz ぐらいの信号がいるように見えてしまいます。矩形窓の場合は 0 Hz がいるだけです。
波形の方を見てみると、矩形窓は全体にオフセットが乗っているだけなのに対し、カイザー窓の方は時間内で何か周期性があるかのような波形になってしまいます。それでこのような結果になるわけです。
FFT というのは、観測時間内を評価する際それらがそのまま繰り返されるような信号だと思って計算する、と説明しました。つまり今回の DC 成分を含んだ信号をカイザー窓を通して評価すると云うことは次のような波形を評価していることになります。
なるほどこれでは低い周波数が重畳したような波形になってしまいますので、上記のような結果が出てしまうのもやむなし、ということです。
低い周波数も見たいからランプ関数対策はしたくない、だけどこういうのも困る、という用途もあると思いますので、カイザー窓を使う場合は、実は先に DC カット処理（観測時間内の平均を計算して、あらかじめ全部のデータから差し引いておく）をしておく必要があるわけです。
LTspice：「雑魚と思っていた DC 成分がカイザー窓を使った瞬間、迷惑な存在になった、ということか？」
画蔵：「そう。うっかりしがちだが、ちゃんと対策しておかないと面倒なことになる」
窓関数として観測時間が十分あって、サンプリング数も十分あれば何も考えずに矩形窓を使えば良いですし、ほとんどの応用ではそれを意識せずに使っていると思います。しかし、組み込み機器などでやすでに採られている限られた時間幅のデータを分析しようとした場合、端数の周波数成分が混ざっているとノイズフロアが上がってしまい、観測したい信号レベルが小さいと見えなくなってしまうことがあるので、窓関数としてカイザー窓を使って分解能を犠牲にしても S / N を良くする方法がある、ということを解説して来ました。
一方ランプ関数などがあった場合、その影響などを取り除くために HPF などを入れて低い周波数成分をカットして 1 秒間の観測時間で精度測定良く出来るのは 10Hz 以上だろう、ということも確認してみました。ですが、そういうものはないから出来るだけ低い周波数まで観測したいということもあるかと思います。そういう場合は DC 成分などは 0 Hz として検出されるだけだから無視しても良さそう、ということにしておきましたが、今回はその DC 成分がカイザー窓を使った場合どのような影響を検出結果に与えるか調べてみます。
スクリプトです。前回のものを加工したので、整理されていません。ご容赦。
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
n=1
phi(n)=(n-1)*45 // 正弦波の位相: 0deg, 45deg, 90deg, 135deg, 180deg
deg=string(phi)
fs1=10 // 信号周波数
amp1=1 // 信号周波数の振幅
dc=1 // 直流成分
frequency=strcat([string(fs1),'Hz'])
Kaiser=window('kr',N,%pi*2) // Kaiser窓
S_sin=zeros(1,N);
i=1
S_sin(i,:)=amp1*sin(2*%pi*(t*fs1+phi(i)/360))+dc
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
plot_title=strcat([frequency,' waveform'])
plot_title1=strcat([frequency,' Kaiser window waveform α=2'])
fig_axis=[0,1,-1,3]
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','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 α=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)
波形はこうなりました。
FFT の結果はこうです。
なんとなく予想されたことですが、DC 成分を含んだ信号をカイザー窓を通してから FFT をかけると、1 Hz、2 Hz ぐらいの信号がいるように見えてしまいます。矩形窓の場合は 0 Hz がいるだけです。
波形の方を見てみると、矩形窓は全体にオフセットが乗っているだけなのに対し、カイザー窓の方は時間内で何か周期性があるかのような波形になってしまいます。それでこのような結果になるわけです。
FFT というのは、観測時間内を評価する際それらがそのまま繰り返されるような信号だと思って計算する、と説明しました。つまり今回の DC 成分を含んだ信号をカイザー窓を通して評価すると云うことは次のような波形を評価していることになります。
なるほどこれでは低い周波数が重畳したような波形になってしまいますので、上記のような結果が出てしまうのもやむなし、ということです。
低い周波数も見たいからランプ関数対策はしたくない、だけどこういうのも困る、という用途もあると思いますので、カイザー窓を使う場合は、実は先に DC カット処理（観測時間内の平均を計算して、あらかじめ全部のデータから差し引いておく）をしておく必要があるわけです。
LTspice：「雑魚と思っていた DC 成分がカイザー窓を使った瞬間、迷惑な存在になった、ということか？」
画蔵：「そう。うっかりしがちだが、ちゃんと対策しておかないと面倒なことになる」