今回は、端数でない周波数にランダムノイズが加わった場合について調べてみます。
ランダムノイズですから、周波数もレベルもまちまちということで、要は端数の周波数も含めば端数でないものも含んでいるわけで、前回の結果から窓関数を通しておけば多少 S/N の改善ができそうな気がしますが、どうでしょうか。
スクリプトは以下のようにしています。
ランダムノイズを観測期間すべてに渡って発生させていますが、このままだとナイキスト周波数を越えたものまで入ってきてしまうので、前に使った 200 Hz の一次 LPF を通してエイリアシングの影響を軽減させておきます。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
rand('seed',0) // 乱数のもと
randomdata=((rand(1,2048,'uniform'))-0.5) // 乱数発生
s=poly(0,'s')
cutoffrate=10 // サンプリング周波数に対する LPF のカットオフ周波数比。
fc=N/cutoffrate // カットオフ周波数。ここでは 204.8Hz
LPF=fc/(fc+s/2/%pi) // LPF 伝達関数式
LPFsys=syslin('c',LPF);
scf(1);clf;bode(LPFsys,10,1e4); // LPF ボーデ線図
title('LPF bode','fontsize',4)
Noise_signal=csim(randomdata,t,LPFsys) // LPF を通したあとの信号
fig_axis=[0,1,-2,2]
scf(2);clf;plot(t,Noise_signal)
title('Noise signal','fontsize',4)
mtlb_axis(fig_axis)
fs1=4 // 低い信号周波数
fs2=13 // 観測したい信号の周波数
amp1=1 // 低い信号周波数の振幅
amp2=0.01 // 観測したい信号の振幅
amp3=0.05 // ノイズのレベル
Kaiser=window('kr',N,%pi*2) // Kaiser窓
S_sin=amp1*sin(2*%pi*t*fs1)+amp2*sin(2*%pi*t*fs2)+amp3*Noise_signal
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
scf(3);clf;plot(t,S_sin,t,WS_sin)
title('4Hz + 13Hz + Noise signal α=2','fontsize',4)
mtlb_axis(fig_axis)
legend('Rectangular window','Kaiser window')
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=50 // ここで表示幅を変える。
ft_axis=[0,fmax,-100,10]
scf(4);plot([0:fmax],20.*log10(M_sin(1:fmax+1)),[0:fmax],20.*log10(WM_sin(1:fmax+1)))
title(' 4Hz + 13Hz + Noise Spectrum α=2','fontsize',4)
mtlb_axis(ft_axis)
legend('Rectangular window','Kaiser window')
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
LPF の特性はこうなっています。2 KHz で -20dB ですから、0 Hz に対するエイリアシングは 1 / 10 ぐらいは出てしまいそうです。
印加したノイズ波形です。
合成波形です。4 Hz の振幅は 1、13Hz は 0.01、先ほどのノイズ波形の利得は 0.05 です。
FFT を通した結果です。
だいたい予想通りで、カイザー窓は周波数のピークを周波数方向にならす効果があるようで、本来の信号も周波数分解能が前回と同様広がってしまいますが、ノイズ成分もピークを抑えながら周波数方向にならしてしまうため、本来の信号が見やすくなっている、、、と、思います．．．。
簡単に言うと、いわゆる時間領域の LPF が時間軸方向に信号をならす働きがあるのに対して、カイザー窓は周波数領域でのスペクトラムをならすフィルターというところでしょうか。
カイザー窓はαを変えることで特性を変化させることが出来るので、見たい信号とノイズとの関係や観測時間などを鑑みながら色々試してみると良いと思います。
ランダムノイズですから、周波数もレベルもまちまちということで、要は端数の周波数も含めば端数でないものも含んでいるわけで、前回の結果から窓関数を通しておけば多少 S/N の改善ができそうな気がしますが、どうでしょうか。
スクリプトは以下のようにしています。
ランダムノイズを観測期間すべてに渡って発生させていますが、このままだとナイキスト周波数を越えたものまで入ってきてしまうので、前に使った 200 Hz の一次 LPF を通してエイリアシングの影響を軽減させておきます。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
N=2048
t=0:1/N:1-1/N
rand('seed',0) // 乱数のもと
randomdata=((rand(1,2048,'uniform'))-0.5) // 乱数発生
s=poly(0,'s')
cutoffrate=10 // サンプリング周波数に対する LPF のカットオフ周波数比。
fc=N/cutoffrate // カットオフ周波数。ここでは 204.8Hz
LPF=fc/(fc+s/2/%pi) // LPF 伝達関数式
LPFsys=syslin('c',LPF);
scf(1);clf;bode(LPFsys,10,1e4); // LPF ボーデ線図
title('LPF bode','fontsize',4)
Noise_signal=csim(randomdata,t,LPFsys) // LPF を通したあとの信号
fig_axis=[0,1,-2,2]
scf(2);clf;plot(t,Noise_signal)
title('Noise signal','fontsize',4)
mtlb_axis(fig_axis)
fs1=4 // 低い信号周波数
fs2=13 // 観測したい信号の周波数
amp1=1 // 低い信号周波数の振幅
amp2=0.01 // 観測したい信号の振幅
amp3=0.05 // ノイズのレベル
Kaiser=window('kr',N,%pi*2) // Kaiser窓
S_sin=amp1*sin(2*%pi*t*fs1)+amp2*sin(2*%pi*t*fs2)+amp3*Noise_signal
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
scf(3);clf;plot(t,S_sin,t,WS_sin)
title('4Hz + 13Hz + Noise signal α=2','fontsize',4)
mtlb_axis(fig_axis)
legend('Rectangular window','Kaiser window')
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=50 // ここで表示幅を変える。
ft_axis=[0,fmax,-100,10]
scf(4);plot([0:fmax],20.*log10(M_sin(1:fmax+1)),[0:fmax],20.*log10(WM_sin(1:fmax+1)))
title(' 4Hz + 13Hz + Noise Spectrum α=2','fontsize',4)
mtlb_axis(ft_axis)
legend('Rectangular window','Kaiser window')
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
LPF の特性はこうなっています。2 KHz で -20dB ですから、0 Hz に対するエイリアシングは 1 / 10 ぐらいは出てしまいそうです。
印加したノイズ波形です。
合成波形です。4 Hz の振幅は 1、13Hz は 0.01、先ほどのノイズ波形の利得は 0.05 です。
FFT を通した結果です。
だいたい予想通りで、カイザー窓は周波数のピークを周波数方向にならす効果があるようで、本来の信号も周波数分解能が前回と同様広がってしまいますが、ノイズ成分もピークを抑えながら周波数方向にならしてしまうため、本来の信号が見やすくなっている、、、と、思います．．．。
簡単に言うと、いわゆる時間領域の LPF が時間軸方向に信号をならす働きがあるのに対して、カイザー窓は周波数領域でのスペクトラムをならすフィルターというところでしょうか。
カイザー窓はαを変えることで特性を変化させることが出来るので、見たい信号とノイズとの関係や観測時間などを鑑みながら色々試してみると良いと思います。