LTspice：「まだ退治できていないんだろう」
画蔵：「攻撃力、守備力共に高いので一筋縄ではいかない」
さて、分析したい信号にランプ関数＝観測時間の中でゆっくり変化するような信号が含まれる場合は、スペクトラム分析を掛けると見たい周波数の真値が分からなくなる恐れがあるということを説明してきました。
ちょっと大げさな信号ではありますが、次のような信号が観測時間に入ってきたとします。
この信号の中で本当に見たい信号は 10Hz で、ゆっくり上昇している成分はたとえば温度変化などによって起きる電気オフセットでこれは無視して良い、と思っていることを想定して分析をしたいのですが、このままスペクトラム分析を掛けると次のような信号だと思って結果を出力してしまいます。
これでは「のこぎり波」を分析しているのか、重畳されている信号を分析しているのか分からなくなります。
そこで、前回触れたような HPF をいれることで観測時間に現れるランプ関数成分を少しでも取り除いておきたいわけです。
ここで注意しておきたいのは、HPF といえども除去したい周波数が低い場合は除去し終わるまでに時間が掛かると云うことです。なので信号が現れてから、ちょっと待って観測を開始することになります。
次のスクリプトを順番に見ながら解説しておこうと思います。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
t_range=2 // 測定時間幅 ←
N=2000 // サンプリング周波数
t=0:1/N:t_range-1/N // 時間内のサンプル数
fig_axis=[0,2,-2.0,2.0]
fs=10 // 正弦波の周波数
amp_sin=0.1 // 正弦波の振幅
Ramp(1:N*t_range)=[-(N*t_range)/2+0.5:1:(N*t_range)/2-0.5]/N*2 // ランプ関数
phi=100/N // 正弦波の位相
S_sin(1,1:N*t_range)=sin(2*%pi*(t+phi)*fs)*amp_sin // 本来の正弦波
Signal=Ramp+S_sin
scf(10);clf;plot(t,Signal)
mtlb_axis(fig_axis)
title('Ramp and 10Hz signal','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
count=4
for i=1:count
Ver_signal(1,1+(i-1)*N*t_range:N*i*t_range)=Signal(1,:)
end
scf(11);plot(0:1/N:count*t_range-1/N,Ver_signal)
title('Virtuall analysis signal','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
// HPF を掛ける
s=poly(0,'s')
dt=1/N
dires=[]
cutoffrate=1/N // サンプリング周波数に対する HPF のカットオフ周波数比
fc=1/dt*cutoffrate // カットオフ周波数
HPF=1-fc/(fc+s/2/%pi) // HPF 伝達関数式
HPFsys=syslin('c',HPF);
clf(20);scf(20);clf;bode(HPFsys,1,1e3); // HPF ボーデ線図
title('HPF bode','fontsize',4)
// t=0:dt:dt/cutoffrate
scf(30);clf
HPF_signal=csim(Signal,t,HPFsys) // HPF を通したあとの信号
plot(t,HPF_signal)
title('Signal with HPF','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=50 // DFT の周波数範囲
tr=t(N/2*t_range:N*t_range-1) // DFT を掛ける時間範囲：HPF の応答が落ち着いたところを取り出している。
for f1=1:fmax
x=sin(2*%pi*tr*f1) // 正弦波
y=cos(2*%pi*tr*f1) // 余弦波
S_combo_sin=x.*S_sin(N/2*t_range:N*t_range-1)
S_combo_cos=y.*S_sin(N/2*t_range:N*t_range-1)
S_sin_level(1,f1)=sum(S_combo_sin)./N*2
S_cos_level(1,f1)=sum(S_combo_cos)./N*2
S_level(1,f1)=((S_sin_level(f1))^2+(S_cos_level(f1))^2)^0.5
HPF_combo_sin=x.*HPF_signal(N/2*t_range:N*t_range-1)
HPF_combo_cos=y.*HPF_signal(N/2*t_range:N*t_range-1)
HPF_sin_level(1,f1)=sum(HPF_combo_sin)./N*2
HPF_cos_level(1,f1)=sum(HPF_combo_cos)./N*2
HPF_level(1,f1)=((HPF_sin_level(f1))^2+(HPF_cos_level(f1))^2)^0.5
combo_sin=x.*Signal(N/2*t_range:N*t_range-1)
combo_cos=y.*Signal(N/2*t_range:N*t_range-1)
sin_level(1,f1)=sum(combo_sin)./N*2
cos_level(1,f1)=sum(combo_cos)./N*2
level(1,f1)=((sin_level(f1))^2+(cos_level(f1))^2)^0.5
end
ft_axis=[0,fmax,-1.0,1.0]
scf(40);clf;plot(1:fmax,level,1:fmax,sin_level,1:fmax,cos_level)
mtlb_axis(ft_axis)
ftdb_axis=[0,fmax,-100,0]
scf(50);clf;plot(1:fmax,20*log10(S_level),1:fmax,20*log10(level),1:fmax,20*log10(HPF_level))
title('Frequency signal level','fontsize',4)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
//xtitle('Frequency sugnal level','Frequency(Hz)','Level(dB)')
legend('only 10Hz','With ramp signal','With ramp and HPF')
mtlb_axis(ftdb_axis)
だいぶ長くなってしまいましたが、少しずつやっていきます。
まず、始めの方で t_range という変数を設定しています。これは今までは観測時間 1 秒としていたのですが、その整数倍だけデータを取得することにします。ここでは 2 秒間データを取得するものとします。
ランプ関数と 10Hz の合成信号がグラフィックウィンド 10 に表示されています。2 秒間です。
これをこのままスペクトラム分析に掛けると、実際には先ほどのような波形を分析した結果が出てしまいます。
そこで「// HPF を掛ける」ところから、元の信号に HPF を掛けてみることにします。
伝達関数は、
(s /(2πfc) / (1 + s /(2πfc) でここでは、fc = 1Hz です。
伝達特性は、このようになります。
1Hz のところで、利得が -3dB、位相が 45度ですので目標通りの値になっています。
使い回しのスクリプトなので読みづらいと思いますが、ご容赦下さい。
取得されたデータにこの HPF を掛けた結果がこのようになります。
最初の方はどうしても急峻な立ち上がりを持ってしまいますが、1 秒後ぐらいなら落ち着いています。そこで点線が示している 1 秒後からデータを観測したと云うことにして、そこを分析することにします。直流成分は残ってしまいますが、ここでは影響は出ないので無視できます。
結果はこうなりました。（対数表示）
青線は、10Hz だけの信号だった場合で、当然 10Hz 以外は 0 でこれが本来得たい結果です。
ところが前回やったように、ランプ関数が乗り込んでいると緑線のような結果が出て、これはランプ関数のスペクトルの内 10Hz の成分がどういう位相かで、本来の信号のレベルが振り回されることになります。ここでは少しずらすことで一応ピークが見えるようにしていますが、運が悪いと逆に凹んでしまいます。
赤線は HPF で 1Hz 以下を減衰させたので、影響は減って本来値に近づいています。
ちなみに HPF のカットオフ周波数を 2Hz にした場合は 10Hz 以下の成分はいなくなります。ただし、今度は HPF の影響で 10Hz の信号レベルが下がってしまいますので、観察したい周波数を考えながら HPF の特性を決めて下さい。
たいていの用途では気にする必要はないと思いますが、組み込み機器で測定時間などに制限がある場合は注意が必要です。
scilab：「一応中ボス退治に成功したわけ？」
画蔵：「退治しきれなかったが、悪行をある程度は押さえ込むことには成功したと云うことだ」
画蔵：「攻撃力、守備力共に高いので一筋縄ではいかない」
さて、分析したい信号にランプ関数＝観測時間の中でゆっくり変化するような信号が含まれる場合は、スペクトラム分析を掛けると見たい周波数の真値が分からなくなる恐れがあるということを説明してきました。
ちょっと大げさな信号ではありますが、次のような信号が観測時間に入ってきたとします。
この信号の中で本当に見たい信号は 10Hz で、ゆっくり上昇している成分はたとえば温度変化などによって起きる電気オフセットでこれは無視して良い、と思っていることを想定して分析をしたいのですが、このままスペクトラム分析を掛けると次のような信号だと思って結果を出力してしまいます。
これでは「のこぎり波」を分析しているのか、重畳されている信号を分析しているのか分からなくなります。
そこで、前回触れたような HPF をいれることで観測時間に現れるランプ関数成分を少しでも取り除いておきたいわけです。
ここで注意しておきたいのは、HPF といえども除去したい周波数が低い場合は除去し終わるまでに時間が掛かると云うことです。なので信号が現れてから、ちょっと待って観測を開始することになります。
次のスクリプトを順番に見ながら解説しておこうと思います。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
t_range=2 // 測定時間幅 ←
N=2000 // サンプリング周波数
t=0:1/N:t_range-1/N // 時間内のサンプル数
fig_axis=[0,2,-2.0,2.0]
fs=10 // 正弦波の周波数
amp_sin=0.1 // 正弦波の振幅
Ramp(1:N*t_range)=[-(N*t_range)/2+0.5:1:(N*t_range)/2-0.5]/N*2 // ランプ関数
phi=100/N // 正弦波の位相
S_sin(1,1:N*t_range)=sin(2*%pi*(t+phi)*fs)*amp_sin // 本来の正弦波
Signal=Ramp+S_sin
scf(10);clf;plot(t,Signal)
mtlb_axis(fig_axis)
title('Ramp and 10Hz signal','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
count=4
for i=1:count
Ver_signal(1,1+(i-1)*N*t_range:N*i*t_range)=Signal(1,:)
end
scf(11);plot(0:1/N:count*t_range-1/N,Ver_signal)
title('Virtuall analysis signal','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
// HPF を掛ける
s=poly(0,'s')
dt=1/N
dires=[]
cutoffrate=1/N // サンプリング周波数に対する HPF のカットオフ周波数比
fc=1/dt*cutoffrate // カットオフ周波数
HPF=1-fc/(fc+s/2/%pi) // HPF 伝達関数式
HPFsys=syslin('c',HPF);
clf(20);scf(20);clf;bode(HPFsys,1,1e3); // HPF ボーデ線図
title('HPF bode','fontsize',4)
// t=0:dt:dt/cutoffrate
scf(30);clf
HPF_signal=csim(Signal,t,HPFsys) // HPF を通したあとの信号
plot(t,HPF_signal)
title('Signal with HPF','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
fmax=50 // DFT の周波数範囲
tr=t(N/2*t_range:N*t_range-1) // DFT を掛ける時間範囲：HPF の応答が落ち着いたところを取り出している。
for f1=1:fmax
x=sin(2*%pi*tr*f1) // 正弦波
y=cos(2*%pi*tr*f1) // 余弦波
S_combo_sin=x.*S_sin(N/2*t_range:N*t_range-1)
S_combo_cos=y.*S_sin(N/2*t_range:N*t_range-1)
S_sin_level(1,f1)=sum(S_combo_sin)./N*2
S_cos_level(1,f1)=sum(S_combo_cos)./N*2
S_level(1,f1)=((S_sin_level(f1))^2+(S_cos_level(f1))^2)^0.5
HPF_combo_sin=x.*HPF_signal(N/2*t_range:N*t_range-1)
HPF_combo_cos=y.*HPF_signal(N/2*t_range:N*t_range-1)
HPF_sin_level(1,f1)=sum(HPF_combo_sin)./N*2
HPF_cos_level(1,f1)=sum(HPF_combo_cos)./N*2
HPF_level(1,f1)=((HPF_sin_level(f1))^2+(HPF_cos_level(f1))^2)^0.5
combo_sin=x.*Signal(N/2*t_range:N*t_range-1)
combo_cos=y.*Signal(N/2*t_range:N*t_range-1)
sin_level(1,f1)=sum(combo_sin)./N*2
cos_level(1,f1)=sum(combo_cos)./N*2
level(1,f1)=((sin_level(f1))^2+(cos_level(f1))^2)^0.5
end
ft_axis=[0,fmax,-1.0,1.0]
scf(40);clf;plot(1:fmax,level,1:fmax,sin_level,1:fmax,cos_level)
mtlb_axis(ft_axis)
ftdb_axis=[0,fmax,-100,0]
scf(50);clf;plot(1:fmax,20*log10(S_level),1:fmax,20*log10(level),1:fmax,20*log10(HPF_level))
title('Frequency signal level','fontsize',4)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
//xtitle('Frequency sugnal level','Frequency(Hz)','Level(dB)')
legend('only 10Hz','With ramp signal','With ramp and HPF')
mtlb_axis(ftdb_axis)
だいぶ長くなってしまいましたが、少しずつやっていきます。
まず、始めの方で t_range という変数を設定しています。これは今までは観測時間 1 秒としていたのですが、その整数倍だけデータを取得することにします。ここでは 2 秒間データを取得するものとします。
ランプ関数と 10Hz の合成信号がグラフィックウィンド 10 に表示されています。2 秒間です。
これをこのままスペクトラム分析に掛けると、実際には先ほどのような波形を分析した結果が出てしまいます。
そこで「// HPF を掛ける」ところから、元の信号に HPF を掛けてみることにします。
伝達関数は、
(s /(2πfc) / (1 + s /(2πfc) でここでは、fc = 1Hz です。
伝達特性は、このようになります。
1Hz のところで、利得が -3dB、位相が 45度ですので目標通りの値になっています。
使い回しのスクリプトなので読みづらいと思いますが、ご容赦下さい。
取得されたデータにこの HPF を掛けた結果がこのようになります。
最初の方はどうしても急峻な立ち上がりを持ってしまいますが、1 秒後ぐらいなら落ち着いています。そこで点線が示している 1 秒後からデータを観測したと云うことにして、そこを分析することにします。直流成分は残ってしまいますが、ここでは影響は出ないので無視できます。
結果はこうなりました。（対数表示）
青線は、10Hz だけの信号だった場合で、当然 10Hz 以外は 0 でこれが本来得たい結果です。
ところが前回やったように、ランプ関数が乗り込んでいると緑線のような結果が出て、これはランプ関数のスペクトルの内 10Hz の成分がどういう位相かで、本来の信号のレベルが振り回されることになります。ここでは少しずらすことで一応ピークが見えるようにしていますが、運が悪いと逆に凹んでしまいます。
赤線は HPF で 1Hz 以下を減衰させたので、影響は減って本来値に近づいています。
ちなみに HPF のカットオフ周波数を 2Hz にした場合は 10Hz 以下の成分はいなくなります。ただし、今度は HPF の影響で 10Hz の信号レベルが下がってしまいますので、観察したい周波数を考えながら HPF の特性を決めて下さい。
たいていの用途では気にする必要はないと思いますが、組み込み機器で測定時間などに制限がある場合は注意が必要です。
scilab：「一応中ボス退治に成功したわけ？」
画蔵：「退治しきれなかったが、悪行をある程度は押さえ込むことには成功したと云うことだ」