昨日の続きで、17段移動平均フィルタ。
17段移動平均フィルタのゲイン特性を計算する。横軸は周波数[Hz]。１ｍｓサンプリングとしている。
縦軸は比率。他の移動平均と比較するために直流でのゲインを１に正規化している。
以下は、横軸を拡大。
１７回移動平均では、1000[Hz]/17=58.8[Hz]でゲインがゼロになる。
２０回移動平均では、1000[Hz]/20=50[Hz]でゲインがゼロになる。
詳しく知りたい方は以下の記事参照。
周波数が異なる場合についての考察：
http://robotcontroller.cocolog-nifty.com/blog/2008/11/post-671a.html
ゲイン特性計算式のまとめ：
http://robotcontroller.cocolog-nifty.com/blog/2008/11/post-7dee.html
Link：
プライムモーション社（Windowsで手軽にリアルタイムIO制御）
ひさびさに技術の話。
とあるテーマで、センサのアナログ信号を16ビットADカードでパソコンに取り込むことになった。
当然、WindowsリアルタイムIO制御「MOS Bench（モスベンチ）」でサポートされているCONTEC社16bitADカード「AD16-16U(PCI)EV」を使う。
パラメータ設定で該当チャンネルを許可すれば、勝手に１ｍｓサンプリングをしてくれる。（上の画像の「アナログ入力チャネルの有効(1)無効(0)」の部分）
パラメータ設定で入力範囲を-10Vから+10Vの指定をすれば、16bitデータを［V］に変換してくれる。（上の画像の「アナログ入力チャネルの入力レンジの下限値-10.00、および、アナログ入力チャネルの入力レンジの上限値10.00」の部分）
また、パラメータ設定で移動平均回数を指定すれば、勝手に移動平均をしてくれて、平均後のデータを下図の画面表示しながら、ロギングデータファイルに１２万行保存できる。
以下は、20回移動平均後の出力。この波形は、ロギングデータのcsvファイルをエクセルに読み込んで、編集したもの。
200ms間に12個の波がある。200ms/12=16.7ms。この逆数は、60Hz。ACラインのノイズが入っている。約2.5mVp-p。
1ms、すなわち1000Hzサンプリングなので、1000Hz/60Hz=16.7で、17回の移動平均をしてみる。
ちなみに、1000Hz/17=58.8Hz。
例によって、「MOS Bench」のパラメータ設定画面で、移動平均回数を17回に設定する。（下の画像の「アナログ入力チャネルのフィルタ段数」の部分）
以下は、１７回移動平均後の出力。
ノイズの振幅は約1mVp-p。3LSBにわたって動いている。前述の波形に対して、1/2.5のノイズになった。
目標とする仕様から判断すれば、これで十分。
しかし、200ms間に5個の波がある。200ms/5=40ms。この逆数は、25Hz。
Link：
プライムモーション社（Windowsで手軽にリアルタイムIO制御）
インパルス不変法で一次IIRローパスフィルタを設計して、
y(nT) = x(nT) + e^(-T/τ)*y{(n-1)T}
(式8)
のフィルタ計算式を得た。
詳細の導出については、以下の記事を参照。
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-8a20.html
このブロック図は、以下となる。
これを置き換えると、
となる。係数1-Aを右シフト演算で置き換えると
となる。
この右シフト回数nと、時定数τ、最終値の関係を表に示す。この表では、サンプリング周期T=1msとしている。
例えば、右シフト回数が４回、したがって、1/2^4=1/16=1-Aのとき、時定数τは約15.5msとなる。
この係数のときのインパルス応答は、０に漸近しないで、１５に落ち着いてしまうことに注意したい。これは1/16の整数演算で、15が桁落ちで後段に伝搬しないことによる。詳細は、以下の記事を参照されたい。
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-ad59.html
また、右シフト回数が４回、1/16演算をすると、フィルタ出力の最終値が分母の値16になることにも注意したい。
上記表のエクセルはこちら
「motion_oyaji_time_constant_table.xls」をダウンロード
次に、サンプリング周期Ｔを変えたら、どのような時定数τが得られるか表に示す。
右シフト回数n=4にて、サンプリング周期T=2msならば、時定数τ=31msになる。
エクセルのデータはこちら
「motion_oyaji_time_constant_table_varT.xls」をダウンロード
この表の時定数で満足できないときは、以下の計算ブロック図を試してみてほしい。つまり、Ｂ＝１／８（右シフト３回）、Ｃ＝１／１６（右シフト４回）とすれば、中間の時定数が得られる。、
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
先の記事に引きつづき、簡単１次IIRローパスフィルタの波形を示す。今回は、周波数特性。
ディジタルフィルタ計算式は、
y(nT) = x(nT) - 1/8*y{(n-1)T} + y{(n-1)T}
(式24)
サンプリング周期Ｔ=1msとしたとき、時定数は、約7msになります。
周波数特性の計算式導出方法について、知りたい方は、以下の記事を参照願います。
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms-7ef.html
周波数特性のグラフは、以下になります。同じ時定数のアナログ一次ローパスフィルタのの周波数特性と比較しています。
振幅特性は似ているけれど、位相特性は違いますね。
■振幅特性（リニア周波数軸）
■振幅特性[dB]（リニア周波数軸）
■位相特性[rad]（リニア周波数軸）
■位相特性[deg]（リニア周波数軸）
エクセルのファイルはこちら
「motion_oyaji_freq_IIIR1st_1_8_linear_comp_081123.xls」をダウンロード
つぎに周波数軸をlog軸にした波形も合わせて示します。
■振幅特性（log周波数軸）
■振幅特性[dB]（log周波数軸）
■位相特性[rad]（log周波数軸）
■位相特性[deg]（log周波数軸）
エクセルのデータはこちら
「motion_oyaji_freq_IIIR1st_1_8_log_comp_081123.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
先の記事に引きつづき、簡単１次IIRローパスフィルタの波形を示す。今回は、正弦波入力に対する応答波形。
ディジタルフィルタ計算式は、
y(nT) = x(nT) - 1/8*y{(n-1)T} + y{(n-1)T}
(式24)
時定数は、約7msになる。
以下に示す波形の横軸はサンプリングｎ。サンプリング周期T=1msなら、横軸は時間[ms]。ディジタルフィルタの入力x(nT)は振幅127の正弦波である。ディジタルフィルタ出力y(nT)とグラフ中で比較するディジタルフィルタ入力は、振幅127正弦波に対して、8倍の振幅、すなわち127*8=1016の振幅にしている。8倍の「8」は、「最終値の定理」によって求められる値である。最終値の定理については、以下の記事を参照してください。
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms-820.html
■100ms周期の正弦波に対する出力（整数計算）
出力y(nT)の振幅は 918。入力との比は918 /(127*8) = 0.90になっている。
位相は7ms程度遅れているので、 7ms/100ms*360度 = 25度程度の遅れになっている。
■20ms周期の正弦波に対する出力（整数計算）
出力y(nT)の振幅は 396。入力との比は396 /(127*8) = 0.39になっている。
位相は約3.5ms程度遅れているので、 3.5ms/100ms*360度 = 63度程度の遅れになっている。
■10ms周期の正弦波に対する出力（整数計算）
出力y(nT)の振幅は 212。入力との比は212 /(127*8) = 0.21になっている。
位相は1.5ms程度遅れているので、 1.5ms/10ms*360度 = 54度程度の遅れになっている。
■5ms周期の正弦波に対する出力（整数計算）
出力y(nT)の振幅は 105。入力との比は105 /(127*8) = 0.10になっている。
位相は1.5ms程度遅れているので、 1.5ms/5ms*360度 = 108度程度の遅れになっている。
エクセルのデータはこちら
「motion_oyaji_IIR_1_8_time_sin_081128.zip」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
一次IIIRフィルタの簡単な実装の記事を続ける。以前の
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-ecd2.html
の記事から
一次遅れのIIRフィルタの計算式は
y(nT) = x(nT) + e^(-T/τ)*y{(n-1)T} (式8)
となることがわかった。この式を変形すると
y(nT) = x(nT) - {1-e^(-T/τ)}*y{(n-1)T} + y{(n-1)T}
(式9)
となる。ブロック図で書くと
となる。今までの記事では、係数{1-e^(-T/τ)}を
{1-e^(-T/τ)} = 1/16
すなわち、右シフト演算４回として、実装した場合の各種波形グラフを見てきた。ディジタルフィルタの計算式は、
y(nT) = x(nT) - 1/16*y{(n-1)T} + y{(n-1)T}
(式11)
である。このとき、たとえばサンプリング周期T=1msとすると、時定数は約15msになる。
◆インパルス応答の時間波形（整数計算）
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-ad59.html
◆ステップ応答の時間波形（整数計算）
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms-820.html
◆正弦波入力応答の時間波形（整数計算）
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms-c85.html
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms2-ac.html
◆周波数応答（リニア周波数軸、整数計算）
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-58eb.html
◆周波数応答（ログ周波数軸、整数計算）
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-000f.html
今回は、時定数を変えてみよう。たとえばサンプリング周期T=1msのとき、時定数は約7msになるようにする。係数としては、1/16を1/8にする。右シフト３回で実装することを目論んでいる。
{1-e^(-T/τ)} = 1/8
ディジタルフィルタの計算式は
y(nT) = x(nT) - 1/8*y{(n-1)T} + y{(n-1)T}
(式24)
となる。
正確な時定数τは
e^(-T/τ) = 1-1/8 = -0.875
にて、サンプリング周期T=1msとすれば、時定数
τ=7.488876[ms]
である。
さて、波形を示す。
■インパルス応答の時間波形（整数計算）
127のステップ入力の場合の応答である。時定数τ=7.488876[ms]のアナログローパスフィルタの応答と比較している。横軸は、サンプリングｎ。サンプリング周期T=1msとすれば、横軸は時間[ms]。
ディジタルフィルタ出力y(nT)の場合は、ゼロに近づかない。７で飽和する。途中の1/8で7が桁落ちでゼロになってしまい、演算結果に反映されず、残ってしまう。
■ステップ応答の時間波形（整数計算）
ディジタルフィルタは127のステップ入力を入れると、最終値は127*8=1016になる。最終値がなぜ8倍されるかというのは、「最終値の定理」から決まる。最終値の定理については、
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms-820.html
の記事を参照。
グラフでは、1016のステップ入力のアナログ１次ローパスフィルタ出力と比較している。
ほぼ同じになる。
エクセルデータはこちら
「motion_oyaji_IIR_1_8_time_step_081121.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
超簡単な一次のIIRフィルタ
y(nT) = x(nT) - 1/16*y{(n-1)T} + y{(n-1)T}
(式11)
に連続したパルスを入力してみる。計算の簡単化のために、係数1/16を使用している。サンプリング周期Ｔ=1msの場合、時定数τは約15msになる。設計の詳細については、前の記事
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-8a20.html
および、係数1/16の背景については、
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-ad59.html
を参照してほしい。
さて、整数演算のため、途中の1/16で桁落ちするので、パルスは127の大きさのパルスを入力する。横軸はサンプリングｎ、サンプリング周期T=1msの場合は時間[ms]になる。
■40ms周期でパルスを入力した場合の出力y(nT)
40ms毎に1msの127の値が入力されるので、平均値127*1ms/40msに、最終値16を掛けた値と比較している。
最終値については、記事
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms-820.html
の中ほどを参照してほしい。
■20ms周期でパルスを入力した場合の出力y(nT)
平均値127*1ms/20msに対して、最終値16を掛けた値と比較している。
■10ms周期でパルスを入力した場合の出力y(nT)
平均値127*1ms/10msに対して、最終値16を掛けた値と比較している。
■5ms周期でパルスを入力した場合の出力y(nT)
平均値127*1ms/5msに対して、最終値16を掛けた値と比較している。
出力y(nT)は、ローパスフィルタが効いて、連続パルス入力の平均値を得ていることがわかる。読者のみなさん、使えそうですか？
エクセルのデータはこちら
「motion_oyaji_IIR_1_16_time_conti_pulse_081121.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
簡単一次IIRローパスフィルタに±パルスを入力した応答波形を調べてみよう。
モータ制御していて、位置制御でサーボロックをかけたとき、このような現象が起こる。
フィルタの計算式は
y(nT) = x(nT) - 1/16*y{(n-1)T} + y{(n-1)T}
(式11)
この式の背景は、
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-8a20.html
および
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-ecd2.html
の記事を参照してください。
整数演算で1/16を行っているため、桁落ちが発生する。この影響を小さくするために、+127、-127のパルスを入力している。
サンプリング周期Ｔ=1msの場合、時定数τは約15msになる。
■40ms毎、+127、-127のパルスを入力した場合の出力（整数演算）
横軸はサンプリング。サンプリング周期Ｔ=1msの場合、時間[ms]になる。
■80ms毎、+127、-127のパルスを入力した場合の出力（整数演算）
飽和している値は１５。1/16の整数計算をしているので、15が1/16でゼロになってしまい、残ってしまう。詳細の記事は、インパルス応答の記事
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-ad59.html
を参照してください。
エクセルのデータはこちら
「motion_oyaji_IIR_1_16_time_both_pulse_081121.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
先の記事で一次IIRローパスフィルタとアナログフィルタの周波数特性を比較した。横軸の周波数はリニアであった。
今回は、周波数をlog軸にしてみる。アナログフィルタの周波数特性に慣れている人は、横軸がlogのほうがわかりやすいと思う。
■ゲイン特性(1Hzから1000Hz)
■ゲイン特性[dB](1Hzから1000Hz)
■位相特性[rad](1Hzから1000Hz)
■位相特性[deg](1Hzから1000Hz)
エクセルのデータはこちら
「motion_oyaji_freq_IIIR1st_linear_comp_081123.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
先の記事で、インパルス不変法で設計した一次IIRローパスフィルタの振幅特性や位相特性のグラフを紹介した。
ディジタルフィルタの計算式は、
y(nT) = x(nT) - 1/16*y{(n-1)T} + y{(n-1)T}
(式11)
で、係数を1/16として、シフト演算を使えるようにしている。
サンプリング周期T=1msの場合、一次IIRフィルタの時定数は、15.4946[ms]となる。詳細は
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-ad59.html
を参照。
今日は、そのグラフを同じ時定数のアナログローパスフィルタと比較してみる。さて、一次ローパスフィルタの伝達関数は、時定数をτ[s]として
1
Ｇ(s) = ------ (式20)
1+sτ
s=jωを代入して、ωの関数にする。
1
Ｇ(jω) = -------
1+jωτ
1
= -------------- * e^{jφ(ω)} (式21)
√ {1+(wτ)^2}
ここで、位相は
φ(ω) = - arctan(wτ) (式22)
絶対値の部分をとって、ゲインは、
1
g(ω) = ----------- (式23)
√ {1+(wτ)^2}
時定数τは、前述のように
τ = 15.4946E-3
同じ時定数のアナログローパスフィルタの位相特性(式22)、ゲイン特性(式23)が得られたので、先の記事
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms-7ef.html
で導出した一次IIRローパスフィルタの振幅特性や位相特性を比較してみる。
グラフは
・細かい部分がわかるように、0から200Hzのグラフ
・サンプリング周波数1000Hzまでがわかるように、0から1000Hzのグラフ
を示す。
■ゲイン特性（0から200Hz）
比較のため、ディジタルのゲイン特性を正規化して（直流でのゲインの値、１６ですべてのｹﾞｲﾝ値を割って）いる。
200Hzまでは、ほとんど差がない。
■ゲイン特性（0から1000Hz）
サンプリング周期T=1msとしているので、1000Hzはサンプリング周波数である。
サンプリング周波数1000Hzの半分、500Hzからずれてくる。
■ゲイン特性[dB]（0から200Hz）
次に縦軸をデシベル[dB]にしてみる。
■ゲイン特性[dB]（0から1000Hz）
サンプリング周期T=1msとしているので、1000Hzはサンプリング周波数である。
■位相特性[rad]（0から200Hz）
次は、位相特性をみてみる。
20Hzくらいまでおなじ。ディジタルはアナログより、位相が遅れない。
■位相特性[rad]（0から1000Hz）
サンプリング周期T=1msとしているので、1000Hzはサンプリング周波数である。
■位相特性[deg]（0から200Hz）
次に縦軸を度[deg]にしてみる。
■位相特性[deg]（0から1000Hz）
サンプリング周期T=1msとしているので、1000Hzはサンプリング周波数である。
エクセルのデータはこちら
「motion_oyaji_freq_IIIR1st_linear_comp_081123.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
今日は、一次のIIRフィルタの周波数特性のグラフを示す。
ゲイン特性の計算式は、先の記事より
1
g(ω) = -----------------------------
√ { (1-AcosωT)^2 +(AsinωT)^2}
(式18)
位相特性の計算式は、先の記事より
φ(ω) = - arctan ( AsinωT / 1-AcosωT) [rad]
(式19)
ここで、
A = e^(-T/τ) (式13)
Ｔはサンプリング周期、τは時定数。
以下のグラフでは、
１－Ａ = 1 - e^(-T/τ) = 1/16
としてあります。
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-ad59.html
の記事の
1-e^(-T/τ) =1/16 (式10)
を意味しています。
ここで、サンプリング周期Ｔ=1msとすれば、時定数τは約15msとなります。
■ゲイン特性
周波数が上がるとゲインが小さくなります。周波数が高い正弦波をこのフィルタに入れると、出力波形の振幅が小さくなります。ローパスフィルタを設計したのだから、当然ですね。実際の入力の正弦波と出力の正弦波については、先の記事
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms-c85.html
を見てください。
■ゲイン特性 [dB]
デシベル[dB]に直しました。
■位相特性 [rad]
周波数が上がると位相が次第に遅れます。しかし、さらに周波数が上がると位相遅れが小さくなります。おもしろいですね。
■位相特性 [deg]
度で計算しました。
エクセルのデータは、こちら
「motion_oyaji_freq_IIIR1st_linear_081123.xls」をダウンロード
今回、サンプリング周期Ｔ=1msとして、時定数τが約15msになっていますが、
サンプリング周期を長くして、Ｔ=10msとすれば、時定数τが約150msとなります。周波数特性のグラフは同じ形になりますが、横軸の周波数最大値は200Hzを20Hzとして考えてください。
逆にサンプリング周期を短くして、Ｔ=0.1msとすれば、時定数τが約1.5msとなります。周波数特性のグラフは同じ形になりますが、横軸の周波数最大値は200Hzを2000Hzとして考えてください。
今回のディジタルフィルタの計算式は、
y(nT) = x(nT) - 1/16*y{(n-1)T} + y{(n-1)T}
(式11)
です。乗算器や浮動小数点演算を使用しないで、シフト演算だけで一次IIRフィルタを実現しました。これなら、安いマイコンやVHDL、verilogで実現できますよね。
これから、今回のサンプリング周期Ｔ=1ms、時定数τ約15ms以外についても考えていきます。
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
さて、今回は、インパルス不変法で設計した一次のIIRローパスフィルタの周波数特性計算式を導きます。
本のディジタルフィルタの計算式は、先の記事
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-ecd2.html
から
y(nT) = x(nT) - {1-e^(-T/τ)}*y{(n-1)T} + y{(n-1)T}
(式9)
である。この式において、サンプリング周期T=1ms、時定数τを約15msとおいて、各種の入力波形に対して、出力波形を計算した。
インパルス入力に対する出力波形
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-ad59.html
ステップ入力に対する出力波形
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms-820.html
正弦波入力に対する出力波形
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iirt1ms15ms-c85.html
を計算してきた。
さて、(式9)の伝達関数は、先の記事
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-8a20.html
から
1
Ｇ(z)=------------------ (式5)
1-e^(-T/τ)*z^(-1)
この式において、z=e^(jωT)、あるいは 、z^(-1)=e^(-jωT)を代入すると、ωTに関わる式
1
Ｇ{e^(jωT)} = ------------------ (式10)
1-e^(-T/τ)*e^(-jωT)
が得られる。オイラーの公式
e^(-jθ) = cosθ - j * sinθ (式11)
より、
e^(-jωT) = cos(ωT) - j * sin(ωT)
であるから、
1
Ｇ{e^(jωT)} = -------------------------
1-A*{ cos(ωT) - j * sin(ωT) }
(式12)
ここで、
A = e^(-T/τ) (式13)
と置いた。さらに、分母の実数部と虚数部(jの付いている項）を分けて
1
Ｇ{e^(jωT)} = -------------------------
1-A*cos(ωT) ＋ j *A*sin(ωT)
(式14)
となる。分母の
実数部 ＋ ｊ 虚数部
の形を
絶対値＊ｅ^(j位相）
の形に直す。この形を
Abs * e^(j *Angle）
と置けば
Abs = √ { (実数部^2) + (虚数部^2) ｝
Angle = arctan ( 虚数部 / 実数部 ) [rad]
で得られる。(式14)の分母は
実数部 = 1-A*cos(ωT)
虚数部 = A*sin(ωT)
なので
Abs = √ { (1-AcosωT)^2 +(AsinωT)^2}
(式15)
Angle = arctan ( AsinωT / 1-AcosωT) [rad]
(式16)
となる。
1 1
-------------- = --- * e^{j(-Angle)}
Abs * e^(j *Angle） Abs
なので、(式14)は、
1
Ｇ{e^(jωT)} = --- * e^{j(-Angle)}
Abs
(式17)
となる。周波数特性のゲイン特性（振幅特性）は、(式17)の絶対値の部分から
1
g(ω) = ----
Abs
1
= -----------------------------
√ { (1-AcosωT)^2 +(AsinωT)^2}
(式18)
デシベル[dB]表示にしたいときは、
20 * log 10 | g(ω) | [dB]
を計算する。
次に、位相特性は、(式17)の位相の部分から
φ(ω) = -Angle
= - arctan { AsinωT / (1-AcosωT) } [rad]
(式19)
度[deg]の表示にしたいときは
φ(ω) / (2π) * 360 [deg]
を計算する。
ただし、
A = e^(-T/τ) (式13)
次の記事では、周波数特性（ゲイン特性、位相特性）を計算してみよう。
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
昨日の続きで、一次IIRローパスフィルタに正弦波を入力して、その入力と出力の比較波形をしめす。ディジタルフィルタの計算式は、
y(nT) = x(nT) - 1/16*y{(n-1)T} + y{(n-1)T}
(式11)
横軸はサンプリングｎ、サンプリング周期T=1msなら、時間[ms]となる。サンプリング周期T=1msのとき、時定数が約15msのローパスフィルタになる。細かい条件は、昨日の記事を見てください。
■周期20msの正弦波入力に対する出力（実数演算）
■周期20msの正弦波入力に対する出力（整数演算）
出力振幅は2032の期待に対して、408くらいであり、約1/5になっている。また、位相は4msほど遅れている。周期20msに対して、位相遅れを計算すると、4ms/20ms*360°=72°程度になる。
後日、アップロード予定のゲイン特性の50Hz（周期20ms）計算値は、0.20であり、上記の約1/5と一致している。位相特性の50Hz計算値は、-69°で、上記の72°程度と一致している。
■周期10msの正弦波入力に対する出力（実数演算）
■周期10msの正弦波入力に対する出力（整数演算）
振幅は2032の期待に対して、313くらいであり、約0.15倍になっている。位相は2msほど遅れている。周期10msに対して、位相遅れを計算すると、2ms/10ms*360°=72°程度になる。
後日、アップロード予定のゲイン特性の100Hz（周期10ms）計算値は、0.10であり、上記の約0.15と一致している。周期が短くなって、波形が粗くなってくるので、これくらいの誤差はある。位相特性の100Hz計算値は、-66°で、上記の72°程度と一致している。
■周期5msの正弦波入力に対する出力（実数演算）
■周期5msの正弦波入力に対する出力（整数演算）
周期5msなので、１周期内に５サンプリングしかない状態で、非常に粗い波形になっている。
振幅は2032の期待に対して、104くらいであり、約0.05倍になっている。位相は0.5msほど遅れている。周期5msに対して、位相遅れを計算すると、0.5ms/5ms*360°=36°程度になる。
後日、アップロード予定のゲイン特性の200Hz（周期5ms）計算値は、0.055であり、上記の約0.05倍と一致している。位相特性の200Hz計算値は、-51°で、上記の36°程度に近い。
エクセルの計算データは、こちら
「motion_oyaji_IIR_1_16_time_sin_081119.zip」をダウンロード
圧縮してあります。
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
最近の記事では、一次のIIRフィルタで簡単なローパスフィルタを設計している。計算式は、しばらく前の記事から、
y(nT) = x(nT) - {1-e^(-T/τ)}*y{(n-1)T} + y{(n-1)T}
(式9)
である。計算の簡単化のために、係数1/16を使用して、時定数τは約15msである。
インパルス応答、ステップ応答は、先の記事で計算してみた。
今回は、正弦波入力x(nT)を入れた応答出力y(nT)を計算する。いずれも、サンプリング周期T=1ms、横軸は時間[ms]である。実数演算と整数演算を比較して示す。整数演算の場合の正弦波振幅は127として、整数に丸めた入力を使用している。同様に整数演算1/16の出力も整数に丸めている。また、いずれの波形においても、入力を１６倍にした出力（理想波形出力）と比較している。１６倍の１６は、先の記事で書いたように最終値の定理で求められる。
■周期1000msの正弦波入力に対する出力（実数演算）
■周期1000msの正弦波入力に対する出力（整数演算）
時定数15msのローパスフィルタに対して、1000ms周期を正弦波を入力した場合、出力の振幅、位相ともに入力に対して大きな変化がない。
実数計算、整数計算の大きな差もない。
差異の詳細については、追ってアップロードするエクセルファイルを見てください。
■周期100msの正弦波入力に対する出力（実数演算）
■周期100msの正弦波入力に対する出力（整数演算）
出力振幅は2032の期待に対して、実際は1454であり、約0.721になっている。位相は12msほど遅れている。周期100msに対して、位相遅れを計算すると、12ms/100ms*360°=43°程度になる。
後日、アップロード予定のゲイン特性の10Hz（周期100ms）計算値は、0.717であり、上記の約0.721と一致している。また、位相特性の100Hz計算値は、-42°で、上記の43°程度と一致している。
別の周波数の正弦波入力に対する出力波形は、次の記事で。
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
まず、おさらい。
１次にIIRフィルタで、サンプリング周期T=1ms、時定数τ=15ms程度のローパスフィルタを実装しようとしている。
計算式は、さきの記事より、
y(nT) = x(nT) - 1/16*y{(n-1)T} + y{(n-1)T}
(式11)
インパルス応答を先の記事で計算したので、今日は、ステップ応答を計算してみよう。
ステップ応答とは、たとえると、スイッチが入ったような応答である。いままでゼロだったのに、突然ある値が入力されたら、どのような出力が出てくるかという波形である。
■ステップ応答（実数計算）
まず、入力x(nT)は、ずっと1。式にすると、
x(nT) = 1 (n=0,1,2,,,)
これを入力として、(式11)のディジタルフィルタ出力y(nT)を計算すると１６に近づく。
y(nT) → 16
この値は、最終値の定理として知られている計算式からも得られる。
フィルタの伝達関数をG(z)として、値が１のステップ入力を入力したとき、出力y(nT)の最終値は、
z
lim{y(nT)} = lim (z-1) * G(z) *----
n->∞ z->1 z-1
= lim z*G(z)
z->1
今計算しているディジタルフィルタの伝達関数は、先の記事の（式５）
http://robotcontroller.cocolog-nifty.com/blog/2008/11/iir-8a20.html
より
1
Ｇ(z)=------------------
1-e^(-T/τ)*z^(-1)
である。これを代入して、y(nT)の最終値は
z
lim{y(nT)} = lim ------------------------
n->∞ z->1 1-e^(-T/τ)*z^(-1)
1
= -------
1-e^(-T/τ)
先の記事で、簡単な演算を目論んで
1-e^(-T/τ) = 1/16
と置いたので、y(nT)の最終値は、１6となる。
グラフの比較波形として、アナログ波形を載せるが、最終値を一致させて比較しやすくするために、値１６のステップ入力を入れたアナログ１次ローパスフィルタ出力波形
y(t) = 16 * {1-e^(-t/15.4946E-3)}
を示している。15.4946E-3は、ディジタルフィルタ設計式から得られた時定数である。約15ms。
■ステップ応答（整数計算）
計算の途中、1/16で桁落ちするので、x(nT)として、127を入力している。式としては、
x(nT) = 127 (n=0,1,2,,,)
最終値は、127の入力に変わったので、127*16=2032である。
比較用のアナログ波形は、最終値を合わせるために
y(t) = 127 * 16 * {1-e^(-t/15.4946E-3)}
を示している。
おおむね同じ波形になる。横軸はサンプリング、T=1msの場合は、時間[ms]。
差異に興味ある方は、エクセルのデータをみてね。データはこちら、
「motion_oyaji_IIR_1_16_time_step_081116.xls」をダウンロード
最終値の定理は、例えば、
横山、佐川、貴家の訳；「ディジタル制御システム 解析と設計」、日刊工業新聞社、1990年、33ページ
に載っている。
Link：
プライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
先の記事のおさらい。
１次のIIRのローパスフィルタは、以下の式で、実現できると書いた。
y(nT) = x(nT) - {1-e^(-T/τ)}*y{(n-1)T} + y{(n-1)T}
(式9)
ここで、例として、係数
1-e^(-T/τ) =1/16 (式10)
としてみよう。4回右シフトで1/16を実現しようともくろんでいる。
1 - 1/16 = e^(-T/τ)
両辺のlnをとって、
-T/τ= ln(1-1/16)
τ = -T /{ln(1-1/16)}
=T * 15.4946
となる。サンプリング周期T=1msとすれば、時定数τ=15.49msのディジタルフィルタが実現できる。 (式9)は、(式10)と置いたので
y(nT) = x(nT) - 1/16*y{(n-1)T} + y{(n-1)T}
(式11)
となる。これを、実際に計算してみよう。
■インパルス応答時間波形（実数計算）
まず、入力x(nT)として、パルスを入力してみよう。
最初の１サンプリングだけ、１。あとは、ずっと0。式にすると、
x(nT) = 1 (n=0)
x(nT) = 0 (n=1,2,3,,,)
この場合の出力y(nT)は、以下のグラフになる。1/16の計算で小数点以下が発生するが、実数で計算している。アナログフィルタ出力の理想波形
y(t) = e^(-t/0.0154946)
と比較している。
横軸はサンプリングｎ。サンプリング周期Ｔ＝１ｍｓの場合は、時間[ms]と考えて良い。アナログ波形とディジタル演算波形に、ほとんど差がない。どのくらい差があるか知りたい方は、追って、エクセルデータをアップロードするのでみてください。
■インパルス応答時間波形（整数計算）
次に整数計算してみる。途中の1/16で桁落ち（小数点以下は無視される）ので、入力x(nT)を大きい値にする。ここでは、127にする。つまり、最初の１サンプリングだけ、127。あとは、ずっと0。式にすると、
x(nT) = 127 (n=0)
x(nT) = 0 (n=1,2,3,,,)
この場合の出力y(nT)は、以下のグラフになる。アナログの理想波形と比較している。横軸はサンプリングｎ。サンプリング周期Ｔ＝１ｍｓの場合は、時間[ms]と考えて良い。
出力y(nT)が0にならないことがわかる。１５で落ち着いてしまう。理由は、途中の1/16演算に15が入力されても、演算結果が0になってしまい、出力に現れないことによる。
つまり、1/16を途中に入れると、16より１少ない15が誤差として残る。
使えるか？
16bit（-32768から+32767）等の有限桁数のデバイスで実装を狙っているので、誤差はやむを得ないという前提である。誤差15が気にならないように１６倍したデータをx(nT)に入力すれば、誤差の割合は減る。つまり、精度は上がる。
あとで説明するかもしれないが、1/32を途中に入れると、32より１少ない31が誤差として残る。同様に1/8を途中に入れると、8より１少ない7が誤差として残る。
なお、今は、サンプリング周期T=1msとして、時定数τ=15.49msとなっているが、
サンプリング周期を1/10にして、T=0.1msなら、時定数は、1/10の1.549ms、
サンプリング周期を10倍にして、T=10msなら、時定数は、10倍の154.9ms、
どんなサンプリング周期でも15.49倍である。
次の記事では、ステップ応答の波形を計算する。
Link：
プライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
前の記事のおさらいから。
一次遅れのIIRフィルタの計算式は
y(nT) = x(nT) + e^(-T/τ)*y{(n-1)T} (式8)
となることがわかった。さて、係数のe^(-T/τ)は、どんな値になるだろうか？たとえば、サンプリング周期Ｔ＝１ｍｓ、設計したいフィルタの時定数τ＝１０ｍｓとすると、e^(-0.1)=0.904837となる。
0.904837をどのように扱うか？
ためしに、この値を１から引いてみよう。
1 - 0.904837 = 0.095163
約1/10.5。1/16としてしまえば、右シフト演算４回で行ける。割り算よりシフト演算なら格段に速い。VHDLやVerilogでも簡単に記述実装できる。
そうか！
係数をe^(-T/τ)ではなくて、１－e^(-T/τ)の形になるようにすれば簡単な演算になるのだ。（式8）を係数１－e^(-T/τ)になるよう変形してみる。
y(nT) = x(nT) - {1-e^(-T/τ)}*y{(n-1)T} + y{(n-1)T}
(式9)
これをブロック図に書くと、
実装例、入出力時間波形、周波数特性については、次回以降で説明します。お楽しみに！
Link：
プライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
さて、FIRフィルタの次は、IIRフィルタを設計してみる。
FIR（Finite Impulse Response）フィルタ。有限インパルス応答フィルタ。
IIR（Infinite Impulse Response）フィルタ。無限インパルス応答フィルタ。
例えば、IIRフィルタの入力に１サンプリングだけ１、その次からずっと０を入力した場合、IIRフィルタ出力は小さな値になっていくが無限に続く。その波形を示す。一次遅れの例である。
これがＩＩＲ（無限のインパルス応答）なる名前のゆえん。
さて、IIRフィルタ設計法にもいろいろあるが、わかりやすい「インパルス不変法」で一次遅れフィルタを設計してみる。
「インパルス不変法」というのは、アナログフィルタ入力に細いパルスを入れた場合のアナログフィルタ出力と上記のように、最初のみ１で、次からずっと０のデータを入力したIIRフィルタ出力波形とを一致させようという方法である。
例えば、アナログフィルタが一次遅れでその入力に細いパルスを入れた時の出力が上記波形のように、指数関数的で減少し、その時定数がτ[s]であるとすれば、出力波形は以下の式になる。
g(t) = e^(-t/τ) ( t>=0 ) (式１)
この波形をT[s]ごと、サンプリングした場合、g(t)の値は、上式にt=0, t=T, t=2T, t=3T,,,を代入して、その値は
1, e^(-T/τ), e^(-2T/τ), e^(-3T/τ) , , ,
と変化していく。 一般には
g(nT) = e^(-nT/τ) ; n=0,1,2,,, （式２）
と表現できる。(式2)のz変換する。右辺は、１サンプリング遅れるごとにその値にz-1をかけて、その総和を求める。数列のz変換は級数になる。左辺は単純にG(z)とおく。
G(z) = 1 + e^(-T/τ)*z(-1) + e^(-2T/τ)*z^(-2) + e^(-3T/τ)*z^(-3) + ,,,
(式3)
これは、公比e^(-T/τ)*z(-1)の等比級数である。（式３）の両辺にe^(-T/τ)*z(-1)をかける。
G(z)*e^(-T/τ)*z(-1) =
e^(-T/τ)*z(-1) + e^(-2T/τ)*z^(-2) + e^(-3T/τ)*z^(-3) + ,,,
(式4)
(式3)の両辺から(式4)の両辺を引くと
G(z) - G(z)*e^(-T/τ)*z(-1) = 1
G(z)について整理して
1
G(z) = ----------------- （式5）
1 - e^(-T/τ)*z(-1)
これが、時定数τの一次遅れになるIIRフィルタである。
G(z) = Y(z) / X(z) を代入して、フィルタ入力Ｘ(z）とフィルタ出力Y(z)の関係式に変形する。
Y(z) 1
--- = ----------------- （式6）
X(z) 1 - e^(-T/τ)*z(-1)
これをY(z)=の形に変形すると、
Y(z) = X(z) + e^(-T/τ)*Y(z)*z(-1) （式7）
これを逆ｚ変換して、実時間の式へ戻す。Y(z) => y(nT) , X(z) => x(nT) , Y(z)*z^(-1) => Y{(n-1)T}と代入して、
y(nT) = x(nT) + e^(-T/τ)*y{(n-1)T} (式8)
この計算式で一次遅れのIIRフィルタができる。ブロック図は、
図としては簡単だが、Ａの値が１と０．９の間の値になる。実際に、Ｔ＝１ｍｓとして、時定数τ＝１０ｍｓとすると、e^(-0.1)=0.904837となる。これを浮動小数点演算ユニット（FPU）や乗算器（かけ算を行うハードウェア）がない安い１６ビットマイコンや３２ビットマイコン、あるいはディジタル回路用ハードウェア記述言語のVHDLやverilogでどのように実装したらよいのか、工夫が必要。工夫については、次の記事で。
Link：
プライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
モータの軸が振動しているとき、すなわち、ロータリーエンコーダを入力がアップダウンするときの移動平均フィルタによる速度検出フィルタの動きをみてみたい。
・周期４０ｍｓで、±１の振動をしている場合
・周期２０ｍｓで、±１の振動をしている場合
いずれの周期の振動も、振幅が小さく、入力データｘ（nT）の値が小さい領域での移動平均フィルタ動作なので、入力ｘ（nT）と移動平均フィルタ出力y(nT)の明確な差異がわからない。
エクセルのデータはこちら「10081112.xls」をダウンロード
シミュレーションのブロック図は、こちら
http://robotcontroller.cocolog-nifty.com/blog/2008/11/post-8678.html
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
今回は、移動平均フィルタへの入力を変えてみる。いままでの記事では、正弦波入力波形を示してきた。今回は、周波数検出、速度検出を意識して、パルスを入力してみる。
一定周波数のパルスをカウントして、一定周期Ｔ毎、その差分を取る(引き算をする)。
x(n) = c(nT) - c{(n-1)T}
これは、周波数検出になる。もし、そのパルス発生源がロータリエンコーダやリニアスケールならば、速度検出になる。このカウンタの差分データを１０回移動平均に入力してみる。
サンプリング周期T=1ms、すなわち、サンプリング周波数1000Hzとして、いろいろな周波数を入力し、１０回移動平均入力のx（ｎT)と、移動平均出力のy(nT)を以下に比較する。横軸は時間[ms]、縦軸はx（ｎT)もしくはy(nT)。
・10010Hzのパルス入力に対するx（ｎT)とy(nT)
入力x(nT)、出力y(nT)ともの１のずれがある。入力x(nT)では平均10に対して１、出力y(nT)では平均100に対して１、なので、フィルタを通したことによりS/N比（信号とノイズの比）が１０倍改善されている。
・1010Hzのパルス入力に対するx（ｎT)とy(nT)
入力x(nT)、出力y(nT)ともの１パルスのずれがある。入力x(nT)では平均1に対して１、出力y(nT)では平均10に対して１、なので、フィルタを通したことによりS/N比（信号とノイズの比）が１０倍改善されている。
・110Hzのパルス入力に対するx（ｎT)とy(nT)
入力x(nT)、出力y(nT)ともの１パルスのずれがある。入力x(nT)では0に対してたまに１、出力y(nT)では1に対してたまに１、なので、フィルタを通したことによりS/N比（信号とノイズの比）が改善されている。
・11Hzのパルス入力に対するx（ｎT)とy(nT)
入力x(nT)、出力y(nT)ともの１パルスのずれがある。入力x(nT)では0に対してたまに１、出力y(nT)では0に対してたまに１、なので、フィルタを通したことによりS/N比（信号とノイズの比）改善効果は、ない。
・まとめ
”１０回”移動平均にすると、サンプリング周波数１０００Ｈｚに対して、その”１／１０”の１００Ｈｚくらいまで、ノイズ改善の効果がある。”１／１００”の１０Ｈｚ前後になると、効果がない。
エクセルの元データはこちら「10081111.xls」をダウンロード
本記事では、定速時の動作について計算した。＋１、－１、＋１と振動している場合の動作については、次の記事で。
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
移動平均回数をどのように決めたらよいか？
移動平均回数を増やすメリット：
データの分解能が増える。
シャープなローパスフィルタが得られる。
移動平均回数を増やすデメリット：
高い周波数の信号（速い動き）を検出できない。
演算回数が増える。
必要メモリが増える。
位相特性が遅れる。（動きが遅れて検出される）。
どのように段数決めるか？
答えは、バランスですな！
「答えになってない！」と言われそうですが、、、
IIRフィルタについても、別の機会に説明していく予定。期待してね！
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
単純移動平均フィルタ位相特性について、移動平均回数をパラメータとしてまとめる。先の記事の周波数特性の式から、位相部分、すなわち、e^（ｊφ）の形のφの部分を取り出す。
・２回平均
φ(ω) = - ωT / 2
・４回平均
φ(ω) = - ωT 3 / 2
・８回平均
φ(ω) = - ωT 7 / 2
・１０回平均
φ(ω) = - ωT 9 / 2
これをグラフにして比較する。サンプリング周期Ｔ＝１ｍｓとしている。横軸は、周波数［Ｈｚ］。縦軸は[rad]である。
次に、縦軸を[deg]にした位相特性グラフを示す。上のグラフの縦軸に対して、*360/2πを行っている。横軸は同じように周波数[Hz]。サンプリング周期Ｔは１ｍｓ、すなわち、サンプリング周波数は１０００Ｈｚである。
サンプリング周波数が２倍の2000Hzになる場合は、上図横軸右端の500Hzを1000Hzに置き換えてほしい。逆に、サンプリング周波数が１／２の５00Hzになる場合は、上図横軸右端の500Hzを250Hzに置き換える。
エクセルの元データは、こちら「motion_oyaji_freq_MA_Conc_Phase_081107.xls」をダウンロード
移動平均段数増加のメリット、デメリットについては、別の記事で。
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
単純移動平均ゲイン特性のまとめ
先の記事の周波数特性の式から、振幅の部分を取り出して
・２回平均
g(ω) = | 2* cos(ωT/2) |
・４回平均
| sin(2ωT) |
g(ω) = ------------
| sin(ωT/2) |
・８回平均
| sin(4ωT) |
g(ω) = -----------
| sin(ωT/2) |
・１０回平均
| sin(5ωT) |
g(ω) = ----------
| sin(ωT/2) |
これをグラフにする。サンプリング周期T=1msとしている。横軸は周波数[Hｚ］。
２回平均は、周波数ゼロ、すなわち、直流にて、ゲインは２となる。４回平均は、４。８回平均は、８。１０回平均は１０。直流でのゲインを１になるように規格化して、比較しやすくする。つまり、２回平均のゲイン特性は、２で割る。４回平均のゲイン特性は、４で割る。他も同様。
規格化したグラフは、
サンプリング周期Ｔ＝１ｍｓ、すなわちサンプリング周波数は1000Hzである。このとき、
２回平均の場合は、サンプリング周波数の１／２の５００Ｈｚでゲインゼロ
４回平均の場合は、サンプリング周波数の１／４の２５０Ｈｚでゲインゼロ
８回平均の場合は、サンプリング周波数の１／８の１２５Ｈｚでゲインゼロ
１０回平均の場合は、サンプリング周波数の１／１０の１００Ｈｚでゲインゼロ
となる。
例えば、他のサンプリング周波数で、２０００Ｈｚにしたら、同様に
２回平均の場合は、サンプリング周波数の１／２の１０００Ｈｚでゲインゼロ
４回平均の場合は、サンプリング周波数の１／４の５００Ｈｚでゲインゼロ
８回平均の場合は、サンプリング周波数の１／８の２５０Ｈｚでゲインゼロ
１０回平均の場合は、サンプリング周波数の１／１０の２００Ｈｚでゲインゼロ
となる。１／ｎも同様。
次に、縦軸の振幅をｄB（デシベル）単位にする。つまり、上記グラフの縦軸の１０底のlogをとって、20をかける。すると
横軸は、周波数[Hz]。サンプリング周期は１ｍｓ。平均回数を増やしたほうが、シャープなローパスフィルタになる。
エクセルの元データはこちら「motion_oyaji_freq_MA_Conc_Amp_081107.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
単純移動平均周波数特性のまとめ
先の記事の伝達関数Ｇ（ｚ）のｚにz=e^(jωT)を代入して、オイラーの公式を利用して、式を整理する。
・２回平均
G(ωT) = 2* cos(ωT/2) * e^(-jωT/2)
・４回平均
sin(2ωT)
G(ωT) = --------- * e^(-jωT3/2)
sin(ωT/2)
・８回平均
sin(4ωT)
G(ωT) = --------- * e^(-jωT7/2)
sin(ωT/2)
・１０回平均
sin(5ωT)
G(ωT) = --------- * e^(-jωT9/2)
sin(ωT/2)
一般の形のn回平均もわかるよね！
２回移動平均のオイラーの公式を使った変形については、こちらを参照
http://robotcontroller.cocolog-nifty.com/blog/2008/10/post-0e41.html
４回移動平均のオイラーの公式を使った変形については、こちらを参照
http://robotcontroller.cocolog-nifty.com/blog/2008/10/post-2781.html
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
単純移動平均の伝達関数のまとめ
移動平均にはローパスフィルタ（LPF）の役目の他に、分解能を増やすという機能もある。そのため、よく使う、情報を1ビット増やすための2回平均、情報を２ビット増やすための４回平均、情報を３ビット増やすための８回平均、情報を10倍にする10回平均について、まとめる。それ以外の平均回数については、本ブログ記事から推測できますよね！
・２回平均
Ｇ(z) = 1 + z^(-1)
・４回平均
1 - z^(-4)
G(z) = ----------
1 - z^(-1)
・８回平均
1 - z^(-8)
G(z) = ----------
1 - z^(-1)
・１０回平均
1 - z^(-10)
G(z) = ----------
1 - z^(-1)
単純２回移動平均の伝達関数の求め方はこちら
http://robotcontroller.cocolog-nifty.com/blog/2008/10/post-0e41.html
単純４回移動平均の伝達関数の求め方はこちら
http://robotcontroller.cocolog-nifty.com/blog/2008/10/post-2781.html
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
1000Hzサンプリング（制御周期1ms）の４回移動平均フィルタに100Hz,200Hz,250のsin波を入力した場合の、
入力x(nT)、出力y(nT)、内部状態のx{(n-1)T}、x{(n-2)T}、x{(n-3)T}を示す。横軸の単位はms。
x(nT)を正弦波として、
y(nT) = x(nT) + x{(n-1)T} + x{(n-2)T} + x{(n-3)T}
を計算したグラフである。
サンプリング周波数1000Hzの1/10の100Hzの正弦波を4回移動平均フィルタへ入力すると、出力振幅が入力の約3倍になる。
サンプリング周波数1000Hzの1/5の200Hzでは、出力振幅が入力の約１倍。
サンプリング周波数1000Hzの1/4の250Hzでは、出力はゼロ。フィルタリングされる。
たとえば、サンプリング周波数が２倍の2000Hz（制御周期500us）のとき、やはり、サンプリング周波数2000Hzの”1/4”の500Hzでは、出力がゼロになる。
読者のみなさんのシステムのサンプリング周波数が1000Hzのｎ倍だったら、特性はｎ倍の周波数に拡大されると考えてください。1/nでも同様。
”２回”移動平均の場合は、サンプリング周波数の”1/2”の周波数の正弦波を入力すると出力がゼロ。
”４回”移動平均の場合は、サンプリング周波数の”1/4”の周波数の正弦波を入力すると出力がゼロ。
エクセルの元データはこちら「motion_oyaji_wave_MA4_081029.xls」をダウンロード
”２回”移動平均の時間波形はこちらhttp://robotcontroller.cocolog-nifty.com/blog/2008/10/post-1979.html
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
前の記事のおさらいから。
単純４回移動平均の周波数特性は、下式であることがわかった。
sin(2ωT)
G{ｅ^(jωT)} = -------- * e^(-jωT3/2) (式７)
sin(ωT/2)
この式の振幅部分（e＾以外の部分）を切りだして、絶対値にして、ゲイン特性（振幅特性）は、
｜sin (2ωT)｜
g(ω) = ------------- (式８)
｜sin (ωT/2)｜
ω=0のときは、右辺は0/0になるので、その時の値は、分子と分母を各ωで微分して、ω=0を代入して
2T * cos(2ωT)
lim ---------------- = 4 （式９）
ω→0 T/2 * cos(ωT/2）
になる。
詳細は、森口、宇田川、一松：「数学公式Ⅰ－微分積分・平面曲線－」、岩波全書、４２ページ（1981年17刷）を参照してください。いまは、新しい版になっています。
あるいは、公式
sin4θ = 8sinθ(cosθ)^3 - 4sinθcosθ （式１０）
を（式８）の分子に用いて、ゲイン特性を
g(ω) = 8 * {cos(ωT/2)}^3 - 4 * cos(ωT/2) （式１１）
としてもよい。
上述のsin4θの公式は、森口、宇田川、一松：「数学公式Ⅱ－級数・フーリエ解析－」、岩波全書、１８６ページ（1979年17刷）を参照してください。いまは、新しい版になっています。
位相特性は、（式７）の位相部を切りだして
φ(ω) = - 3ωT/2 (式１２）
具体的な例として、サンプリング周期1000Hz（制御周期1ms）、横軸周波数のグラフにする。グラフからサンプリング周期の1/4の250Hzでゲインがゼロ、つまり、250Hzの正弦波信号は出力に現れない。フィルタリングされる。サンプリング周期の1/10の10０Hzの正弦波信号は入力の約３倍になって、出力現れる。0Hz、すなわち、直流は、４倍になって現れる。
エクセルの元データはこちら 「motion_oyaji_freq_MA4_081029.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
次は、４回の単純移動平均をしてみよう。
■図の入出力関係は
y(nT) = x(nT) + x{(n-1)T} + x{(n-2)T} + x{(n-3)T} (式1)
■伝達関数を求める。(式１)の両辺をまずＺ変換する。y(nT) => Y(z) , x(nT) => X(z) , x{(n-1)T)} => X(z)*z^(-1) , x{(n-2)T)} => X(z)*z^(-2) , x{(n-3)T)} => X(z)*z^(-3)と置き換える。すると
Y(z) = X(z) + X(Z)*z^(-1) + X(Z)*z^(-2) + X(Z)*z^(-3)
= X(z) * { 1 + z^(-1) + z ^(-2) + z^(-3) }
次に出力Y(z)とX(z)の比の形に変形する。
Y(z)/X(z) = 1 + z^(-1) + z ^(-2) + z^(-3) (式２)
これが伝達関数である。簡単のため、Y(z) / X(z) = G(z)とおくと、
G(z) = 1 + z^(-1) + z ^(-2) + z^(-3) (式３)
これは、等比級数である。これをまとまった形に変形する。（式３）の両辺にz^(-1)をかけると、
G(z)*z^(-1) = z^(-1) + z ^(-2) + z^(-3) + z^(-4) (式4)
(式３)の両辺から、（式４）の両辺を引いて
G(z) - G(z)*z^(-1) = 1 - z^(-4)
G(z)について整理すると、
1 - z^(-4)
Ｇ(z) = --------- (式５)
1 - z^(-1)
■次に周波数特性を求める。伝達関数の式において、z^(-1)をe^(-jωT)に置き換える。左辺のZは、e^(jωT)に置き換える。ωは角周波数。Tはサンプリング周期。
1 - e^(-jω4T)
G{ｅ^(jωT)} = ------------ (式６)
1 - e^(-jωT)
（式６）の分子は、オイラーの公式
sinθ = { e^(jθ) - e^(-jθ) } / 2j
を使って、
1 - e^(-jω4T) = e^(-jω2T) * { e^(jωT) - e^(-jωT) }
= 2j * e^(-jω2T) * { e^(jωT) - e^(-jωT) } / 2j
= 2j * e^(-jω2T) * sin(2ωT)
となる。
(式６)の分母も同様に、分子に対して4TがTになっただけなので
1 - e^(-jωT) = 2j * e^(-jωT/2) * sin(ωT/2)
となる。したがって、（式６）は、
sin(2ωT)
G{ｅ^(jωT)} = -------- * e^(-jωT3/2) (式７)
sin(ωT/2)
これが周波数特性である。
■参考
２回移動平均の周波数特性はこちら
http://robotcontroller.cocolog-nifty.com/blog/2008/10/post-0e41.html
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
２回移動平均の入出力波形を示す。サンプリング周波数は1000Hz。横軸はms。入力信号は、100Hz、900Hz、1100Hzのsin波である。このとき、入力x(z)、出力y(z）ともに同じ波形になる。
これは、
900 = 1000 - 100
1100 = 1000 ＋ 100
のようにサンプリング周波数1000Hzを基準にして、100Hzが±された周波数なので、このような現象が発生する。
1000Hzだけでなく、その整数倍でも発生する。たとえば、1000の２倍の
1900 = 2000 - 100
2100 = 2000 ＋ 100
でも同じになる。
この現象は、ここで例に挙げた移動平均処理とは関係なく、信号をサンプリングすることによって発生する一般的な現象である。
エクセルの元データはこちら「motion_oyaji_wave_MA2_ov_081025.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
1000Hzサンプリング（制御周期1ms）の２回移動平均フィルタに100Hz,200Hz,333Hz,500Hzのsin波を入力した場合の、
入力x(nT)、出力y(nT)、内部状態のx{(n-1)T}
を示す。横軸の単位はms。
x(nT)を正弦波として、
y(nT) = x(nT) + x{(n-1)T}
を計算したグラフである。
サンプリング周波数1000Hzの1/10の100Hzの正弦波を２回移動平均フィルタへ入力すると、出力振幅が入力の約２倍になる。
サンプリング周波数1000Hzの1/5の200Hzでは、出力振幅が入力の約１．５倍。
サンプリング周波数1000Hzの1/3の333Hzでは、出力振幅が入力の約１倍。
サンプリング周波数1000Hzの1/2の500Hzでは、出力はゼロ。フィルタリングされる。
たとえば、サンプリング周波数が２倍の2000Hz（制御周期500us）のとき、やはり、サンプリング周波数2000Hzの”1/2”の1000Hzでは、出力がゼロになる。
読者のみなさんのシステムのサンプリング周波数が1000Hzのｎ倍だったら、特性はｎ倍の周波数に拡大されると考えてください。1/nでも同様。
エクセルの元データはこちら「motion_oyaji_wave_MA2_081025.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (1) | トラックバック (0)
前の記事のおさらいから。
２回移動平均の周波数特性は
Y{ｅ^(jωT)} / X{ｅ^(jωT)} = 2 * cos(ωT/2 ) * e^(-jωT/2) (式３)
であることがわかった。この式の振幅部分（e＾以外の部分）を切りだして、絶対値にして、ゲイン特性（振幅特性）は、
g(ω) = ｜2 * cos (ωT/2)｜ (式4)
具体的な例として、サンプリング周期1000Hz（制御周期1ms）、横軸周波数のグラフにする。グラフからサンプリング周期の半分500Hzでゲインがゼロ、つまり、５００Hzの正弦波信号は出力に現れない。フィルタリングされる。サンプリング周期の1/10の10０Hzの正弦波信号は入力の1.9倍になって、出力現れる。また、0Hz、すなわち、直流は、２倍になって現れる。
同じく位相部分（eのカッコの中のjを除いた部分）を取り出して
φ(ω) = -ωT/2 [rad] (式5)
同じく、サンプリング周期1000Hz（制御周期1ms）、横軸周波数のグラフにする。サンプリング周期の1/10の10０Hzの信号は位相が18°遅れて、出力に現れる。
エクセルの元データはこちら「motion_oyaji_freq_MA2_081025.xls」をダウンロード
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
仕事で、移動平均の周波数特性の検討書を作ることになったので、昔勉強した本を引っ張りだして、書いている。「移動平均」は、株価の分析にも出てくる。ここでは、ディジタル信号処理のローパスフィルタとして使用する。まず、簡単な２回移動平均から。FIRフィルタのもっとも簡単な形。
FIR： Finite Impulse Response (有限インパルス応答)
■図の入出力関係は
y(nT) = x(nT) + x{(n-1)T} (式1)
■伝達関数を求める。(式１)の両辺をまずＺ変換する。y(nT)=>Y(z),x(nT)=>X(z),x{(n-1)T)}=>X(z)*z^(-1)と置き換える。すると
Y(z) = X(z) + X(Z)*z^(-1)
次に出力Y(z)とX(z)の比の形に変形する。
Y(z)/X(z) = 1 + z^(-1) (式２)
これが伝達関数である。
■次に周波数特性を求める。伝達関数の式において、z^(-1)をe^(-jωT)に置き換える。左辺のZは、e^(jωT)に置き換える。ωは角周波数。Tはサンプリング周期。株式市場なら、１日。速度制御なら1msから数百usあたり。
Y{ｅ^(jωT)} / X{ｅ^(jωT)} = 1 + e^(-jωT)
これを振幅と位相の式に変形する。複素数に関わる知識が必要になる。 オイラーの公式
e^(jθ) + e^(-jθ) = 2*cosθ
を使えるように、変形する。
{ e^(jωT/2) + e^(-jωT/2)} * e^(-jωT/2 )
上記のオイラーの公式を使って,
Y{ｅ^(jωT)} / X{ｅ^(jωT)} = 2 * cos(ωT/2 ) * e^(-jωT/2) (式３)
これが、２回移動平均の周波数特性である。
Link：
「モーションおやじ」のプライムモーション社（みんなで手軽にWindowsリアルタイムIO制御）
| 固定リンク | コメント (0) | トラックバック (0)
最近のコメント