前回に続き、端数の周波数を検出する方法を考えてみます。
なぜ端数になってしまうかというと、観測時間の逆数の周波数の整数倍でしかきれいに分離できないということですが、では観測時間を変えてしまえば良い、というのが今回の発想の元です。
すでに 1 秒 = 2,048 ポイントのデータがあるとして、これをすべて使って FFT を掛ければ確かに 1 Hz 単位ですが、少し減らして見たい周波数付近で観測時間の逆数の周波数の整数倍になる様にしたらどうでしょうか。
FFT 自体はサンプリング数に 2 べき乗でないといけないという制約はない（遅くなるだけ、除く Excel）ので、たとえば 10.3Hz ぐらいを見たいのならば 10.3Hz がちょうど観測時間に整数回、収まる様にしてみるということです。
もっとも組み込み FFT はプログラムの都合上、2 のべき乗でしか扱えない可能性もあるかもしれませんから、だいたい辺りを付けてこのテーマの最初の方にやったような正弦波と余弦波の組み合わせでやることになるでしょう。周波数が絞り込まれていれば計算量はたいしたことありません。
ここでは、scilab が 2 のべき乗以外のサンプル数も扱えるので、その機能を使ってやることにします。
スクリプトはこうです。サンプル数を 2048 から 20 サンプルずつ減らして、端数の周波数のレベルがもっとも高くなるところを探すわけです。今回は 10Hz と 11Hz の間にあると云うことで、観測期間の中にだいたい 10 サイクル入るように数字を選ぶことになります。
10Hz なら 2048 サンプルで 10 サイクルで、11Hz なら 1861 サンプルでほぼ 10 サイクルになるだろうというわけです。
clear // Data をクリア
N0=2048
t=0:1/N0:1-1/N0
fs=10.3
signal=sin(2*%pi*(t*fs))
j=1
for f_step=0:20:180
N=N0-f_step
t=0:1/N0:1-1/N0
alpha=2
n=1
Kaiser=window('kr',N,%pi*alpha) // Kaiser窓
WS_sin(j,1:N)=Kaiser.*signal(1:N)
f_sin(j,1:N)=fft(signal(1:N)) // 関数 fft を呼び出す。
M_sin(j,1:N)=20.*log10(abs(f_sin(j,1:N))*2 ./N) // 振幅を計算
fn(j,1:N)=[0:N-1]*N0/N
disp(max(M_sin(j,:)))
wf_sin(j,1:N)=fft(WS_sin(j,1:N))
WM_sin(j,1:N)=20.*log10(abs(wf_sin(j,1:N))*2 ./N)
disp(max(WM_sin(j,:)))
j=j+1
// end
end
fig_axis=[0,1,-2,2]
scf(10);clf;plot(t,signal)
title('Non integer frequency signal','fontsize',4)
mtlb_axis(fig_axis)
//legend(frequency)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=20
ft_axis=[6,15,-60,10]
scf(20);plot(fn(1:j-1,1:fmax+1)',M_sin(1:j-1,1:fmax+1)')
title('Non integer frequency Rectangular window spectrum','fontsize',4)
mtlb_axis(ft_axis)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
fig_axis=[0,1,-2,2]
scf(30);clf;plot(t',WS_sin')
title('Non integer frequency Kaiser window α=2 waveform','fontsize',4)
mtlb_axis(fig_axis)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
ft_axis=[6,15,-60,10]
scf(40);plot(fn(1:j-1,1:fmax+1)',WM_sin(1:j-1,1:fmax+1)')
title('Non integer frequency Kaiser window α=2 spectrum','fontsize',4)
mtlb_axis(ft_axis)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
矩形窓ではこうなりました。
10.3Hz 付近のところでもっとも振幅が出ているので、検出できていると云っていいと思います。
カイザー窓だとこうなります。
分解能が下がってしまうのでわかりにくいですね。
ここまで精度の良い周波数検出が必要なのかどうかは応用次第ですが、収集したデータを分析して「あとちょっと精度が欲しい」と思ったときに試してみる価値があるといいな、と思います。特に周波数が低い領域では、比率として分解能が落ちてしまうのでこういった方法も価値があるでしょう。
なぜ端数になってしまうかというと、観測時間の逆数の周波数の整数倍でしかきれいに分離できないということですが、では観測時間を変えてしまえば良い、というのが今回の発想の元です。
すでに 1 秒 = 2,048 ポイントのデータがあるとして、これをすべて使って FFT を掛ければ確かに 1 Hz 単位ですが、少し減らして見たい周波数付近で観測時間の逆数の周波数の整数倍になる様にしたらどうでしょうか。
FFT 自体はサンプリング数に 2 べき乗でないといけないという制約はない（遅くなるだけ、除く Excel）ので、たとえば 10.3Hz ぐらいを見たいのならば 10.3Hz がちょうど観測時間に整数回、収まる様にしてみるということです。
もっとも組み込み FFT はプログラムの都合上、2 のべき乗でしか扱えない可能性もあるかもしれませんから、だいたい辺りを付けてこのテーマの最初の方にやったような正弦波と余弦波の組み合わせでやることになるでしょう。周波数が絞り込まれていれば計算量はたいしたことありません。
ここでは、scilab が 2 のべき乗以外のサンプル数も扱えるので、その機能を使ってやることにします。
スクリプトはこうです。サンプル数を 2048 から 20 サンプルずつ減らして、端数の周波数のレベルがもっとも高くなるところを探すわけです。今回は 10Hz と 11Hz の間にあると云うことで、観測期間の中にだいたい 10 サイクル入るように数字を選ぶことになります。
10Hz なら 2048 サンプルで 10 サイクルで、11Hz なら 1861 サンプルでほぼ 10 サイクルになるだろうというわけです。
clear // Data をクリア
N0=2048
t=0:1/N0:1-1/N0
fs=10.3
signal=sin(2*%pi*(t*fs))
j=1
for f_step=0:20:180
N=N0-f_step
t=0:1/N0:1-1/N0
alpha=2
n=1
Kaiser=window('kr',N,%pi*alpha) // Kaiser窓
WS_sin(j,1:N)=Kaiser.*signal(1:N)
f_sin(j,1:N)=fft(signal(1:N)) // 関数 fft を呼び出す。
M_sin(j,1:N)=20.*log10(abs(f_sin(j,1:N))*2 ./N) // 振幅を計算
fn(j,1:N)=[0:N-1]*N0/N
disp(max(M_sin(j,:)))
wf_sin(j,1:N)=fft(WS_sin(j,1:N))
WM_sin(j,1:N)=20.*log10(abs(wf_sin(j,1:N))*2 ./N)
disp(max(WM_sin(j,:)))
j=j+1
// end
end
fig_axis=[0,1,-2,2]
scf(10);clf;plot(t,signal)
title('Non integer frequency signal','fontsize',4)
mtlb_axis(fig_axis)
//legend(frequency)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=20
ft_axis=[6,15,-60,10]
scf(20);plot(fn(1:j-1,1:fmax+1)',M_sin(1:j-1,1:fmax+1)')
title('Non integer frequency Rectangular window spectrum','fontsize',4)
mtlb_axis(ft_axis)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
fig_axis=[0,1,-2,2]
scf(30);clf;plot(t',WS_sin')
title('Non integer frequency Kaiser window α=2 waveform','fontsize',4)
mtlb_axis(fig_axis)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
ft_axis=[6,15,-60,10]
scf(40);plot(fn(1:j-1,1:fmax+1)',WM_sin(1:j-1,1:fmax+1)')
title('Non integer frequency Kaiser window α=2 spectrum','fontsize',4)
mtlb_axis(ft_axis)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
矩形窓ではこうなりました。
10.3Hz 付近のところでもっとも振幅が出ているので、検出できていると云っていいと思います。
カイザー窓だとこうなります。
分解能が下がってしまうのでわかりにくいですね。
ここまで精度の良い周波数検出が必要なのかどうかは応用次第ですが、収集したデータを分析して「あとちょっと精度が欲しい」と思ったときに試してみる価値があるといいな、と思います。特に周波数が低い領域では、比率として分解能が落ちてしまうのでこういった方法も価値があるでしょう。