インパルス応答をサンプリング時刻毎に数列化して、FIR フィルタの係数を決めるという方法で一次 HPF を考えてみましたが、そもそも HPF のインパルス応答を求めること自体が無理な話のようです。そこで今度は一次 LPF を用いて HPF をつくるというアプローチを試してみます。
一次 HPF は源信号から LPF の出力を引き算することで得られます。
s / (1 + Ts) = 1 - 1 / (1 + Ts)
このカラクリの応用例はこちらにもあります。
ピーク検出をしたいけど、微分器は難しい。
http://blogs.yahoo.co.jp/susanoo2001_hero/7543384.html
スクリプトはこれです。
s=poly(0,'s')
z=poly(0,'z')
dt=1d-4
dires=[]
cutoffrate=1/100
fc=1/dt*cutoffrate
wc=2*%pi*fc
LPF=fc/(fc+s/2/%pi)
LPFsys=syslin('c',LPF);
clf(8);scf(8);clf;bode(LPFsys,0.01,1e4);
clf(9);scf(9);clf;bode(1-LPFsys,0.01,1e4);
t=0:dt:dt/cutoffrate
scf(10);clf
plot(t,csim('impuls',t,LPFsys))
dires=[csim('impuls',t,LPFsys)]
DLPF=0
param=dires'
n=size(param,'r')
//param(1)=10000-sum(param)
for i=0:n-1
DLPF=DLPF+z^(-i)*param(i+1)*dt
end
DLPFsys=syslin(dt,DLPF);
scf(11);clf;bode([1-DLPFsys],1e-2,4990)
scf(12);clf;bode([DLPFsys],1e-2,4990)
結果は以下のようになりました。
連続領域の LPF ボーデ線図 連続領域の HPF ボーデ線図
連続領域の LPF インパルス応答
離散領域の LPF ボーデ線図
離散領域の HPF ボーデ線図
FIR による HPF の応答がやっぱり変ですね。1Hz 以下では利得がフラットになってしまい、位相進みもありません。1 - LPF が上手く計算できていないのでしょうか。
と思って LPF の特性をよく見てみると低い周波数での利得が 0dB よりやや低くなっています。完全な 0dB なら 1 - LPF の結果は 0 レベル = - 無限大 のなるはずですので、どうやらそもそも LPF の応答の精度が悪いようです。これもインパルス応答の誤差かも知れない、ということで応答波形をよく見てみると、時刻 0 の出力がやっぱりゼロ。微分方程式から解くと本来は何か数値があるはずです。
実は LPF は波形こそなまらせますが、入力に与えられたエネルギーは後ろにばらまかれるだけで維持されます。従って係数として採用されたインパルス応答数列の総和は入力インパルスのエネルギーと等しくなくてはいけません。ではインパルスのエネルギーとはなんぞや、というとインパルスは時間幅ゼロ、面積 1 ということですから、ここではサンプリング時間が 0.0001 秒ということで、10,000 のエネルギーと言えそうです。
そこで上述のスクリプトの //param(1)=10000-sum(param) を有効にして（ // を外す）計算し直してみます。ここでは param(1) は 328.9626 になっていました。
離散領域の LPF ボーデ線図
離散領域の HPF ボーデ線図
LPF の高域が変な感じになりましたが、低域は 0dB になったようです。ですので HPF ももっともらしくなりました。
これも scilab のインパルス応答の計算誤差を補正することで解決はしました。
scilab：だ～か～ら
では最初からラプラス逆変換を使って伝達関数からいきなりインパルス応答を求めてみたらどうなるでしょうか。
スクリプトはこれです。
s=poly(0,'s')
z=poly(0,'z')
dt=1d-4
dires=[]
cutoffrate=1/100
fc=1/dt*cutoffrate
wc=2*%pi*fc
LPF=fc/(fc+s/2/%pi)
LPFsys=syslin('c',LPF);
clf(8);scf(8);clf;bode(LPFsys,0.01,1e4);
clf(9);scf(9);clf;bode(1-LPFsys,0.01,1e4);
t=0:dt:dt/cutoffrate
implsres=wc*exp(-wc*t)
//ラプラス逆変換による波形の計算
sumres=10000
//sumres=sum(implsres)
//応答成分の総和
disp(sumres)
scf(10);clf
plot(t,implsres)
DLPF=0
param=implsres'
n=size(param,'r')
for i=0:n-1
DLPF=DLPF+z^(-i)*param(i+1)/sumres
end
DLPFsys=syslin(dt,DLPF);
scf(11);clf;bode([1-DLPFsys],1e-2,4990)
scf(12);clf;bode([DLPFsys],1e-2,4990)
逆変換によるインパルス応答はこれです。
もっともらしいです。
離散領域の LPF ボーデ線図
離散領域の HPF ボーデ線図
でもダメですね。インパルス応答のサンプリングの問題なんでしょうか。先ほどの説明の通り FIR フィルタの係数の総和が 10,000 になっていないといけないのですが、調べてみると 10,299.355 とオーバーしていました。
時間軸が有限なので不足するのは理解しやすいのですが、このようなオーバーする原因は直感的にはたとえば時刻ゼロで 100、時刻 1 で 90 となる波形を離散領域で扱う際、 0 ～ 1 の間は 100 と見ているため、このような減衰特性の場合はオーバー気味に誤差を持つためではないかと想像しています。この辺りの量子化誤差の検討方法はよくわかりません。
で、例によって今度は総和が 10,000 になるように係数全体に補正を掛けます。
スクリプトの //sumres=sum(implsres) を有効にして、係数全体を割ります。
離散領域の LPF ボーデ線図
離散領域の HPF ボーデ線図
こんどは良さそうですね。
なんだか手間ばっかり掛かって大変でしたが、もともとデジタルフィルタというのは連続領域では実現が難しい（できない）フィルタが作れるというのが売りだと思うので、こんな簡単なフィルタだと労多くして何とやらです。
サーボ特性の位相補正やゲイン補正などの一次系のフィルタは IIR で作った方が簡単で世話なしだと思います。
オーディオや通信信号波形などの信号処理のように急峻なカットオフ特性を持たせたり、位相直線性が必要だったりなど波形歪みを問題にしながら、望む特性を得たいなら FIR の方が素直な特性が得られるというのは一般的に云われていることのようです。
一次 HPF は源信号から LPF の出力を引き算することで得られます。
s / (1 + Ts) = 1 - 1 / (1 + Ts)
このカラクリの応用例はこちらにもあります。
ピーク検出をしたいけど、微分器は難しい。
http://blogs.yahoo.co.jp/susanoo2001_hero/7543384.html
スクリプトはこれです。
s=poly(0,'s')
z=poly(0,'z')
dt=1d-4
dires=[]
cutoffrate=1/100
fc=1/dt*cutoffrate
wc=2*%pi*fc
LPF=fc/(fc+s/2/%pi)
LPFsys=syslin('c',LPF);
clf(8);scf(8);clf;bode(LPFsys,0.01,1e4);
clf(9);scf(9);clf;bode(1-LPFsys,0.01,1e4);
t=0:dt:dt/cutoffrate
scf(10);clf
plot(t,csim('impuls',t,LPFsys))
dires=[csim('impuls',t,LPFsys)]
DLPF=0
param=dires'
n=size(param,'r')
//param(1)=10000-sum(param)
for i=0:n-1
DLPF=DLPF+z^(-i)*param(i+1)*dt
end
DLPFsys=syslin(dt,DLPF);
scf(11);clf;bode([1-DLPFsys],1e-2,4990)
scf(12);clf;bode([DLPFsys],1e-2,4990)
結果は以下のようになりました。
連続領域の LPF ボーデ線図 連続領域の HPF ボーデ線図
連続領域の LPF インパルス応答
離散領域の LPF ボーデ線図
離散領域の HPF ボーデ線図
FIR による HPF の応答がやっぱり変ですね。1Hz 以下では利得がフラットになってしまい、位相進みもありません。1 - LPF が上手く計算できていないのでしょうか。
と思って LPF の特性をよく見てみると低い周波数での利得が 0dB よりやや低くなっています。完全な 0dB なら 1 - LPF の結果は 0 レベル = - 無限大 のなるはずですので、どうやらそもそも LPF の応答の精度が悪いようです。これもインパルス応答の誤差かも知れない、ということで応答波形をよく見てみると、時刻 0 の出力がやっぱりゼロ。微分方程式から解くと本来は何か数値があるはずです。
実は LPF は波形こそなまらせますが、入力に与えられたエネルギーは後ろにばらまかれるだけで維持されます。従って係数として採用されたインパルス応答数列の総和は入力インパルスのエネルギーと等しくなくてはいけません。ではインパルスのエネルギーとはなんぞや、というとインパルスは時間幅ゼロ、面積 1 ということですから、ここではサンプリング時間が 0.0001 秒ということで、10,000 のエネルギーと言えそうです。
そこで上述のスクリプトの //param(1)=10000-sum(param) を有効にして（ // を外す）計算し直してみます。ここでは param(1) は 328.9626 になっていました。
離散領域の LPF ボーデ線図
離散領域の HPF ボーデ線図
LPF の高域が変な感じになりましたが、低域は 0dB になったようです。ですので HPF ももっともらしくなりました。
これも scilab のインパルス応答の計算誤差を補正することで解決はしました。
scilab：だ～か～ら
では最初からラプラス逆変換を使って伝達関数からいきなりインパルス応答を求めてみたらどうなるでしょうか。
スクリプトはこれです。
s=poly(0,'s')
z=poly(0,'z')
dt=1d-4
dires=[]
cutoffrate=1/100
fc=1/dt*cutoffrate
wc=2*%pi*fc
LPF=fc/(fc+s/2/%pi)
LPFsys=syslin('c',LPF);
clf(8);scf(8);clf;bode(LPFsys,0.01,1e4);
clf(9);scf(9);clf;bode(1-LPFsys,0.01,1e4);
t=0:dt:dt/cutoffrate
implsres=wc*exp(-wc*t)
//ラプラス逆変換による波形の計算
sumres=10000
//sumres=sum(implsres)
//応答成分の総和
disp(sumres)
scf(10);clf
plot(t,implsres)
DLPF=0
param=implsres'
n=size(param,'r')
for i=0:n-1
DLPF=DLPF+z^(-i)*param(i+1)/sumres
end
DLPFsys=syslin(dt,DLPF);
scf(11);clf;bode([1-DLPFsys],1e-2,4990)
scf(12);clf;bode([DLPFsys],1e-2,4990)
逆変換によるインパルス応答はこれです。
もっともらしいです。
離散領域の LPF ボーデ線図
離散領域の HPF ボーデ線図
でもダメですね。インパルス応答のサンプリングの問題なんでしょうか。先ほどの説明の通り FIR フィルタの係数の総和が 10,000 になっていないといけないのですが、調べてみると 10,299.355 とオーバーしていました。
時間軸が有限なので不足するのは理解しやすいのですが、このようなオーバーする原因は直感的にはたとえば時刻ゼロで 100、時刻 1 で 90 となる波形を離散領域で扱う際、 0 ～ 1 の間は 100 と見ているため、このような減衰特性の場合はオーバー気味に誤差を持つためではないかと想像しています。この辺りの量子化誤差の検討方法はよくわかりません。
で、例によって今度は総和が 10,000 になるように係数全体に補正を掛けます。
スクリプトの //sumres=sum(implsres) を有効にして、係数全体を割ります。
離散領域の LPF ボーデ線図
離散領域の HPF ボーデ線図
こんどは良さそうですね。
なんだか手間ばっかり掛かって大変でしたが、もともとデジタルフィルタというのは連続領域では実現が難しい（できない）フィルタが作れるというのが売りだと思うので、こんな簡単なフィルタだと労多くして何とやらです。
サーボ特性の位相補正やゲイン補正などの一次系のフィルタは IIR で作った方が簡単で世話なしだと思います。
オーディオや通信信号波形などの信号処理のように急峻なカットオフ特性を持たせたり、位相直線性が必要だったりなど波形歪みを問題にしながら、望む特性を得たいなら FIR の方が素直な特性が得られるというのは一般的に云われていることのようです。