Presentation Title

Size: px

Start display at page:

Download "Presentation Title"

みさえあくや
5 years ago
Views:

1 MATLAB で信号処理 ~ 各種センサー信号を題材として ~ MathWorks Japan シニアアプリケーションエンジニア竹本佳充 213 The MathWorks, Inc. 1

2 Agenda 1. センサーシステム設計の概要 2. センサーデータの解析周波数解析信号の遅延類似性の評価等 3. センサーデータのフィルタリングフィルタの種類と特徴フィルタの設計 4. システムモデリングから実機テスト固定小数点化システムモデルへの展開実装環境 5. まとめ 2

3 Agenda 1. センサーシステム設計の概要 2. センサーデータの解析周波数解析信号の遅延類似性の評価等 3. センサーデータのフィルタリングフィルタの種類と特徴フィルタの設計 4. システムモデリングから実機テスト固定小数点化システムモデルへの展開実装環境 5. まとめ 3

4 概要 : センサーとは? ( 出展 : フリー百科事典ウィキペディア (Wikipedia) ) センサまたはセンサー ( 英 : sensor 英語発音 : /ˈsensɔːr/ センソー /ˈsensə/ センサ ) は自然現象や人工物の機械的電磁気的熱的音響的化学的性質あるいはそれらで示される空間情報時間情報を何らかの科学的原理を応用して人間や機械が扱い易い別媒体の信号に置き換える装置のことをいい概要センサーは今日日常生活のあらゆるところで使用されている特に自動車機械航空宇宙医学産業およびロボット工学などで大いに使用される 4

概要 : センサーの種類物理特性力ストレインゲージ ( ひずみゲージ ), ロードセル, 半導体圧力センサー

非接触式, 放射温度計音声マイクロフォン化学イオン濃度ガス濃度時空間特性変位ポテンショメータ,

(PSD), リミットスイッチ距離超音波距離計, 静電容量変位計, 光学式測距, 電磁波測距速度

地震センサー一次元画像リニアイメージセンサー二次元画像 CCD イメージセンサー, CMOS イメージセンサー

5 概要 : センサーの種類物理特性力ストレインゲージ ( ひずみゲージ ), ロードセル, 半導体圧力センサー光光センサー, 光電素子, フォトダイオード磁気磁気センサー温度接触式, サーミスタ, 抵抗測温体, 熱電対, 非接触式, 放射温度計音声マイクロフォン化学イオン濃度ガス濃度時空間特性変位ポテンショメータ, 差動トランス, 回転角センサー, リニアエンコーダ回転数タコジェネレータ, ロータリエンコーダ位置光位置センサー (PSD), リミットスイッチ距離超音波距離計, 静電容量変位計, 光学式測距, 電磁波測距速度レーザードップラー振動速度計, レーザドップラー流速計角速度ジャイロセンサー加速度加速度センサー, 地震センサー一次元画像リニアイメージセンサー二次元画像 CCD イメージセンサー, CMOS イメージセンサーその他液硬度電場電流電圧電力赤外線放射線湿度におい流量傾斜振動時間出展 : フリー百科事典ウィキペディア (Wikipedia) 5

6 概要 : センサーシステムアナログ回路測定物センサ ADC デジタル信号処理 DAC センサシステムは信号処理技術と密接に関係信号処理の基礎知識は必要最適な信号処理システム設計環境が重要 6

7 概要 : 課題 1 本来の周波数成分が再現されない測定物センサアナログ回路 ADC デジタル信号処理 DAC 1 信号の打ち切りによる影響 cos 波 ( 時間軸 ) 1 cos 波 ( 時間軸 ) 打ち切り 5 折り返し歪による影響時間 [s] cos 波 ( 周波数軸 ) 時間 [s] cos 波 ( 周波数軸 ) 周波数 [Hz] サイドローブなし周波数 [Hz] サイドローブ発生回避する方法は? 7

8 パワー / 周波数 (db/hz) 概要 : 課題 2 所望成分の取得測定物センサアナログ回路所望成分 ADC デジタル信号処理 DAC 不要成分ピリオドグラムパワースペクトル密度推定特定の周波数成分のみ取り出すには? 8

9 概要 : 課題 3 ADC 前処理と DAC 後処理測定物センサアナログ回路 ADC デジタル信号処理 DAC ADC 周波数の 1/2 以上の成分はアンチエイリアシングフィルタで除去 DAC 後のイメージ成分はアンチイメージングフィルタで除去最適なフィルタは? 9

10 概要 : 課題 4 システムシミュレーション実機への展開測定物センサアナログ回路 ADC デジタル信号処理 DAC 要素設計からアーキテクチャ評価へデジタル処理とアナログ処理の混載フィードバックループシステムの安定性解析固定小数点化がシステム性能に与える影響実機評価環境の整備これらを効果的に実現する環境は? 1

11 Agenda 1. センサーシステム設計の概要 2. センサーデータの解析周波数解析信号の遅延類似性の評価等 3. センサーデータのフィルタリングフィルタの種類と特徴フィルタの設計 4. システムモデリングから実機テスト固定小数点化システムモデルへの展開実装環境 5. まとめ 11

12 信号の離散化量子化最小ステップサンプリング周期 Ts[s] time[s] サンプリング周波数 fs[hz] =1/Ts[s] 時間軸方向の離散化 : 標本化 ( サンプリング ) 振幅方向の離散化 : 量子化 12

13 サンプリング定理アナログサンプリングデジタル ( サンプリング密 ) デジタル ( サンプリング粗 ) F time[s] F time[s] F time[s] fs1 fs2 エイリアシング frequency[hz] fs1/2 frequency[hz] fs2/2 frequency[hz] 元信号の fs/2[hz] 以上の成分は折り返し ( エイリアシング ) に影響を受ける 2 倍以上のサンプリング or アンチエイリアシング処理が必要 13

14 離散高速フーリエ変換 :FFT N FFT : X ( k) x( j) where e N ( 2 i )/ N j 1 k 1 ( j 1)( k 1) N IFFT : x( j) (1/ N) X ( k) N ( j 1)( k 1) N time[s] FFT DC 成分周波数分解能 : fs/ データ数折り返し fs/2 fs frequency[hz] データタイプは複素数 ( 振幅と位相の情報を持つ ) fs/2( ナイキスト ) で折り返し変換前と後で同じデータ点数データ点数が多いほど周波数分解能は上がる 14

15 窓関数最終値初期値不連続 time[s] window 最終値初期値初期値と最終値の不連続性 ( 打ち切り誤差 ) を低減 (FFT の前処理 ) 周波数成分の再現性向上 15

16 主なスペクトル推定法カテゴリ特徴メソッド名関数名パラメトリック法ノンパラメトリック法部分空間法スペクトル推定法の種類と特徴信号自体から PSD を求める比較的長いデータに有効整形システムの係数推定により PSD を求める比較的短いデータに有効自己相関行列の固有解析固有分解をベースとする低 S/N 正弦波信号のスペクトル推定に有効ピリオドグラムウェルチマルチテーパユールウォーカーバーグ共分散 MUSIC 固有ベクトル periodogram pwelch pmtm pyuler pburg pcov pmusic peig 16

Ex1_1) 信号の FFT( エイリアシングの確認 ) 右クリックナイキスト以上の成分は折り返す 2

8 1 5 2 8 1 5 2 8 1 5 エイリアシング -5-5 -5-1 1 2 3 4 5-1 1

17 Ex1_1) 信号の FFT( エイリアシングの確認 ) 右クリックナイキスト以上の成分は折り返すエイリアシング f1 = 1, f2 = 5 f1 = 1, f2 = f1 = 1, f2 = 55 17

18 Ex1_2) 信号の FFT( 窓関数の効果 ) Sig >>wintool 拡大 2 ハニングサイドローブの減衰大 ( ハニング窓の特徴 ) メインローブが狭い ( ハミング窓の特徴 ) ハミング 18

19 時間パワー / 周波数 (db/hz) Ex1_3) スペクトログラム FFT ピリオドグラムパワースペクトル密度推定 Sig 周波数成分は求まるが時間の情報が欠落周波数 (khz) x 1-3 スペクトログラム時変信号を効果的に可視化時間 x 1-3 ノッキングセンサの例周波数 (Hz) x

20 遅延量の推定 time[s] x サンプル time[s] 相互相関 (xcorr) x サンプル時間軸信号の類似性の指標信号の遅延量推定に使用可能 2

21 周波数成分の類似性共通成分 frequency[hz] コヒーレンス (mscohere) frequency[hz] frequency[hz] クロススペクトル / 位相推定 (cpsd) 位相遅れ frequency[hz] 周波数軸信号の類似性の指標位相遅れ推定 21

C31 C21 s3 s2 s1 Ex1_4) 遅延量測定相互相関関数の計算.1.5 -.

3 Samples 1 s1 s2 s3 振動センサー信号 Cross-Correlations Sig -35 サンプルずれ 15

22 C31 C21 s3 s2 s1 Ex1_4) 遅延量測定相互相関関数の計算 Samples 1 s1 s2 s3 振動センサー信号 Cross-Correlations Sig -35 サンプルずれ 15 サンプルずれ相互相関 (C21) 相互相関 (C31).5 s2 s1 s3 振動センサの例 Samples 22

s2 P2 s1 P1 Ex1_5) 周波数成分類似度推定 Sig コヒーレンス / クロススペクトルの計算 4 2-2 Time Series.2.15.1 Power Spectrum -4.5-6.2.4.6.8 1 6 4 2-2 -4.2.4.6.8 1 Time (secs) 1 2 3 4 5 共通成分.

23 s2 P2 s1 P1 Ex1_5) 周波数成分類似度推定 Sig コヒーレンス / クロススペクトルの計算 Time Series Power Spectrum Time (secs) 共通成分 Frequency (Hz) Coherence Estimate.75 コヒーレンス Cross Spectrum Phase (deg) 位相 Frequency (Hz) -6 23

24 補足 : 信号解析アプリ (sptool) Sig 時間軸表示スペクトル推定結果 24

25 補足 : File Exchange e キーワード検索でデモをサーチ 25

26 Agenda 1. センサーシステム設計の概要 2. センサーデータの解析周波数解析信号の遅延類似性の評価等 3. センサーデータのフィルタリングフィルタの種類と特徴フィルタの設計 4. システムモデリングから実機テスト固定小数点化システムモデルへの展開実装環境 5. まとめ 26

27 フィルタの種類ローパスフィルタ f(hz) 種類フィルタバンドパスフィルタ f(hz) 実現方法ハイパスフィルタ f(hz) アナログフィルタデジタルフィルタパッシブアクティブ FIR IIR 27

28 デジタルフィルタ :FIR 差分方程式 : N y(n) = Σ b k x(n-k) k= x(n) b FIR フィルタの構成 + y(n) Z -1 x(n-1) b 1 Z -1 b 2 x(n-2) Z -1 y(n) = b *x(n) + b 1 *x(n-1) + b 2 *x(n-2) + + b N *x(n-n) Z -1 x(n-n) b N 28

29 デジタルフィルタ :FIR (cont d) FIR フィルタの伝達関数 k k H(z) = = N b k Z 1 b(1)+b(2)z -1 +b(3) Z b(n+1)z -N 1 ( 次数 :N) 入力 X(z) システム伝達関数 H(z) H( z) Y ( z) X ( z) 出力 Y(z) Z 変換 X ( z) x( n) Z n x( n) X ( z) x( n 1) z 1 X ( z) n 29

30 デジタルフィルタ :FIR/MA 5 点移動平均フィルタの例 x(n) 1/5 + y(n) x(n) + y(n) 1/5 x(n-1) Z -1 1/5 x(n-1) Z -1 Z -1 1/5 Z -1 x(n-2) x(n-3) Z -1 Z -1 1/5 x(n-2) x(n-3) Z -1 Z -1 差分方程式 y(n) = 1/5*{x(n) +x(n-1) + x(n-2) +x(n-3)+x(n-4)} x(n-4) 1/5 x(n-4) 1 1/ 5* 伝達関数 z 2 z z 1 z

31 デジタルフィルタ :IIR 差分方程式 : N y(n) = Σ b k x(n-k) k= N Σ k=1 a k y(n-k) x(n) b IIR フィルタの構成 + + y(n) x(n-1) Z -1 b 1 -a 1 Z -1 y(n-1) Z -1 b 2 -a 2 Z -1 x(n-2) y(n-2) Z -1 Z -1 Z -1 Z -1 x(n-n) b N -a N y(n-n) 31

32 デジタルフィルタ :IIR (cont d) IIR フィルタの伝達関数 k k H(z) = N = 1 N k 1 b k Z a k Z k b(1)+b(2)z -1 +b(3) Z b(n+1)z -N 1+a(2)Z -1 +a(3) Z a(n+1)z -N 因数分解 ( 次数 :N) 零点 / 極 H( z) H ( z ( z 1)( z )( z 1 2) ( z ) ( z 2 N ) ) N 零点極 32

33 FIR フィルタと IIR フィルタとの比較項目 FIR (Finite Impulse Response) IIR (Infinite Impulse Response) インパルス応答有限無限次数多い少ないフィードバック構造なしあり安定性安定不安定群遅延特性線形非線形 H/W コスト高い低い速度遅い早い 33

34 IIR フィルタの種類と特徴項目バタワースチェビシェフ Ⅰ 型チェビシェフ Ⅱ 型楕円通過域リップルなしありなしあり遮断域リップルなしなしありあり群遅延特性特徴波形イメージリップル特性を持たないため平坦な周波数特性が得られる通過域リップルを大きくすることで急峻なロールオフ特性が得られる遮断域の減衰量を小さく指定することで急峻なロールオフ特性が得られる通過域と遮断域双方のリップルを制御可能なため最も急峻なロールオフ特性が得られる 34

35 振幅 (db) Ex2_1) デジタルフィルタ設計手順関数ベース ( 従来手法 ) 既述がシンプル処理の可読性 Sig 振幅応答 (db) フィルタオブジェクトベース正規化周波数 ( ラジアン / サンプル ) フィルタの属性管理豊富な解析機能固定小数点化仕様ベースアルゴリズムの絞り込み試行錯誤削減 35

36 振幅 (db) 振幅 (db) 群遅延 ( サンプル ) 振幅 (db) 振幅 (db) Ex2_2) 各種フィルタの特性比較 Sig 振幅応答 (db) 振幅応答 (db) バタワースチェビシェフ I 型バタワースチェビシェフ I 型チェビシェフ II 型楕円 45 群遅延正規化周波数 ( ラジアン / サンプル ) バタワース vs チェビシェフ Ⅰ 振幅応答 (db) バタワースチェビシェフ II 型正規化周波数 ( ラジアン / サンプル ) 振幅応答 (db) バタワースチェビシェフ I 型チェビシェフ II 型楕円バタワースチェビシェフ I 型チェビシェフ II 型楕円正規化周波数 ( ラジアン / サンプル ) バタワース vs チェビシェフ Ⅰ/ チェビシェフ Ⅱ/ 楕円 ( 群遅延特性 ) 正規化周波数 ( ラジアン / サンプル ) バタワース vs チェビシェフ Ⅱ 正規化周波数 ( ラジアン / サンプル ) バタワース vs チェビシェフ Ⅰ/ チェビシェフ Ⅱ/ 楕円 ( 振幅特性 ) 36

37 振幅 (db) Ex2_3) フィルタ仕様と手法選択 Sig フィルタ仕様 Ap IIR 型サンプリング周波数: 1k[Hz] フィルタ次数:1 次以下左側遮断周波数(Fst1): 3k[Hz] 右側遮断周波数(Fst2): 9k[Hz] 通過帯域(Fp1, Fp2):5~7k[Hz] 通過帯域リップル(Ap):.5[dB] 遮断帯域減衰量(Ast1, Ast2): 4[dB] Ast1 Fst1 Fp1 Fp2 Fst2 Ast2 frequency[hz] 振幅応答 (db) バタワースチェビシェフ I 型チェビシェフ II 型楕円周波数 (khz) バタワース (12 次 ):NG チェビシェフ Ⅰ Ⅱ 楕円 (8 次 ):OK 37

38 振幅 (db) 振幅 (db) Ex2_3) フィルタ仕様と手法選択 Sig フィルタ仕様 Ap IIR 型サンプリング周波数: 1k[Hz] フィルタ次数:1 次以下左側遮断周波数(Fst1): 3k[Hz] 4k[Hz] 右側遮断周波数(Fst2): 9k[Hz] 8k[Hz] 通過帯域(Fp1, Fp2):5~7k[Hz] 通過帯域リップル(Ap):.5[dB] 遮断帯域減衰量(Ast1, Ast2): 4[dB] Ast1 Fst1 Fp1 Fp2 Fst2 Ast2 frequency[hz] 振幅応答 (db) 振幅応答 (db) バタワースチェビシェフ I 型チェビシェフ II 型楕円バタワースチェビシェフ I 型チェビシェフ II 型楕円周波数 (khz) 周波数 (khz) バタワース (12 次 ):NG チェビシェフ Ⅰ Ⅱ 楕円 (8 次 ):OK バタワース (22 次 ) チェビシェフ Ⅰ Ⅱ(12 次 ) :NG 楕円 (8 次 ):OK 38

39 フィルタ設計アプリ : filterbuilder Filter Builder 環境の起動 39

40 フィルタ設計アプリ : filterbuilder(cont d) フィルタの種類を選択周波数仕様振幅仕様の設定フィルタの方式を選択 4

41 フィルタ設計アプリ : filterbuilder(cont d) 振幅特性位相特性群遅延特性各種情報 ( フィルタの安定性使用する H/W コスト等 ) 41

Ex2_4)フィルタ処理(Ex1_3の続き) 4 Sig ピリオドグラムパワースペクトル密度推定 -2 ピリオドグラムパワースペクトル密度推定 -2-3 -3-4 -4 2 パワー/周波数

周波数 (khz) 35 4 45 2.5 3 周波数 (Hz) 3.5 4 4.5 5 不要成分の除去 x 1-3 x 1 5 5 4.5 4.5-3 -1 4 4 3.5 3.5 3 時間時間 -2 2.

42 Ex2_4)フィルタ処理(Ex1_3の続き) 4 Sig ピリオドグラムパワースペクトル密度推定 -2 ピリオドグラムパワースペクトル密度推定パワー/周波数 (db/hz) パワー/周波数 (db/hz) 周波数 (khz) 周波数 (khz) 周波数 (Hz) 不要成分の除去 x 1-3 x 時間時間時間 x 周波数 (Hz) x x

43 Voltage [mv] Voltage [mv] Voltage [mv] Voltage [mv] Voltage [mv] Ex2_5) 適応フィルタの使用例 Sig/DSP 4 Maternal Heartbeat Signal 4 Measured Signal MHB( 母親 ) 3 2 f( ) 回 / 分 1 1 測定データ mV Time [sec] Time [sec] FHB( 胎児 ) 139 回 / 分.25mV Fetal Heartbeat Signal Convergence of Adaptive Noise Canceller Measured Signal Error Signal f() を推定し MHB 成分を除去 Time [sec] Time [sec].5 Steady-State Error Signal MHB FHB f(mhz) + FHB + noise 残差 ( 復元された FHB) Time [sec] 43

44 補足 :Phased Array System Toolbox レーダー波形解析レーダ方程式センサーアレイ解析 44

45 Agenda 1. センサーシステム設計の概要 2. センサーデータの解析周波数解析信号の遅延類似性の評価等 3. センサーデータのフィルタリングフィルタの種類と特徴フィルタの設計 4. システムモデリングから実機テスト固定小数点化システムモデルへの展開実装環境 5. まとめ 45

46 振幅 (db) 変数の固定小数点化 Sig/DSP/Fix 固定小数点長を設定振幅応答 (db) -1-2 フィルター #1: 量子化フィルター #1: リファレンス周波数 (khz) 浮動小数点 ( リファレンス ) と固定小数点との比較 46

47 Simulink ブロックの自動生成 Sig/DSP/Fix/SL 高抽象度ブロックを使用プリミティブブロックを使用 47

48 Ex3_1) Simulink 環境への展開 Sig/DSP/Fix/SL filterbuilder により自動生成されたブロックノッキング信号生成時間軸波形周波数軸波形 48

49 I/O, 可視化系機能強化 Sig/DSP MIDI Control 外部コントローラでパラメータ調整 Time Scope Spectrum Analyzer Peak 検出統計値表示平均値 RMS デューティ比 Max/Min ホールド複素数表示 49

50 I/O, 可視化系機能強化 (cont d) Sig/DSP Logic Analyzer System Object アナログおよびデジタル表示フォーマット Bin, Hex, Oct, Dec 表示カーソルを合わせた箇所の値を表示デジタル 2 値表示デジタル (Hex) 表示アナログ表示 5

51 RoTH : Simulink ハードウエア連携機能 RoTH: Run on Target Hardware Simulink 単体でモデル組み込み実装を実現モデル C コード生成 ( 隠蔽化 ) コンパイル / 実装 Student Version でも可能! R213a~ Gumstix Overo Arduino BeagleBoard/ PandaBoard LEGO MINDSTORMS NXT Digilent Analog Discovery Raspberry Pi 51

52 Ex3_2) Raspberry Pi 赤外線センサモデル SL 音声を再生画像データをセーブモーションセンサが動きを検知したタイミング実行 (Simulink 標準デモとして提供 ) 52

プロトタイピング用浮動小数点 CPU 向けコード Embedded Coder 組込み用 Cコード生成

53 補足 : C コード生成検証オプション MATLAB Coder MATLAB コードから等価な C コードを生成スタンドアローン化アルゴリズムの配布 Simulink Coder SimulinkモデルからCコード生成プロトタイピング用浮動小数点 CPU 向けコード Embedded Coder 組込み用 Cコード生成固定小数点コード対応各社 IDEとリンクデバッグの自動化ターゲット用デバイスドライバブロックデスクトップまたは組込み用 C コード生成および検証 53

補足 : HDL コード生成検証オプション Filter Design HDL Coder MATLAB で設計した IIR/FIR/

Z -1 Z -1 Z -1 Z -1 + + + + HDL Coder Stateflow, MATLAB コード Simulink

ベンダツールと連携 HDL Verifier HDLシミュレータと協調シミュレーション ModelSim / Incisive(NC

54 補足 : HDL コード生成検証オプション Filter Design HDL Coder MATLAB で設計した IIR/FIR/ マルチレートフィルタから HDL 生成ターゲット依存しない VHDL/Verilog テストベンチ (HDL, Model) 生成 Z -1 Z -1 Z -1 Z HDL Coder Stateflow, MATLAB コード Simulink モデルから HDL 生成ターゲット依存しない VHDL/Verilog テストベンチ (HDL, Model) 生成 FPGA ベンダツールと連携 HDL Verifier HDLシミュレータと協調シミュレーション ModelSim / Incisive(NC Verilog) FPGA-in-the-Loop SystemC TLM 2. 生成 ASIC/FPGA に実装可能なコードを生成および検証 54

55 補足 : 信号処理トレーニングコース MATLAB による信号処理 Simulink による信号処理 55

56 Agenda 1. センサーシステム設計の概要 2. センサーデータの解析周波数解析信号の遅延類似性の評価等 3. センサーデータのフィルタリングフィルタの種類と特徴フィルタの設計 4. システムモデリングから実機テスト固定小数点化システムモデルへの展開実装環境 5. まとめ 56

57 まとめ 1 センサーデータの解析 FFT をベースとした周波数解析基本環境窓関数スペクトル推定相互相関等豊富な信号処理関数ライブラリセンサーデータの処理与えられた仕様をもとにフィルタ設計多様なフィルタ手法と最適なフィルタ探索 Filterbuilder 環境システムモデリングから実機テスト変数の固定小数点化簡易実装環境汎用性の高い自動コード生成機能 57

58 まとめ 2 Signal Processing Toolbox DSP System Toolbox 信号解析フィルタリングシステム設計 Fixed Point Designer Simulink 実機評価 Sig: DSP: Fix: SL: Signal Processing Toolbox DSP System Toolbox Fixed Point Designer Simulink 58

Signal Processing Toolbox

Signal Processing Toolbox 信号処理解析およびアルゴリズム開発の実行 Signal Processing Toolbox はアナログおよびデジタル信号処理 (DSP) の業界標準アルゴリズムを提供しますこの Toolbox を使用すると時間領域および周波数領域での信号の可視化スペクトル解析における FFT の計算 FIR および IIR フィルターの設計コンボリューション