1 - Researchmap

ガウス誤差関数を利用した
収束の速いヒルベルト変換ディジタルフィルタ
2013年3月6日
吉川浩†，棟方渚† ，小野哲雄†
†北海道大学大学院情報科学研究科
情報処理学会第75回全国大会
ヒルベルト変換とは
Hilbert Transformation
• 下記の周波数特性 H(ω) を持つ変換
1.5
縦軸は虚数
式で表すと H(ω) = -i * sgn(ω)
1
0.5
-5
-3
0
-1
-0.5
1
3
5
-1
-1.5
• 全ての周波数成分の位相を90度ずらす
• A*cos(ωt) のヒルベルト変換は A*sin(ωt)
情報処理学会第75回全国大会
ヒルベルト変換とは
Hilbert Transformation
• ディジタル信号処理には欠かせない
• 解析信号（複素信号）を得るために用いる
• 信号解析では信号を複素数で扱うと便利
• A e-i(ωt+θ) = A cos(ωt+θ) + i A sin(ωt+θ)
瞬時振幅（A）と瞬時位相（ωt+θ）を正確に取り出せる
• 実世界の信号は実数のみ（虚数がない）
• ヒルベルト変換で虚数部を生成
A cos(ωt+θ)  A sin(ωt+θ)
情報処理学会第75回全国大会
ヒルベルト変換フィルタの設計
*)
• 一般的な手法(Parks-McClellan法
* James McClellan and Thomas Parks. 1972
• 反復法によりFIRフィルタ係数を求める
• 最適解が得られる
• 計算が遅い
• 条件によっては反復計算が収束しないことがある
• 提案手法(ガウス誤差関数を用いる方法)
• ガウス誤差関数の逆フーリエ変換から求める
• 実用解が得られる（最適とは限らない）
• 反復法ではないので計算が速い
• 必ず計算結果が得られる
情報処理学会第75回全国大会
提案手法の説明
情報処理学会第75回全国大会
着想
• HTのインパルス応答からFIRフィルタを設計
(周波数領域での積は時間領域での畳み込み)
H(ω)
F-1

h(t)
y(t )   h( x) f (t  x)dx

N
yi   hk  f i  k
方形窓を掛けた
ヒルベルト変換特性
ヒルベルト変換の
インパルス応答波形
k 0
fN
h
Z-1
無限に続くインパルスを有限タップで打ち切る0h0 fN
⇒ 打切り誤差が生じる
fN-1
h1
f
Z-1 N-2
h2
h1 fN-1 h2 fN-2
f1
hn-1
Z-1
f0
hn
hN-1 f1 hN f0
+
y = Σ hk fN-k
情報処理学会第75回全国大会
着想
• 特に，方形のような不連続窓のインパルス
応答波形はなかなか零に収束しない
• 打切り誤差が大きくなる
打切り
打切り
• 方形窓ではなく連続な窓にする
情報処理学会第75回全国大会
下記特性の窓関数を見つける
• 窓関数の周波数領域の特性
• 通過域はできるだけ平坦
• 連続で滑らか
• 窓関数のインパルス応答の特性
• | t | の増加とともに素早く零へ収束
ガウス誤差関数を使った窓
情報処理学会第75回全国大会
ガウス誤差関数 erf(x)
• 定義式
erf( x) 
2


x
0
e t dt
2
「ガウス関数」を積分したもの．
ガウス関数は時間と周波数の分解能に対して最適
• グラフ
1
0.5
0
-5
-3
-1
1
3
5
erfは平坦な特性を持つ
-0.5
-1
情報処理学会第75回全国大会
ガウス誤差関数によるフィルタ
• 提案フィルタの周波数特性
周波数特性 H(ω)
ω1=0のとき
ω1≠0のとき
   
   
  
    
H ( )  i  K erf  2
  erf  2
  erf  1
  erf  1

  
  
  
   
K は定数（利得を正規化する）
ω1
-ω2
ω2
-ω1
ω1, ω2 によってカットオフ周波数を設定できるようにした
σによって遷移域の傾きを設定できるようにした
• 提案式を逆フーリエ変換
h(t )  M  e

 2t 2
4

cos(1t )  cos(2t )
t
Mは定数
急速に減衰する項
インパルス応答 h(t)
打切り
打切り
情報処理学会第75回全国大会
評価
情報処理学会第75回全国大会
ヒルベルト変換フィルタ設計ツール
• パラメータ変更に対してリアルタイムに表示できる
• PM法は解が収束するまで反復計算するのでリアルタイムは無理
結果はC言語で保存される
double hcoeff[] = {
5.3340066754195e-002,
6.4546659217938e-003,
-1.5383950543097e-002,
1.0909027758063e-001,
2.8823052782103e-001,
2.7244132465002e-001,
0.0000000000000e+000,
-2.7244132465002e-001,
-2.8823052782103e-001,
-1.0909027758063e-001,
1.5383950543097e-002,
-6.4546659217938e-003,
-5.3340066754195e-002 };
情報処理学会第75回全国大会
性能評価
• パラメータ（通過域，タップ数）を変えPM法
と提案手法のフィルタ特性を比較
• 用いたフィルタ設計ソフトウェア
• 提案手法：自作プログラム
• PM法：三上直樹氏のフィルタ設計プログラム
• はじめて学ぶディジタル・フィルタと高速フーリエ変換
CQ出版. 2005
情報処理学会第75回全国大会
性能評価 – 実験1
• PM法と提案手法の性能比較 (1)
• TAP数=27，通過域 = [ 0.10， 0.40 ]
方形窓
（偏差
6.4×10-2）
Parks-McClellan法
（偏差
3.4×10-5）
提案手法
（偏差 9.9×10-4）
情報処理学会第75回全国大会
性能評価 – 実験2
• PM法と提案手法の性能比較 (2)
• TAP数=51，通過域 = [ 0.10， 0.40 ]
PM法は
収束せず
方形窓
（偏差
3.6×10-2）
Parks-McClellan法
提案手法
（偏差 1.3×10-5）
情報処理学会第75回全国大会
性能評価 – 実験3
• PM法と提案手法の性能比較 (3)
• TAP数=101，通過域 = [ 0.03， 0.47 ]
方形窓
（偏差
5.8×10-2）
Parks-McClellan法
（偏差
2.3×10-5）
提案手法
（偏差 5.7×10-4）
情報処理学会第75回全国大会
性能評価 – 実験4
• バンドパス特性を持つHTフィルタ
通過域 = [ 0.1125, 0.1375 ]
Tap数 = 15
Parks-McClellan法
通過域
提案手法
通過域
• PM法は上記通過域で収束せず（通過域を広げて収束）
• 収束する条件では利得が１を大きく超える場所がある
情報処理学会第75回全国大会
応用例
情報処理学会第75回全国大会
応用例
音声に低周波データを載せる
• 脈波や心電図などを音声の振幅に載せて
マイク端子から取得する装置を開発
(第16回日本バーチャルリアリティ学会大会，情報処理学会インタラクション2012)
スマートフォンを使った脈波測定器
スマートフォンを使った心電測定器
応用例
AM復調の原理
• ヒルベルト変換を利用して振幅を求める
振幅変調された音声
A(t) cos(ωt+θ)
ヒルベルト変換対を求める
A(t) cos(ωt+θ) ，A(t) sin(ωt+θ)
絶対値を計算する
A(t)
• ディジタル信号処理によって正確な振幅が分かる
• (従来)包絡線検波  サンプリング点が頂点でないと誤差
• ヒルベルト変換  サンプリング点が頂点でなくてもよい
(実信号を解析信号化し瞬時振幅を得る)
情報処理学会第75回全国大会
応用例
復調ブロックダイアグラム
A(t)cos(ωt) + Noise
BPF
|A(t)|2
A(t)sin(ωt)
HTF
pow2
+
↓
A(t)
8kHz  200Hz
A(t)cos(ωt)
pow2
LPF
(Sinc)
HTF
ヒルベルト変換フィルタ（実験４で作ったもの．タップ数=15）
BPF
バンドパスフィルタ（入力信号からノイズとDC成分を除去）
pow2
入力を２乗する
↓
デシメーション・フィルタ（サンプリング周波数を下げる）
入力の平方根を計算する
情報処理学会第75回全国大会
応用例
実際の音声波形と復調データ
振幅変調された音声波形
復調されたデータ
情報処理学会第75回全国大会
まとめ
情報処理学会第75回全国大会
まとめ／Future Work
• まとめ
• erf(x) を利用したHTフィルタ設計法を提案
• 反復法ではないため計算が高速／必ず解が得られる
• 得られるフィルタの性能はPM法より若干劣る
• 実用解が得られる
• スマホで生体情報を取得する装置へ応用
• Future Work
• erf(x)のヒルベルト変換以外への適用
• ウェーブレット変換 etc.
情報処理学会第75回全国大会
ご静聴ありがとうございました

Download Report