Comments
Description
Transcript
ディジタル信号処理の基礎 - 東京電機大学 音響信号処理研究室
講習概要 （基礎編） 「ディジタル信号処理の基礎」 １．アナログ信号とディジタル信号 ２．時間領域と周波数領域 2008. 11. 6,7 ３．線形システム 東京電機大学 ４．ディジタルフィルタ 金田 豊 [email protected] 音の伝送 と ディジタル信号処理 講習概要 （発展編） ５．伝達関数による対象系のモデル化 ディジタル アナログ Ａ/Ｄ ６．インパルス応答の測定法 電気 ７．逆フィルタ ８．適応フィルタとその応用 マイク 表示 処理 伝送 Ｄ/Ａ アナログ 圧縮 雑音除去 認識 音 電気 蓄積 スピーカ 例えば、携帯電話 音 ９．受音系における信号処理の応用 １．アナログ信号とディジタル信号 Ａ／Ｄ（アナログ－ディジタル）変換の手順 ① Ａ/Ｄ（アナログ－ディジタル）変換 ② 標本化定理 ③ Ｄ/Ａ（ディジタル－アナログ）変換 ④ 最近のＡ/Ｄ，Ｄ/Ａ 変換 アナログ信号 → 標本化（サンプリング） → 量子化 ⑤ 標本化定理再考 ⑥ ディジタル世界とアナログ世界 → ディジタル信号 1 標本化（サンプリング） アナログ信号（連続信号） 一定の周期 Ｔｓ で アナログ信号ｘ(t) の値を求めること ｘ(t) Ｔｓ： 標本化周期 （ｆｓ＝1/Ts： 標本化周波数） Ｔｓ ｔ 0 ● ｘ(t) 時間 01 Ｔｓ 2Ｔｓ 3Ｔｓ ・ 時間とともに値が連続的に変化 標本化した信号 ＝ 離散（時間) 信号 離散（時間）信号 ｘ(k) ｘ(2) １ 離散信号からディジタル信号 ｘ(5) ● ● ｘ(1) 0 ３ ２ ｔ ● ・ どの時刻においても値を持つ ｔ：実数 Ｔｓ 4Ｔｓ ● ４ k ５ 離散信号 ＝ 実数の数列 ｛x(0), x(1), x(2), ……｝ ● ● ｘ(3) ｘ(4) ｘ(0) ・ 時間が整数値 （Ｔｓを単位) → 整数 ｋで表す。 ・ アナログ時間との対応は ｔ＝k･Ｔｓ 「ｎ」など も使う ・ 振幅は実数値であるので 離散信号は実数の数列 量子化 ディジタル信号 ＝ 整数の数列 ｛x(0) ’, x(1) ’, x(2) ’, ……｝ まとめて ｛x(0), x(1), x(2), x(3), ……｝ ｘ(k) と表す 例）｛-5.9, 0.1, 3.29, -2.333, ……｝ デジタル信号の実用的データ形式 （離散信号）実数値振幅 → 整数値振幅（ディジタル信号） 量子化単位 Δ の整数倍の値 ４ 例） （ コンピュータ や ＣＤ の内部で、ディジタルデータは、 Ｎ個の１と０の組み合わせで表現されている） 量子化単位 Δ 量子化 ● ● ３ ● 離散 ディジタル ２ ● １ 例 ） 語長 Ｎ が 16 （＝16ビット) のデータ ＝ ディジタルオーディオ（ＣＤなどの）データ 0 1 0 0 1 1 10 1 0 1 1 0 1 0 0 ２バイト整数形式 (short int) で ０ 時間 離散 ｘ（ｋ） ： ｛1.41Δ，2.62Δ， 3.3Δ，… ｝ Δで割って ↓ ↓ ↓ 四捨五入 ディジタル ｘ‘（ｋ） ：｛ 1， 3， 3， ……… ｝ - ２１５ ～ (２１５－１) → -３２７６８～＋３２７６７ を表現 その他） ８ビット（電話）、ほか 2 クイズです クイズです 第２問 第１問 アナログ量は ディジタル量は ２進数 で表され、 １０進数 で表される 次の量のうち、ディジタル量とみなしても よいものに○を、みなせないものに × をつけよ。 (1) 200 (2) 0.13 正しい（○） 誤り（×） (3) 3/7 標本化 誤差はない（標本化定理） 量子化 量子化誤差 ディジタル信号 量子化誤差 量子化単位 Δ アナログ信号 （ 量子化ステップ幅 ） ＡＤ変換時の情報損失（誤差) ● ４ ● ３ 離散 ディジタル ● ２ Δ ● １ 時間 ０ 一般的な信号に対する量子化誤差はランダムな雑音 → -Δ/２ ～ Δ/２の間で一様分布 （量子化雑音） 0 1 0 0 1 1 10 1 0 1 1 0 1 0 0 雑音パワー Δ２/１２ -Δ/２ 量子化単位Δ と 量子化雑音の具体例 Δ/２ → 導出は付録 量子化雑音とＳＮ比 例 ） 16-bit 量子化 信号の ＋１０Ｖ 想定最大値 ２15 等分 物理量 （電圧） 例） １６ビットデータで表される、最大振幅正弦波 ｘ（ｋ） ＋３２７６７ ディジタル量 ２16 等分 ０ ２15 等分 信号の 想定最低値 －１０Ｖ Δ ＝１０Ｖ／３２７６８ ≒ ３００μＶ －３２７６８ 量子化雑音の実効値 √ Δ２/１２ ≒ １００μＶ ｘ（ｋ）＝ ２15 Δ sin ωｋ 信号 ｘ（ｋ）のパワーは （振幅の２乗の1/2） Ｐｘ ＝ ２29 Δ2 一方、量子化雑音のパワー Ｐn＝Δ2/１２ ＳＮ比（Ｐｘ /Ｐn）は 約 ９８ ｄＢ 右下に、関連する テキストのページを 表示します。 Ｌ ビット最大振幅正弦波のＳＮ比 ≒ ６L + ２ [ｄB］ （ p.69） 3 ビット数に関する補足 ［ Ａ/Ｄ変換時の注意 １ ］ 過大入力 過大入力 （オーバーフロー） ・ １６ビットでＡＤしても、計算機内部での演算は ２４ビット，３２ビット、と、ビット数を増したり また、浮動小数点で行う。 （桁落ちや切り捨てなどの 演算誤差の影響を 軽減するため） → ハードクリップされる ＋最大 －最大 ●処理結果に重大な影響 （よくやるミスの一つ） ⇒ 対策例： 測定結果をグラフ表示して、視覚で確認 ●PCのA/Dでは、最大値より小さな値でクリップ される場合がある！ （アナログ保護回路） 直流成分の影響の例 ［ Ａ/Ｄ変換時の注意 ２ ］ 直流成分（バイアス） n2 n2 b 0 0 n1 n1 b ● 処理結果に重大な影響 ● 例えば， ・ 信号の２乗和をとってパワーを計算する場合 ・ 信号間の相関を計算する場合、など n1，n2 はお互いに 関係の無い値をとるので 無相関になる しかし、n1，n2 に直流値 ｂ が加わると、無相関に はならない。 n1が大きければ n2 も大きいという関係がついてしまう 信号処理暗黙の前提 バイアス（直流成分）＝ゼロ は、 信号処理の暗黙の前提 であることが多い 離散信号 と ディジタル信号 標本化 標本化 ＋ 量子化 時間方向での 離散化 時間 ＋ 振幅方向 での離散化 例えば、音も直流はゼロ Ａ/Ｄ したデータに信号処理を行う場合、 まず直流成分を除去しておくほうが無難 （平均値を全体から引く、などで） 4 標本化と量子化 離散信号 ディジタル信号 （注意） 通常、 ディジタル信号処理理論で扱うのは、 ディジタル信号ではなく、離散信号 ・ ディジタル信号処理の理論解析は、通常 離散信号で行う。（暗黙の慣行） ｘ(k) 整数では割り算の取り扱いが面倒なので 実数値の離散信号を使用 → 厳密な本は 「離散時間信号処理」 と呼んでいる （ p.2，p.13） ただし、多くの場合、 離散信号≒ディジタル信号 と見なせる 現実的には、 離散信号＝真の信号値＋測定誤差（電気的雑音など） である。 一般には、量子化雑音 ≪ 測定誤差であるので、 量子化雑音 16-bit 8-bit 0.5 0 0 -0.5 -0.5 -1 0 5 ディジタル信号 ≒ 離散信号 例外）有限語長演算の問題 （桁落ちなど） フィルタ設計への影響 （極・零位置の誤差） 本講習では、ディジタル信号と離散信号の区別をしない 4-bit 0.5 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 10 2-bit ms. 15 20 -1 0 5 10 15 20 0 5 10 ms. 15 20 1.5-bit zero-cross 0.2 （３値） 0.1 0 -0.1 0 5 10 ms. 15 20 -0.2 Ａ/Ｄ 変換 （まとめ） ｘ(t) ｘ(0) ｘ(1) Ａ／Ｄ 変換器 ｘ(2) ｘ(3) 標本化＋量子化 ・ ・ ◇ 連続信号を入力すると（整数）数列を出力 ◇ 信号処理を考えるときは実数値数列（離散信号）を 用いる ＝量子化雑音を無視 ☆ 過大入力と直流成分には気をつけよう 5 標本化周波数 １．アナログ信号とディジタル信号 ① Ａ/Ｄ（アナログ－ディジタル）変換 ｘ(t) Ｔｓ ② 標本化定理 ｔ 0 1 ③ Ｄ/Ａ（ディジタル－アナログ）変換 ④ 最近のＡ/Ｄ，Ｄ/Ａ 変換 Ｔｓ： 標本化周期 （または、サンプリング周期） ⑤ 標本化定理再考 1/Ｔs ＝ ｆｓ： 標本化周波数 （または、サンプリング周波数） ⑥ ディジタル世界とアナログ世界 １秒間のデータ数を表す 標本化定理 （サンプリング定理） 信号の帯域幅 標本化定理 （周期で表せば） ｆmax fmax ＜ fs/2 0 Hz からの場合、信号に含まれる最大周波数 標本化周波数 fs 1/Tmax ＜ 1/（Ｔs・2） fs ＞ 2･ｆmax または Tmax Ｔs 2＜ fmax ＜ fs/2 この条件を満たせば、（理論的には）元の情報は 失わない。（ ＝ 原信号を再現できる） 周期信号の１ 周期に２回以上、 標本化すれば、元の情報は失わない。 標本化周波数の例 信号の パワースペクトル 信号が、どのような周波数成分を どの程度含んでいるか、を表す図 ｆmax ＜ fs/2 を満足 60 例） 60 「あ」 50 50 40 40 30 30 振幅 [dB] ２）オーディオ 信号周波数の上限 ｆmax ＝ 20 kHz 標本化周波数 fs ＝ 44.1kHz，48kHz 振幅 [dB] ｆmax ＜ fs/2 を満足 パワー （ 成分の大きさ） １）電話 信号周波数の上限 ｆmax ＝ 3.4 kHz 標本化周波数 fs ＝ ８kHz 20 10 20 10 0 0 -10 -10 -20 0 「い」 500 1000 1500 2000 周波数 [Hz] 2500 周波数 ｆ 3000 3500 4000 -20 0 500 1000 1500 2000 周波数 [Hz] 2500 3000 3500 4000 (p.10, 23) 6 標本化定理が満たされないと 折り返しひずみ の例 折り返し歪み（エリアシング）が発生 （３/４）fｓ → （１/４）ｆｓ パワー パワー パワー 信号の パワー スペクトル 1 ―ｆｓ 2 0 周波数 f 0 パワー f 0 ｆ 1 ３ ―ｆｓ ―ｆｓ 2 ４ 0 fｓ → ０ １ ―ｆｓ 2 1 1 ―fs 以上の周波数が含まれると、０～ ― fs の区間に折り返される 2 2 0 1 ―ｆｓ 2 ｆｓ ｆ 1 ― ｆｓ 2 0 ｆ ｆｓ 不自然な雑音になる （ p.7） 折り返しひずみの発生（１） 周波数が （3/4）ｆｓ の信号 ｆｓ： 標本化周波数 標本化周期 Ｔｓ ｆ 1 1 ３ ―ｆｓ ―ｆｓ ―ｆｓ 4 2 ４ （３/４）ｆｓ 標本化 ｔ ｔ 1 区別がつかない （１/４）ｆｓ ｔ 1 信号の周期 Ｔ 周期Ｔが 周波数が (3/4)ｆｓ 1 4 Ｔs ＝ (3/4)ｆｓ 3 3 Ｔ Ｔs＝ 4 折り返しひずみの発生（２） ３ 1 ― ｆｓの信号と ― ｆｓの記号が同じディジタル信号になってしまう ４ ４ 折り返しひずみ アナログ世界から見た折り返しひずみ 周波数が ｆｓ の信号 1 1 ｔ ｔ 周波数が ｆｓ の信号は、直流と同じディジタル信号になってしまう Ｄ/Ａ ＬＰＦ Ｄ/Ａ ＬＰF 1 同じ 1 ― ｆｓ ４ 区別がつかない ｔ Ａ/Ｄ ｔ 1 （直流） 1 ― ｆｓ ４ ３ ― ｆｓ ４ （ ｆｓ ） Ａ/Ｄ 同じ 1 1 ― ｆｓ ４ ｆ ＞（１/２）ｆｓ の信号が、Ａ／Ｄされ、Ｄ／Ａされると （１/２）ｆｓ で折りかえった信号に変化する 折り返しひずみ 7 折り返し防止フィルタ 入力 折り返し防止 （低域通過）フィルタ 現実的な折り返し防止フィルタ ゲ 1 イ ン 折り返しひずみを防ぐためには、 Ａ/Ｄ変換器の前に、 （ｆｓ/2） 以上の成分を除去する 「折り返し防止フィルタ」を設置 ゲ 1 イ ン 周波数 0 遮断周波数は ｆｓ/２より小さい 例) ｆｓ＝8kHz 遮断周波数＝3.4kHz Ａ/Ｄ ｆｓ/２ ｆｓ/２ 急峻なフィルタは 時間応答が長い 最近は、AD変換器に 一体化されている場合が多い 内部処理による折り返し歪 発生の例 折り返し歪発生原因 ◇ 入口 周波数 0 スイープ正弦波（ディジタル信号） 1 ◇ 内部処理 ！ 120 100 0.5 周波数 ＡＤ時の帯域オーバー → 最近は、AD変換器に一体化されて いる場合が多いので、考えなくて良い （というより、選択の余地なし → 特性の良いものを選ぶ必要 ） 0 -0.5 80 60 40 20 -1 0 10 20 30 40 50 60 70 80 90 100 20 時間 波頭をディジタル処理でクリップ 60 80 100 120 時間 折り返しひずみ 120 100 周波数 0.5 非線形処理による高調波成分の発生 → 後述 （ディジタル信号の掛け算、 波形の非線形加工、ほか） 40 高調波ひずみ 1 0 -0.5 80 60 40 20 -1 0 10 20 30 40 50 時間 60 70 80 90 100 20 40 60 80 100 120 時間 標本化定理 と 折り返しひずみ （まとめ） fs ＞ ２・fmax が、A/D→D/A で原信号を復元するための 必要十分条件（理論的には）。 これ満たされないと 折り返しひずみが発生する ２つのアナログ信号が 同じディジタル信号に なってしまう現象 8 Ｄ/Ａ変換の理論的しくみ １．アナログ信号とディジタル信号 ① Ａ/Ｄ（アナログ－ディジタル）変換 ディジタル信号 （数字の列） ② 標本化定理 ③ Ｄ/Ａ（ディジタル－アナログ）変換 アナログ 信号 パルス列 （波形） ｛x(0),x(1),x(2)……｝ ④ 最近のＡ/Ｄ，Ｄ/Ａ 変換 理想 ローパス フィルタ 理想的 Ｄ/Ａ ⑤ 標本化定理再考 ⑥ ディジタル世界とアナログ世界 パルスを なめらかに つなぐ Ｄ/Ａ （広義) ｓｉｎｃ 関数 （予備知識） どうやって 「なめらかに」 つなぐ？ ｙ（ｘ）＝ ｓin ｘ ｘ ＝ ｙ 1 ｘ ｓin ｘ （ 振幅が 1/x の正弦波） 1 ｘ ｓin ｘ 元のアナログ波形を 復元するためには ？ （予備知識） π 2π 3π ｘ＝0 では ? ｓｉｎｃ 関数 ｘ ｙ（ｘ）＝ ｓin ｘ ピークの 包絡 波形の復元（1） 1 ｘ ｘ（0） ｙ １ ｘ 0 ｘ（1） ｔ＝0 以外のサンプル点 で、0 となる ｘ＝π，２π，３π，．．． で、0 となる ｔ 原点で高さは１ ｘ（0）の高さを持ち、 0 ｘ π 2π 3π 標本化周期の２倍の周 期を持ったｓｉｎｃ 関数 ｘ（2） ｘ（3） → ｘ＝0 の値は付録 9 波形の復元（３） 波形の復元（２） ｘ(1) ｘ(0) ｘ（0） ｘ（1） ｘ（1）の高さを持った ｔ ｓｉｎｃ 関数 ｔ ｘ（k），k＝0,1,2,... の高さの ｘ（3） ｘ（2） ｓｉｎｃ 関数を重ね足し合わせ ると、原波形が復元される。 ｘ(3) ｘ(2) ｓｉｎｃ 関数による 「補間 （内挿）」 理想ローパスフィルタ (p.53) 波形復元の数学的表現 ゲイン パルス列波形 ｓｉｎｃ関数 １ 0 0 fｓ/２ ◇ 理想ローパスフィルタとは、 0 ～ (ｆｓ/２) の周波数 → ゲイン １ (ｆｓ/２) 以上の周波数 → ゲイン 0 t （通過） （阻止） ◇ パルス列を理想ローパスフィルタに通せば、 ｓｉｎｃ 関数で補間される 現実のＤ/Ａ変換器の問題点 パルス列 （波形） Ｄ/Ａ 理想パルスが 出力できない ・ アパーチャ効果 ( 保持効果 ） (p.54) たたみ込み ｘ(0) ｘ(1) ｆ アナログ信号 理想 ローパス フィルタ 1 理想ローパス が存在しない ・ ｆs/2 付近の特性の乱れ ＊ t ｘ(2) ｘ(3) 波形の復元 （補間） は、 パルス列波形と ｓｉｎｃ関数との たたみ込み （＝ 理想ローパスフィルタ に通す ） 最近の変換器では ◇ 最近の変換器 （オーバーサンプリング、∑Δ） では特性が改善された。 （理想ローパスフィルタをディジタルフィルタ として近似実現したため） ＤＡとフィルタは一体化されたので、 フィルタを意識しなくても良くない場合が多い。 ・ 折り返し誤差 10 しかし、使用するDAは、一度チェックを しておいたほうが better 例） ＤＡ （低域フィルタ を含む） ディジタル 正弦波 オシロ レベルメータ AD など ◇ 特に ｆs/2 付近の周波数特性 ◇ 特に、パソコン（オーディオ用）の変換器 ◇ sinc 関数によるプリエコー発生の問題 （聴覚への影響） （計測用・ 高級品でも） ◇ 効果的なチェック方法については、現在、検討中 近似理想ローパスフィルタの問題 パルス列 （波形） アナログ信号 近似理想 ローパス フィルタ パルス 近似理想 ローパス フィルタ 1 定常信号に 対しては OK 非定常信号 （過渡的・パルス的） に対しては？ プリエコーが聞こえる 場合がある、らしい 右側応答のみを利用したフィルタもある、らしい ＤＡ変換 （まとめ） 理論的には、 ディジタル信号 （数値列） 理想 ＤＡ ｓｉｎｃ関数の 理想 たたみ込み ＬＰＦ アナログ信号 パルス列 の再現 現実的には、 ・ 古い型では、ｆs/2 付近の特性の乱れ。 ・ 最近は、問題が少なくなったようであるが、 特性をチェックしてから利用することが望ましい。 11 最近のＡ/Ｄ，Ｄ/Ａ 変換 １．アナログ信号とディジタル信号 ① Ａ/Ｄ（アナログ－ディジタル）変換 １） オーバサンプリング方式 ② 標本化定理 ２） Σ⊿ 方式 （＝ ⊿ Σ方式 ） （１ビット Ａ/Ｄ，Ｄ/Ａ） ③ Ｄ/Ａ（ディジタル－アナログ）変換 ④ 最近のＡ/Ｄ，Ｄ/Ａ 変換 ⑤ 標本化定理再考 ⑥ ディジタル世界とアナログ世界 １） オーバサンプリング方式 オーバサンプリングの効果 信号が含む上限周波数を ｆm、 標本化周波数 ｆs が大きくなると、 ディジタル化できる周波数帯域 ｆs/2 が増加。 標本化周波数 を ｆs とする。 ◇ 基本方式： ｆs ＝ 2･ｆm ◇ オーバサンプリング方式： ｆs ＝ 2･ｆm ×Ｎ （Ｎ：2以上の整数） パワー 目的帯域 0 サンプリング定理で必要とさせる周波数を 上回る周波数でサンプリング 目的帯域 パワー 0 基本方式の 量子化雑音 ｆm N倍 ｆs/2 N･ｆm ｆ (p.57) 利点 その１ 目的帯域内の量子化雑音が減少 (ｆs/2) ｆm 利点 その２ 量子化雑音は白色雑 音で、そのパワーは ディジタル化した 帯域内（0～ ｆs/2） に一様に分布するの で、標本化周波数 ｆsを 増加すれば、目的 帯域内に含まれる量 子化雑音は減少する。 N･ｆm ｆ オーバサンプリングの 量子化雑音 (p.71) Ａ/Ｄ，Ｄ/Ａ用のアナログ低域フィルタは低次のものでよい ｆm 以上の成分は、 ディジタル化後、 ディジタルフィルタで 除去 パワー 基本方式： fm と fs/2 が近 いので、急峻な フィルタが必要 ｆm N･ｆm ｆ オーバサンプリング方式： 緩やかな傾斜の 低次のフィルタでＯＫ 12 利点 その３ 基本方式との受け渡し Ｄ/Ａ 時 のフィルタの悪影響が改善 標本化周波数が大きくなったので、下図に示すように、 ローバスフィルタの遮断周波数 ｆc 付近の悪影響 （遅延、リップル、アパーチャ効果など）の影響が小さい ゲ イ ン 1 信号成分 0 ｆｃ ｆｓ/２ 周 波 数 オーバサンプリ ング方式 基本方式 間引き １秒間に N･ｆs0 個の １秒間に ｆs0 個の データ データ 補間 CD規格データなど 例えば、fs0＝44.1kHz ○ オーバサンプリング周波数の例 20kまでの信号 → ｆｓ=48k → 192kHz (p.79) 量子化の表現 ２） Σ⊿方式 （１ビット方式） １） ＭＨｚ 以上もの 高速オーバサンプリング ２） 微分効果で量子化雑音を大幅に低減 ３） Ａ/Ｄ （量子化器）は 簡単な １ビット 高速 アナログ ＋ － 遅延 （⊿Σ方式とも呼ぶ） ＋ － ＤＡ 等価回路 離 散 化 積分 量子化雑音が微分される 1bit ディジタル AD （量子化） 積分 3.43＋（-0.43） → 3 3.88＋（0.12） → 4 (p.78～) 一つの解釈 アナログ （四捨五入） （量子化雑音が加算されて整数になった、 と考える） （量子化） ＤＡ 微分効果 ⊿ 3.43 → 3 3.88 → 4 ディジタル 1bit AD 積分 Σ 例） （信号の流れ） ∫Ｓ Ｓ Ｓ ＋ 積分 ＋ － ∫Ｓ･z-1 遅延 積分 Ｓ 元と同じ Ｓ･z-1 遅延 微分効果 ＋ ＋ 量子化雑音 ＋ － 積分 遅延 微分効果 （雑音の流れ） Ｎ －Ｎ･z-1 ＋ ＋ Ｎ’ 「’」は微分 を表す Ｎ’ 微分波形 － Ｎ･z-1 量子化雑音 積分 Ｎ’･z-1 遅延 Ｎ’ 13 微分の効果 ・ 微分（差分）は高周波域強調 （低周波域抑圧） d sin ω t = ω cos ω t 微分結果はωに比例 dt d2 sin ω t = ω 2 (− sin ω t ) ω2に比例 dt 2 量子化雑音微分の有効性 目的帯域内の量子化雑音が大幅に減少 １ビットADでOK 量子化雑音が微分さ れた結果、雑音のパワ ーは、帯域内（0～ ｆs/2） の高周波域に集中する ので、目的帯域内の 量子化雑音は減少 する。 パワー ｆｓ はすごく 大きい値 ｆ ｆｓ/２ ノイズシェーピング オーバサンプリングの 量子化雑音 Σ ⊿方式の画期的な点 Σ⊿方式の 量子化雑音 (p.80, 85) １．アナログ信号とディジタル信号 ◇変換器が１ビットでできるため、 高精度・低雑音のＡＤ，ＤＡを、 安く・小さく実現できる（実用的観点） ◇ 振幅方向ではなく、時間方向での分解能向上 （＝サンプリング周波数向上）で 量子化雑音が低減できる （理論的観点） ① Ａ/Ｄ（アナログ－ディジタル）変換 ※ ただし、Ａ/Ｄ 変換器の出力やデータ保存の形式は ⑥ ディジタル世界とアナログ世界 ② 標本化定理 ③ Ｄ/Ａ（ディジタル－アナログ）変換 ④ 最近のＡ/Ｄ，Ｄ/Ａ 変換 ⑤ 標本化定理再考 １６ビット（現在の標準的データフォーマット） に変換されて出力・保存される場合が多い。 ディジタル信号と視覚化 ディジタル信号は数列 ｛x(0), x(1), x(2), x(3), x(4), ……｝ 元のアナログ信号は どんな信号？ 通常のPC上の描画プロットは、ディジタル信号（数列）を直線で結ぶ ｘ(k) ｘ(3) ｘ(0) 例えば、 ｛0.85, -0.34, -0.34, 0.85, -0.98, 0.65, ……｝ である。 しかし、これでは、直観的にわかりにくいので、 グラフ化 （アナログ化）することが多い。 ｘ(5) k ｘ(1) ｘ(2) ｘ(4) 14 元のアナログ信号は正弦波 1.5 標本化定理再考（その１） ◇ 標本化周波数 ｆｓ は大きい程良い？ 1 とりあえず NO： → さきほどの例は、グラフの描き方が不適切 0.5 0 グラフ化とは、時間離散のディジタル信号を、 擬似的なアナログ信号として表示して、 視覚的な理解を得ようとするもの -0.5 -1 -1.5 1 周期に １８ 回の標本化 1 周期に ２.６ 回の標本化 （＞２回） 補間 ｘ(0) ｘ(3) ｘ(1) ｘ(2) ｘ(5) 補間 正確な表示のためには、表示の前に 擬似的なアナログ化（DA）のプロセスが必要 → 補間 Ｄ/Ａ変換の理論的しくみ （再掲） ディジタル信号 （数字の列） アナログ 信号 パルス列 （波形） ｘ(4) ・ 離散時間の間の値を計算により求め、 標本化周波数を上げて、アナログ信号に近づける ・ 「内挿」 「Interpolation」 「オーバーサンプリング」とも言う ｛x(0),x(1),x(2)……｝ 理想的 Ｄ/Ａ ・ 標本化周波数の増加率に応じて、 N倍の補間（オーバーサンプリング）、と言う 補間（内挿）の方法 ◇ オーバサンプルしたパルスとｓｉｎｃ関数をたたみ込む （理想LPFによる D→A の補間のシミュレーション） 例） Ｋ倍の補間 のための ｓｉｎｃ関数 ｆ (k)＝ ｓｉｎ （πｋ/Ｋ）／ （πｋ/Ｋ） ｘ(3) sinc 関数との たたみ込み 標本化定理 再考（その２） ◇ 補間によって 全ての離散信号は元のアナログ信号に 正確に戻るのか？ → 補間を行うための注意事項 （例：Ｋ＝２） ｘ(0) 理想 ローパス フィルタ ｘ(5) ＊ ｘ(1) ｘ(2) ｘ(4) サンプルの間にK-1個の零 15 補間（10倍）をしてみた ｆs/2 近くの正弦波のプロット サンプル値(○)を結ぶと、「うなり」のように見える。 若干の補間効果はあるが、「うなり波」のまま f =3821 H z の 正 弦 波 の プ ロ ット , Fs=8000H z 1 1 0 .8 0 .8 0 .6 0 .6 0 .4 0 .4 0 .2 0 .2 0 0 - 0 .2 - 0 .2 - 0 .4 - 0 .4 - 0 .6 - 0 .6 - 0 .8 - 0 .8 -1 補間にも限界があって、 周波数の高い信号には、 より高い標本化周波数が 必要なのか？ -1 0 5 10 15 20 25 30 [s am p le s ] 35 40 45 50 100 問題は補間の品質 補間には品質があるので注意 1 松 0.8 (MATLAB) 1 0 倍 で 補 間 z z = r e s a m ple (y y ,1 0 ,1 ) 無限長のsinc関数を どこまで利用するか？ による差 0.2 0.3 0 -100 100 0 0 π ｓｉｎｃ 関数を長く利用すると-0.1 -0.2 計算時間が大変になる → 組み込み関数は、 -30 安・近・短 （デフォルト） 梅 200 130 135 140 145 150 ◇ 組み込み関数よりは自作のほうが 安心な場合もある （ ⇔ 欠点は速度） ※ sinc 関数に窓（Hanning など）をかけてから 補間するのも有効 -20 -10 0 π 10 20 30 少し品質の高い補間をしてみた （竹） まだ「うなり波」のように見える 念入りな補間 （松） ほぼ一定振幅の正弦波に「見える」ようになった。 s in c 関 数 を 片 側 1 0 周 期 に 。 1 0 倍 で 補 間 z z = r e s a m ple (y y ,1 0 ,1 ,2 0 ) s in c 関 数 を 片 側 1 0 0 周 期 に 。 1 0 倍 で 補 間 z z = r e s a m ple (y y ,1 0 ,1 ,2 0 0 ) 1 1 0 .8 0 .8 0 .6 0 .6 0 .4 0 .4 0 .2 0 .2 0 0 - 0 .2 - 0 .2 - 0 .4 - 0 .4 - 0 .6 - 0 .6 - 0 .8 - 0 .8 -1 100 125 ◇ 最低限、ヘルプを見て、パラメータを 確認しておくことが望ましい。 （私の反省です） 0.1 -0.4 -200 120 ◇ 理想的には、内容を理解して使う 竹 0.2 -0.2 115 ◇ 安易に組み込み関数を利用して良しとしない 0.5 0.4 110 組み込み関数の利用に当たって 0.6 0.4 105 -1 105 110 115 120 125 130 135 140 145 150 100 105 110 115 120 125 130 135 140 145 150 16 標本化定理 再考（その３） 標本化定理 再考（その４） ◇ 毎回「松」では計算が大変そうだが、 → 大事なところでは念入りに補間すべき ◇ 勘違いや補間の手間を避けたい場合などには いっそ高い周波数で標本化したほうが 楽な場合もありそう。 ◇ 標本化定理の言う 「元のアナログ信号を、 ディジタルで近似的に復元」 するためには、 補間が必要 ◇ 補間をするには、 左右にある程度のデータ数 （sinc 関数長程度）が、必要。 ⇒ 「少数のデータ」や、「信号の端付近」は、補間不良 ⇒ 「少数のデータ」や、「信号の端付近」には、 標本化定理は成立しないので、 これらを正確に表示するためには、 標本化周波数を、標本化定理より上げる必要がある。 標本化定理 再考（まとめ） 注意が必要な例： 最大値の時刻の検出 ◇ サンプルを結んでプロットした時に不自然な 波形であっても、標本化定理を満足していれば 補間によって元のアナログ信号に近づけてプロット することができる。 ◇ しかし、補間は、十分な長さの sinc関数を使わないと 結果が不正確なことがある。 10 5 0 ここが最大値 の時刻だと 判定 -5 ◇ 補間がめんどうな場合は、 標本化周波数を上昇させるのも選択肢 -10 ◇ 「少数のデータ」や、「信号の端付近」では、 補間が行えないので標本化定理は成立せず、 標本化周波数は適切に上昇させる必要がある。 波形の高さを比較する場合には補間をする 15 15 -15 2 4 6 8 時間 10 12 [samples] 14 16 18 20 正確なピーク値、ピーク時刻 正確なピーク値とピーク時刻を求めるには 補間が必須 ピーク値 真の最大値 10 補間あり 5 正解は こちら 0 補間しないと 誤差が発生 -5 補間なし -10 -15 2 4 6 8 時間 10 12 [samples] 14 16 18 20 ピークを与える時刻 17 １．アナログ信号とディジタル信号 ディジタル信号処理 が扱う二つの世界 ① Ａ/Ｄ（アナログ－ディジタル）変換 実世界 (Ｒｅａｌ Ｗｏｒｌｄ） ② 標本化定理 アナログの世界 ③ Ｄ/Ａ（ディジタル－アナログ）変換 ④ 最近のＡ/Ｄ，Ｄ/Ａ 変換 ⑤ 標本化定理再考 ⑥ ディジタル世界とアナログ世界 D/A ディジタル から アナログ世界 を操作する A/D 計算機の中の世界 ディジタルの世界 ディジタル を通して、 アナログ世界 を見る、 理解する 注意しなければいけない事柄 アナログ世界とディジタル世界 （まとめ） ◇ 同じ信号が違って見える。 ・ ディジタルをそのままプロットしたのでは、 元のアナログ信号には見えない ⇒ アナログに近づけた（補間した）表示が、 （場合によっては） 必要となる。 実世界 (Ｒｅａｌ Ｗｏｒｌｄ） アナログの世界 ◇ 演算操作、分析手法などにも注意 ・ 違っているようで同じもの (同じ操作) ・ 似ているようで違うもの (違う操作) を理解しておく。 ＡＤ ＤＡ ディジタルの世界 処理 操作する、制御・応答 見る、計測・分析 ディジタル信号処理 18 フーリエ変換 ２．時間領域と周波数領域 時間信号 ① フーリエ変換 ・ 時間と周波数を関係付ける 主要な性質を復習 （証明抜） ・ 主としてアナログの話 → ディジタルも類似なので 周波数スペクトル フーリエ変換 ∞ X（ｆ）＝∫ｘ(t) ｅ ｘ(t) - j 2πｆ t -∞ 注） ｅ - j 2πｆ t （＝ dt ｅ -jωt） ＝ cos(2πｆ t)＋ ｊ･sin(2πｆ t) ② 関数の直交変換 （講習省略） は 「正弦波」 （ 複素正弦波 ） ③ ディジタル領域の周波数変換 ④ 窓関数 フーリエ変換結果の図示に対して 時間信号とスペクトル （時間信号波形） （周波数スペクトル） ｘ(t) 振 幅 ∞ X（ｆ）＝∫ｘ(t) ｅ t 周波数 ｆ フーリエ変換結果 Ｘ(f) は、 視覚的にわかりやすい |Ｘ(f)|2 で図示す ることが多い （特にことわらない） フーリエ変換 補足：ここだけの話 （スペクトルの縦軸の呼び名について） （周波数スペクトル） 周波数 ｆ 重要な場合もあり 物理量に厳密な先生から お叱りを受けそうですが・・・ 私の屁理屈 信号のスペクトルの相対的な形に情報が含まれており、 縦軸（物理量）は何でも良い。 ｘ(t) 振幅 強度 |Ｘ(ｆ)|2 ◇ 周波数スペクトルの 縦軸の量名称については 世の中でも統一がとれて いません。 強度、振幅、パワー、エ ネルギ、エネルギ密度、 etc. etc. ◇ 私もかなりいいかげん です。（信号の場合は「パ ワー」、フィルタは「ゲイ ン」とする場合が多い） dt ◇ フーリエ変換結果 Ｘ(f) は複素数。 ◇ 各周波数成分の振幅（大きさ）は |Ｘ(f)| （ 強度は |Ｘ(f)|2 ） ◇ 各周波数成分の位相は、偏角 ａｒｇ（Ｘ(f)） Ｘ(ｆ) 時間 - j 2πｆ t -∞ t 時間 もちろん、物理量が意味 を持つ場合もあります。 例えば、左図の音声波形 の縦軸を何とするか？ よく 迷うことがある。 物理量とし ては、マイク出力の電圧か、 音圧が妥当だが、どちらを選 んでも、音声信号処理を行う 際に、全く影響は無い。 情報（信号）を表す物理量 が意味を持たないという例 である。 19 時間信号とスペクトル （時間信号） （周波数スペクトル） t 時間信号 |Ｘ(f)|2 パワー ｘ(t) 振 幅 フーリエ変換の定義 周波数スペクトル フーリエ変換 ∞ X（ｆ）＝∫ｘ(t) ｅ ｘ(t) - j 2πｆ t -∞ 時間 周波数 dt ｆ フーリエ変換 対称な性質が 存在 フーリエ逆変換 ∞ ｘ（t）＝∫X(f) ｅ j 2πf t dｆ 注）この項はアナログの話です。 ただし、 類似の性質がディジタルにも成立します。 時間信号 と スペクトル との関係 -∞ 代表的な フーリエ変換対（１） ◇ 時間信号 と スペクトルは、 一つの信号の表裏の関係をなす。 ◇ 時間領域表現（時間信号） と 周波数領域表現（スペクトル） とのうち 場合に応じて、 見やすい（理解しやすい）領域で 信号処理理論が展開される (p.35) X(f) ｘ(t) Ｆ -ｆc 0 t 1 対称性 Ｆ t 1 X(f) 1 0 sinc 関数 方形パルス 例）理想低域フィルタの 周波数スペクトル ｆc ｆ 理想低域フィルタ sinc 関数 0 ゲイン ｓｉｎｃ フーリエ （周波数スペクトル） 変換対 （時間信号） 時間信号とスペクトルの関係把握は重要 補足： 負の周波数 方形 ｆ → 証明は付録 補足の補足： 負の周波数 ・ 複素正弦波で考えれば e j 2π f t → e − j 2π f t 複素平面の単位円上を逆周り -ｆc 0 ｆc ｆ 実時間信号のスペクトルを考える時には、 正の周波数 （0 ～）を考えれば十分であるが、 時間領域との対称性を見るために、 負の周波数に対する Ｘ(f) の値も表している ・ sin で考えれば、周波数が ｆ の正弦波 sin 2πｆ ｔ に対して、 周波数が - f の 正弦波は、 sin 2π(- ｆ ) ｔ （ ＝ - sin 2π ｆ ｔ ） 正負が逆なだけ。 周期（周波数）は同じ 20 フーリエ変換対（２） （時間信号） Ｆ t （時間信号） ｆ 0 （調波構造） Ｆ t 0 厳密には 離散的信号 離散（時間)信号 周期 の例 ｆ 0 実数スペクトル 離散 ｆ 0 離散（線)スペクトル 周期関数 Ｆ ｆ １/Tｓ Ｆ X(f) t （共役）対称 な 複素関数 X(- f)＝ X＊(f) 対称関数 ｘ(t) ＝ｘ(-t) 周期 （周波数スペクトル） ｘ(t) 0 t 離散 Tｓ X(f) 実関数 0 フーリエ変換対（３） 対称 （周波数スペクトル） ｘ(t) 1 実 フーリエ変換対（４） 等パルス列 （時間信号） ｘ(t) 周期スペクトル 等パルス列 （周波数スペクトル） X(f) １/Tｓ Tｓ Ｆ t ｆ 等間隔・等振幅 線スペクトル列 （パルス列） 等間隔・等振幅 パルス列 ☆ パルス間隔は 逆数の関係 演算の関係 （時間信号） たたみ込み 例） 理想ローパスフィルタ 乗算 （周波数スペクトル） 周波数で表すと X(ｆ) ｙ(t)＝∫ｘ(t-u)ｈ(u) du ＝ ｘ(t)＊ｈ(t) たたみ込み 乗算 （窓掛け、振幅変調） ｆ Y(f)＝X(f)・H(f) 乗算 （フィルタリング) 周波数軸上の たたみ込み × Ｆ 時間表すと 対称性 ｙ(ｔ)＝ｘ(ｔ)・ｈ(ｔ) 乗算 ｘ(t) 方形 ｆ 0 ＝ ｆ Ｆ ＊ たたみ 込み t ｓｉｎｃ 理想ローパスフィルタは、ｓｉｎｃ関数のたたみ込み → ｘ（ｔ）がパルス列の時は、補間の操作になる 21 離散化（標本化）した信号のスペクトル 等間隔・等振幅パルス列のたたみ込み 時間表すと ｘ(t) ｆ0 X(ｆ) スペクトル ｐ(t) 離散信号波形 等間隔・等振幅 パルス列 ｆ0 t Ｆ Ｆ Ｆ Ｐ(f) 周波数で表すと 離散信号のスペクトル ＊ X(ｆ) ＝ ＝ t Ts ｆ たたみ込み （波形） パルス列の乗算 t× ＊ ｆ 離散化 ｆ ｆ＝ ｆ ｆ １/Ts＝ｆｓ ｆｓ 周期化 スペクトルの周期化 周期Tsで時間離散化 → ｆs （=1/Ts)でスペクトルは周期化 （波形の周期化） 離散化 時間領域 周波数領域 離散化 周期化 波形 × 等間隔 パルス 離散化によるスペクトルの周期化と 標本化定理 周期化 スペクトル ＊ ｆｓ X(ｆ) （ 周期化 ） 等間隔 パルス （ たたみ込み ） -ｆｓ/2 0 ｆｓ/2 ｆ - ｆｓ/2 0 ｆs で離散化 → ｆs で周期化 離散化 等間隔 パルス （ たたみ込み ） 波形 ＊ スペクトル × ｆ ｆｓ：標本化周波数 ｆｓ X(ｆ) 周期化 ｆｓ/2 （ 周期化 ） 等間隔 パルス -ｆｓ/2 0 ｆｓ/2 ｆ - ｆｓ/2 0 ｆｓ/2 ｆ ｆｓ/２ を越えたＸ(f)の成分はスペクトルが重なる。 → 「折り返される」ように見える → 「折り返しひずみ」 折り返しの解消 折り返しというより，「回り込み」 非対称なスペクトルの場合を考える ｆｓ1 ｆｓ1 - ｆｓ1/2 0 折り返しが生じたら ↓ サンプリング周波数 を増加させる ｆｓ1 → ｆｓ2 ｆ ｆｓ ｆｓ2 - ｆｓ/2 - ｆｓ2/2 0 ｆ 0 ↓ 折り返しの解消 ｆｓ1/2 ｆｓ2 ｆｓ で離散化 ｆｓ2/2 0 ｆｓ/2 ｆ ｆ ｆｓ/2 からはみ出た成分が、-ｆｓ/2から回り込む 22 Ａ/Ｄの前、Ｄ/Ａの後の 理想低域フィルタの働き 時間軸の反転 （時間信号） Ａ/Ｄの前： ｆｓ/２ を越えたＸ(f)の成分を除去 X(ｆ) -ｆｓ/2 0 ｆｓ/2 ｆ （入力） ｘ(t) 1 × -ｆｓ/2 0 ｆｓ/2 ｆ -ｆｓ/2 0 ｆｓ/2 ｆ （連続化＝補間） × ｆ Ｆ X(f) 時間軸の反転 ｘ(-t) Ｄ/Ａの後： 離散信号に含まれる ｆｓ/２ 以上の成分を除去 - ｆｓ/2 0 ｆｓ/2 t 0 （理想低域フィルタ） （離散信号） （周波数スペクトル） t 0 （連続信号） Ｆ X(f)の複素共役 X＊(f) 1 -ｆｓ/2 0 ｆｓ/2 ｆ -ｆｓ/2 フーリエ変換対の まとめ 0 ｆｓ/2 ｆ X（ｆ）＝∫ｘ(t) ｅ - j 2πｆ t dt 演算操作の関係 方形 ｓｉｎｃ たたみこみ 演算 乗算 実 対称 離散化 周期化 離散 周期 連続化 方形窓の 乗算 等パルス列 （sinc 関数の 畳み込み） 等パルス列 軸反転 複素共役 23 ２．時間領域と周波数領域 関数の直交変換とは 関数ｇ（ｔ）を 直交な成分に分解 ① フーリエ変換 ② 関数の直交変換 （講習省略） ③ ディジタル領域の周波数変換 ④ 窓関数 幾何ベクトルの内積 フーリエ変換は 直交変換となっている 直交基底ベクトル ｅ2 ｂ ｅ1 θ ａ ベクトルの長さ 内積 ： （ａ，ｂ）=|ａ| ・ |ｂ| ｃｏｓθ θ=90°→ 直交 → （ａ，ｂ）=０ （ａ，ａ） → θ=０°→ |ａ|２ ・単位ベクトル : 長さが１ |ｅ１| = １ →内積が１ （ｅ１，ｅ１）= |ｅ１|２ = １ ・直交単位ベクトル ｅ１，ｅ２ （ｅ１，ｅ２ ）= ０ 直交した単位ベクトルの集合 数ベクトルと内積 ベクトルの直交分解 ｇ ｇ２ ｅ２ ０ ｅ1 ｇ１ ベクトルｇを直交した軸の成分で表すこと ａ= ａ１ ａ２ ･ ･ ａＮ ｂ= ｇ = ｇ１＋ｇ２ ＝ｇ１･ｅ１ ＋ ｇ２･ｅ２ 例）２次元ベクトル ｇ１はｅ１とｇの内積で求まる （ｅ１，ｇ ）=（ｅ1，ｇ１･ｅ1＋ｇ２･ｅ２） =ｇ１（ｅ１，ｅ１）＋ｇ２ （ｅ１，ｅ２） １ ０ =ｇ１ （正規）直交基底ベクトル ｂ１ ｂ２ ･ ･ ｂＮ a* b 内積 （ａ,ｂ）= ａ１*ｂ１ + ａ２ *ｂ２ + ･･･ Ｎ =Σ ａi* bi = a*ｂ i=1 *：共役転置 例） 幾何ベクトルの成分表示 → 数ベクトル 24 直交基底関数 関数の直交分解 ｇ2（ｔ） ｅ2（ｔ） 0 ｇ（ｔ） ・ 関数の内積は、「積の積分」で定義される （ｅｉ（ｔ）,ｅｊ（ｔ） ）= ｇ1（ｔ） ｅ1（ｔ） ∫ｅｉ*（ｔ）ｅｊ（ｔ）ｄｔ ・ 直交する *：共役 内積がゼロ 連続関数は無限次元のベクトル ・ 直交基底関数系とは…？ 任意の関数ｇ(t) は、 直交基底関数 ｅ１（ｔ）,ｅ２（ｔ）,…により 内積 = 積分がゼロとなる関数の集まり （ｅ1(t), ｅ２(t), ｅ３(t),……） 直交分解することができる。 ｇ（ｔ） = ｇ１・ｅ１（ｔ）＋ｇ２・ｅ（ｔ）＋ ･･･ （予備知識） 複素指数関数（１） （予備知識） 複素指数関数（２） 複素指数関数 ｅ ｊ2πｆ ｔ ＝ ｃｏｓ(2πｆ ｔ) ＋ ｊ ｓｉｎ(2πｆ ｔ) 複素数 実数部 虚数部 Imag 複素平面 長さ１ ｅ ｊ2πｆ Imag j ｅ ｊ2πｆ ｔ 0 単位円上を回転 Real 2πｆ ｔ -1 ｔ t が増加すると ｅ ｊ2πｆ ｔは、 1 ｓｉｎ(2πｆ ｔ) 2πｆ ｔ ｅ ｊ2πｆ ｔは、 ｃｏｓ(2πｆ ｔ) 長さ １ 偏角 2πｆ ｔ 0 Real 直交基底関数の例 （予備知識） 複素指数関数（３） ◇ 負の周波数 –ｆ を持つ複素指数関数 → 単位円上を逆回転 ｅ -ｊ2πｆ ◇ 正（ｅ ｊ2πｆ ｔ ）と 負（ｅ -ｊ2πｆ ｔ ）とがあれば、 ｓｉｎ、ｃｏｓが復元できる。 ｅ ｊ2πｆ ｔ ＝ ｃｏｓ(2πｆ ｔ) ＋ ｊ ｓｉｎ(2πｆ ｔ) ｅ -ｊ2πｆ ｔ ＝ ｃｏｓ(2πｆ ｔ) － ｊ ｓｉｎ(2πｆ ｔ) t が 1/f で一回転 → 周期 1/f -j ｔ 例）複素指数関数 ｅ ｊ2πｆ ｔ の集まり 内積 ＝∫ｅｉ*（ｔ）ｅｊ（ｔ）ｄｔ １ Ｔ ｰ j2πｆｉ ｔ・ｅ ｊ2πｆｊ ｔ ｄｔ ２Ｔ ｰＴ ｅ ∫ １ Ｔ ＝ ２Ｔ∫ ｅ ｊ2π(ｆｊ-fi) ｔ ｄｔ ＝ ｰＴ ｃｏｓ(2πｆ ｔ) ＝（ｅ ｊ2πｆ ｔ ＋ ｅ -ｊ2πｆ ｔ）/２ ｓｉｎ(2πｆ ｔ) ＝（ｅ ｊ2πｆ ｔ － ｅ -ｊ2πｆ ｔ）/２ｊ ◇ 複素正弦波関数とも呼ばれる 1 ｆi ＝ ｆｊ 0 ｆi ≠ ｆｊ 周波数の異なる複素指数関数は直交 複素指数関数の集まりは直交基底関数系 25 直交分解とフーリエ変換 ・関数ｇ(t) の直交分解は、直交基底関数とｇ(t) の 内積（＝ 積の積分）で得られる。 （ｅｉ（ｔ）, ｇ(ｔ） ）= ∫ｅｉ*（ｔ）ｇ（ｔ）ｄｔ 例） 直交基底関数が、複素指数関数 ｅ ｊ2πｆ ｔ （ｅ ｊ2πｆ ｔ, ｇ(ｔ） ）= ∫ｇ（ｔ）ｅ ベクトル合成とフーリエ逆変換 ベクトルの合成 ｇ = ｇ1･ｅ1 + ｇ2･ｅ2 + ｇ3･ｅ3 + ･･･ =Σｇｉ･ｅｉ 大きさ 基底 -j2πｆ ｔｄｔ ＝ G(f) ｇ（ｔ）=∫Ｇ(f) ｅ j2πｆ ｔ ｄｆ フーリエ変換は複素指数関数を基底ベクトル とした関数の直交分解（直交変換) フーリエ変換の直交変換イメージ ｇ（ｔ） Ｇ（ｆ2） ｅ j2πｆ2ｔ ０ ｅ j2πｆ1ｔ Ｇ（ｆ1） ｇ（ｔ） を基底ｅ j2πｆｔ の成分 Ｇ（ｆ） に分解 フーリエ変換 フーリエ逆変換 Ｇ（ｆ） から ｇ（ｔ） を合成 フーリエ逆変換 離散関数と数ベクトル 有限長の離散関数は有限数列 ｅ（ｋ） ＝ ｛ｅ（０），ｅ（１），ｅ（２），．．．ｅ（Ｎ－１）｝ ｇ（ｋ） ＝ ｛ｇ（０），ｇ（１），ｇ（２），．．．ｇ（Ｎ－１）｝ 数ベクトル ｅ = ｅ（０） ｅ（１） ･ ･ ｅ（Ｎ－１） ｇ= ｇ（０） ｇ（１） ･ ･ ｇ（Ｎ－１） 有限長の離散関数は、数ベクトルと等価 有限長の離散関数の内積 内積は数ベクトルと同様に （ｅ（ｋ），ｇ（ｋ）） = ｅ*（０）ｇ（０）＋ｅ*（１） ｇ（１） ＋ｅ*（2） ｇ（2） ＋… Ｎ-１ ｅ*（ｋ）・ｇ（ｋ） =Σ ｋ= ０ 各種ベクトルの内積 幾何ベクトル 数ベクトル 離散関数ベクトル 連続関数ベクトル 内積 は、２つの関数 の『積の総和』 連続関数では『積の積分』に対応 内 ベクトル 確率ベクトル 積 |ａ| |ｂ| ｃｏｓθ * ａ ｂ * ∫ａ （ｔ） ｂ（ｔ）ｄｔ Ｅ [ａ*･ｂ] Ｅ：期待値 ・それぞれは、幾何ベクトルと同じイメージで 考えることができる。 ・内積が０となる場合 「ａとｂは直交している」 と言う。 26 ディジタル領域の周波数変換 ２．時間領域と周波数領域 ・ ｚ変換 ① フーリエ変換 ・ ＤＦＴ（離散フーリエ変換） ② 関数の直交変換 （講習省略） ・ ＦＦＴ（高速フーリエ変換） ③ ディジタル領域の周波数変換 ④ 窓関数 (p.37) ｚ 変換の定義式 ディジタル信号 ＝ 数列 ∞ ｚ 変換 ｚ 変換の例 Ｘ(z) ＝ Σ ｘ(k) ｚ -k ｘ(k) = ｛ｘ(-2), ｘ(-1), ｘ(0), ｘ(1), ｘ(2), ｘ(3) ｝ k = -∞ 時間 → 周波数 ｘ(k)：離散信号 ｚ： 周波数を表す と呼ばれる。 ｚ 変換 × ｚ1 × ｚ0 複素変数 Ｘ(z)： 「 信号 ｘ（ｋ） の ｚ 変換 」 逆 Z × ｚ2 × × ｚ -k を乗じて Σ 総和をとる Ｘ(z) ＝ ｘ(-2) ｚ2＋ ｘ(-1) ｚ1 ＋ ｘ(k) ＝ １ ∮Ｘ(z) ｚ k-1ｄｚ ｊ 2π × ｚ-1 ｚ-2 ｚ-3 ｘ(0) ＋ｘ(1) ｚ-1＋ ｘ(2) ｚ-2＋ｘ(3) ｚ-3 ＝ Σｘ(k) ｚ -k ｚ 変換 ＝ ｚ の多項式 周波数 → 時間 逆 ｚ 変換の実際 フーリエ変換との関係 (1/2) ｘ●(2) ｘ(1) ◇ フーリエ変換 １ ｘ(k) ＝ ∮Ｘ(z) ｚ k-1ｄｚ ｊ 2π X（ｆ）＝∫ｘ(t) ｅ - j 2πｆ t dt ◇ 離散信号のフーリエ変換 実際には複素積分を解くことは少なく、 -k Ｘ(z) を ｚ の級数の形 Σｘ(k) ｚ に表して、 ｘ(t)＝ Ｘ(z) ∞ Σ ｘ(k) ｚ k = -∞ 各ｚ その他 ｘ(3)ｘ(4) ● ｘ(0) ｔ = k/ｆｓ k：整数 k ● ● なので -ｊ2πｆ /ｆｓ -kの係数 -k (k = -∞ ･･･ ∞) 0 ● k Ｘ (ｆ) =Σｘ(k) ｅ k -ｊ2πｆ/ｆｓ -ｊ2π2ｆ/ｆｓ ＝ ｘ(0) ＋ｘ(1)･ｅ ＋ｘ(2)･ｅ その係数を時間信号とする場合が多い。 級数展開 ｘ(k) ｘ(5) ● ｘ(k) ＋ｘ(3)･ｅ -ｊ2π3ｆ/ｆｓ ＋ｘ(4)･ｅ -ｊ2π4ｆ/ｆｓ ＋ｘ(5)･ｅ -ｊ2π5ｆ/ｆｓ ＋ｘ(6)･ｅ -ｊ2π6ｆ/ｆｓ ＋ ・・・ 27 フーリエ変換との関係 (2/2) フーリエ 変換 Ｘ (ｆ)=Σｘ(k) ｅ k 離散系の角周波数 Ω の導入 -ｊ2πｆ k/ｆｓ Ω＝ 2πｆ /ｆs ｅ ｊ2πｆ/ｆｓ ＝ ｚ ここで、 ｅ -ｊ2πk ｆ/ｆｓ と置くと -k ＝ｚ 信号の帯域は ｆ ＝ - ｆs/2 ～ 0 ～ ｆs/2 なので -k Σｘ(k) ｚ k ＝ ｘ(0) ＋ｘ(1) ｚ-1＋ ｘ(2) ｚ-2 それに対応して Ω ＝ -π ～ 0 ～ π ・ 標本化周波数 ｆs に依存しない規格化表現 ・ ｚ ＝ ｅ ｊ2πｆ/ｆｓ → ｅ ｊΩ の簡単表現 ＋ｘ(3) ｚ-3 ＋ｘ(4) ｚ-4 ＋ｘ(5) ｚ-5 ＋ｘ(6) ｚ-6 ＋ ・・・ = Ｘ(z) ｚ 変換 ｚ 変換 と フーリエ変換 （まとめ） Ｘ(z) =Σｘ(k) ｚ k ｚ 変換 ｚ ＝ｅ ｊΩ フーリエ 変換 ｚ 変換の意義 ◇ アナログの ラプラス変換に相当 （フーリエ変換の周波数（実数）を複素数へ拡張） -k ・ 離散系の伝達関数が得られ、 系の特性が把握できる （後述） ｅ ｊΩ ＝ ｚ Ｘ (Ω)＝Σｘ(k) ｅ k （ X（ｆ）＝∫ｘ(t) ｅ ・ 実際的な周波数分析より、フィルタ設計や システム解析などの理論面で有用 -ｊΩ （Ω＝ 2πｆ /ｆs） - j 2πｆ t dt ） ◇ ディジタル信号（＝数列）が多項式として 解析的に扱えるようになる (p.37) ｚ 変換の重要な定理 ｘ(k) ｘ(k-1) （証明） Z Z ＤＦＴ と FFT Ｘ(ｚ) ｚ-1Ｘ(ｚ) 信号の遅れ ＦＦＴ： Fast Fourier Transform （高速フーリエ変換） DFTを早く（少ない計算量で）計算する方法 Ｘ(ｚ)＝Σｘ(k)・ｚ-k k より、ｘ(k-1) の ｚ 変換は、 Σｘ(k-1)・ｚ-k ＝ Σｘ(m)・ｚ-m ｚ-1 k ＝ｚ-1・Ｘ(ｚ) m ＤＦＴ： Discrete Fourier Transform （離散フーリエ変換） 有限長の離散信号に対するフーリエ変換で、 有限個の離散スペクトルを得る k-1 = m k = m+1 計算の結果は、 DFT ＝ FFT 28 ＤＦＴのアナログ的解釈 (1) ＤＦＴの定義式 ◇ ディジタル的には、ＤＦＴは、 Ｎ個の数列（時間信号）と Ｎ個の数列（スペクトル）の関係 N-1 Ｘ (p)＝ Σｘ(k) ・ｅｘｐ（-ｊ2πk･p/N） k=0 k（時間） = ０，１，２，……Ｎ-１ （N個） p（周波数）= ０，１，２，……Ｎ-１ （N個） 逆変換 （I DFT， Ｉ ： Inverse) ｘ(k) ＝ １ N-1 Ｎ ΣＸ(p)・ｅｘｐ（ ◇ アナログ的理解を得るために、 ディジタル信号（＝時間信号数列）が アナログ的パルス波形の列（離散的信号）だと考える。 ｊ2πk･p/N） p=0 周波数スペクトルは周期的。 「アナログ的」 の意味 ＤＦＴのアナログ的解釈 (2) ◇ ＤＦＴ の スペクトルが離散的 ◇ ディジタル信号 ＝ 数列 ｛ｘ(0)，ｘ(1)，ｘ(2)，ｘ(3)，ｘ(4)，ｘ(5)，・・・ ｝ 時間信号 が 周期的 ◇ 時間軸上で表したパルス列 （アナログ的） ｘ(2) ｘ(1) ● ● 0 ◇ アナログ的にはＤＦＴは、 離散 かつ 周期的な信号と 離散 かつ 周期的なスペクトルとの 対応関係と考えられる。 ｘ(5) ● １ ２ ３ ４ ● ● ｘ(3) ● ｘ(0) ｘ(4) ５ k k：連続値 非整数時間はゼロ ◇ ディジタル特有の「回り込み」現象に対しては、 周期的だと考えたほうが、直観的理解が得やすい ＤＦＴのアナログ的解釈 (3) ｘ(ｋ) -N 0 時間信号が離散 → スペクトルが周期的 N ＤＦＴのアナログ的解釈 (4) 2N ｋ ｘ(ｋ) 周期的離散信号 N スペクトルが離散 → 時間信号もＮ点の信号が 周期的に繰り返されていると考える 0 -N N N 周期的離散スペクトル N N この対応関係 がDFT X(p) p (- ｆｓ/2) ｋ N点の時間信号に、 N点のスペクトルが対応 X(p) -N/2 2N 0 p N/2 -N/2 （ｆｓ/2) (- ｆｓ/2) 0 N/2 （ｆｓ/2) 29 スペクトルの表示 DFT と ｚ 変換の対比 X(p) p＝-N/2～N/2-1 時間信号 周波数 スペクトル ｚ 変換 任意個数の 離散信号 連続スペクトル ＤＦＴ Ｎ点の 離散信号 Ｎ点の 離散スペクトル 通常の表示 0 を中心に 正負の周波数 0 -N/2 N/2 (- ｆｓ/2) ＤＦＴの定義式 の計算結果 スペクトルの周期性 0 N/2 0 N (ｆｓ/2) p N/2 以降を 左側に移せば 0 を中心にした 表示になる -N/2 0 N/2 N p -N/2, -N/2+1,…,-2, -1, 0 0, 1, 2, …, N/2, N/2+1, …, N-2, N-1, (N) -N/2 0 N/2 p＝0～N-1 N p ＦＦＴ（DFT) 結果の周波数対応 -N/2, -N/2+1, …, -2, -1, (0) 周波数軸の表示の単位 周波数番号 ｐ＝ 複素共役 アナログ周波数 ＝ Ｘ(0) Ｘ(1) Ｘ(2) ・・・ Ｘ(N/2) 直流 p Nと0は 等しいので、 p ＝0 ～ N-1 の Ｎ個が 周期 X(p) N/2 同じ N スペクトルの巡回関係（数値対応） X(p) p＝-N/2～N/2-1 -N/2 p p＝0～N-1 -N/2 (- ｆｓ/2) N (ｆｓ/2) 2ｆｓ/N ｆｓ/2 ・・・ Ｘ(N-2) X(N-1) Ｘ (-2) Ｘ(-1) = Ｘ*(2) = Ｘ*(1) ｆｓ/N 0 2 ・・・ N/2 Ｘ(0) Ｘ(1) Ｘ(2) ・・・ Ｘ(N/2) 2ｆｓ ・・・ ｆｓ 直流 ｆｓ N N 2 標本間隔を１とした時 ＝ 0 の離散系周波数 (p/N) 1 N 2 N ・・・ 1 2 ◇ ＤＦＴの複素指数関数 ｅｘｐ｛ｊ 2π･(p/N)･k｝ において、 2π･(p/N) ＝ Ω とおくと、 ｅｘｐ｛ｊ Ω･k｝ となる。 離散系角周波数 Ω＝ 0 （N/2） を中心とした共役対称性 1 2π 4π N N ・・・ π ※ 周波数軸を Ω で表示することも多い。 30 離散系の周波数 ２．時間領域と周波数領域 ・・・・・ 周波数１/Ｎ （周期 Ｎ） ０ Ｎ 周波数２/Ｎ （周期 Ｎ/２） ・ ・ ・ ・ ・・・・・ ０ Ｎ ｋ ① フーリエ変換 ② 関数の直交変換 （講習省略） ｋ ③ ディジタル領域の周波数変換 周波数 １/２ （周期 ２） ④ 窓関数 ｋ 窓かけ （Windowing） 窓かけ の定式化 信号は長時間継続している しかし、DFT は、短時間のデータを対象とする 窓関数 ｗ(t) 原信号 ｘ(t) ｘw(t) ＝ 1 × 0 Ｔw Ｔw 長時間のデータから、 短時間のデータを 取り出すことを 「窓かけ」 という Ｔw ◇ 「窓かけ」は、原信号に窓関数 ｗ(t) を乗じたもの と定式化できる。 図では、（方形窓） ⎧1 w(t ) = ⎨ ⎩0 Ｔw：窓長 窓かけ によるスペクトルの変化 t0 ≤ t ≤ t0 + Tw (t0 = 0 の場合）の値。 ◇ 窓かけ（＝時間領域の乗算） は、 周波数領域のたたみ込みなので、 方形窓： w(t ) = 1 ハニング窓：w(t ) = 0.5 − 0.5 cos(2π t / Tw ) ハミング窓：w(t ) = 0.54 − 0.46 cos(2π t / Tw ) 時間 ブラックマン窓：w(t ) = 0.42 − 0.5 cos(2π t / Tw ) + 0.08 cos(4π t / Tw ) 1 方形 0 .8 0 .6 ハミング ハニング ブラックマン 0 .4 0 .2 0 50 その他 を示したが、いろいろな形の窓関数がある。 代表的な窓関数 0 t0 ≤ t ≤ t0 + Tw 100 Ｔw フーリエ変換 周波数 ｘw(t) ＝ ｘ(t) ｗ(t) 原信号 Ｆ Ｘw(f) ＝ × Ｘ(ｆ) ＊ 窓関数 Ｗ(ｆ) たたみ込み 窓かけ後の信号は、原信号のスペクトル Ｘ(ｆ) に、 窓関数のスペクトル Ｗ(ｆ) がたたみ込まれる 31 正弦波のスペクトルに及ぼす窓関数の影響 計算されるスペクトル 窓関数のスペクトル（１） 方形 正弦波 窓関数 F ［ｘw(t)］＝ F ［ｘ(t)］＊ F ［ｗ(t)］ たたみ込み フーリエ変換 0 -20 -20 -40 -40 -60 -80 窓関数のフーリエ変換（スペクトル）を Ｗ(f) とすると、 畳み込んだ結果は、 δ（ｆ－ｆ0）＊Ｗ(f) ＝ W（ｆ－ｆ0） となって、ｆ ＝ ｆ0 における パルス δ（ｆ－ｆ0）が、 窓関数のスペクトルの形 W（ｆ－ｆ0）に変化する。 -60 -50 ハニング -20 -20 -40 -40 ブラックマン 0 2 4 周波数 -60 -50 0 -80 50 周波数 -50 0 50 周波数 ハニング 方形 -70 -2 50 一般的信号に対し、２種の窓関数を使って、 下記スペクトルが得られた。 どちらが正しいスペクトルに近いのか？ -40 -4 0 ブラックマン 「正しい」スペクトル？ -20 -6 -50 方形窓は、正弦波信号の持つ周波数以外の周波数にも、影響を 及ぼす。 ハミング窓も-50ｄB程度であるが、影響を及ぼす。 -10 -80 -80 0 -80 ハミング -60 周波数 -60 方形 -50 50 0 窓関数のスペクトル（２） -30 0 ハニング 周波数 ｆ0 の正弦波のスペクトルは、 ｆ ＝ ｆ0 におけるパルス（δ関数） δ（ｆ－ｆ0） 0 ハミング 0 0 0 -20 -20 -40 -40 -60 -60 6 周波数 方形窓は、主成分が鋭い。ハミングは第１副成分が小さい ブラックマンは主成分が最も広い。主成分が広いと、周波数 分解能が低下する。 窓関数は、スペクトルを見る角度 -80 -80 -50 0 -50 50 0 50 窓の違いによる 周波数分析結果の違いの例（１） 信号は、1000Hz の正弦波と、-50dBの 1500Hz の正弦波の和 ２つの正弦波とは見えない ２つの正弦波が明確 ハニング 方形 -20 -20 [dB] 0 [dB] ◇ 同じ物でも、それを見る角度によって、 見える部分や見える形が違ってくる。 ◇ 角度を変えて見ているだけなので、 どちらも 「正しい」。 0 -40 -40 -60 -60 0 1000 2000 3000 周波数 [Hz] 4000 0 1000 2000 3000 周波数 [Hz] 4000 方形窓は、大きさの異なる周波数成分の分析には不向き 32 窓の違いによる 周波数分析結果の違いの例（２） 信号は、1000Hz と 1030Hz の 同じ大きさの２つの正弦波の和 ２つの正弦波が分離できている 分離できていない ハニング 0 -10 -10 -20 -20 -30 -40 特に、非定常信号の短時間分析・処理において重要 音声分析の例 窓長 短 (8ms) 窓長 長い (64ms) （高い時間分解能） （高い周波数分解能） ・ 大まかなスペクトル形状 ・ 基本周波数の変化 ・ ピッチ周期が見える セ カ イ ノ 1000 120 900 -30 周波数 [dB] [dB] 方形 0 スペクトル形状は窓の長さにも依存 -40 -50 -50 -60 800 -60 800 800 100 700 80 600 500 900 1000 1100 周波数 [Hz] 1200 60 900 1000 1100 周波数 [Hz] 1200 400 40 方形窓は、近接した周波数の分離測定に向いている 300 200 20 100 200 (注） ハ二ングで分離できる場合もある 800 1000 1200 20 40 60 80 100 120 140 時間 時間－周波数変換 （まとめ） ◇ どの窓関数が良いか、一般的な答は無い 用途に応じて、適した窓を選ぶべき アナログ フーリエ変換 ・ 帯域制限 ・ 時間離散化 ◇ 最も一般的なものは、 信号分析ではハニング窓 音声処理ではハミング窓 ◇ 窓関数を用いた分析・処理を行う場合、 何種類かの窓関数を用いて、結果を比較し、 適切なものを選ぶと良い ディジタル （理論向き） ディジタル （実用向き） ◇ あまり差が無い場合も多いが、 一度は比較してみるべき！ ｚ 変換 ［時間］ ： 離散，無限長 ［周波数］： 周期，連続 有限長時間信号 （Ｎ個のデータ） ◇ 窓長も変えたほうが良い （特に非定常信号） ＤＦＴ ［時間］ ： 離散，有限（周期） ［周波数］： 周期，離散 （参考） ＤＦＴの計算は内積 （付録） 各種変換のまとめ N-1 連続 離散 連続 フーリエ変換 フーリエ級数 離散 Ｚ変換 ＤＦＴ 時間 600 時間 窓関数 （まとめ） 周波数 400 フーリエ変換 フーリエ級数 Ｘ (p)＝ Σｘ(k) ・ｅｘｐ（-ｊ2πk･p/N） k=0 ｘ= Ｚ変換 ＤＦＴ 時間 周期 離散 離散・周期 周波数 離散 周期 離散・周期 ｘ（０） ｘ（１） ｘ（２） ･ ･ ･ ｘ（Ｎ-１） ｅp = ｅｘｐ（０） ｅｘｐ（ｊ2πp/N） ｅｘｐ（ｊ4πp/N） ･ ･ ･ ｅｘｐ（ｊ2(N-1)πp/N） * Ｘ (p)＝（ｅp，ｘ）＝ｅp ｘ 33 ＤＦＴの行列表現 （参考） （参考） ＤＦＴ ａ＝ｰ2π/N と置くと → ｅｘｐ（-ｊ2πk p/N）＝ ｅｊ kp･ａ １ ・・・ １ Ｘ（０） = １ １ １ ｅja ｅj2a ・ ・ ・ｅj (N-1)a Ｘ（１） Ｘ（２） １ ｅj2a ｅj4a ･ ・ ・ ･ ・ ・ ･ ・ ・ Ｘ(N-1) 1 ｅj(N-1)a ｅj2(N-1)a DFT 直交行列 （ユニタリ） ｘ（０） ｘ（１） ｘ（２） ･ ･ ･ ｘ(N-1) 周波数 スペクトル ＝ 0 500 1000 離散信号ベクトル 周波数成分ベクトル 1500 2000 周波数 -1 時間 離散信号 -1.5 0 0.01 0.02 0.03 0.04 0.05 ｅp(k)の波形（１） ＤＦＴ の 複素指数関数部 ｅp(k) ＝ ｅｘｐ｛ｊ2π･(p/N)･k｝ の実部波形 ｅｘｐ｛-ｊ2π･（p/N）･ k｝ ｐ=０ ｅｘｐ（ -ｊωｔ） → ｅｘｐ（ -ｊ2π ｆ ｔ） 周波数 同一ベクトルの 異基底軸表現 DFT 直交（ユニタリ）行列 p= ０，１，２，……Ｎ-１ は直交変換 ｅ0(k) ＝ １ ・・・・・ （直流） ０ ｐ=１ に対応して、 （周期 Ｎの複素指数関数） 周波数１/Ｎ ０ １ ２ …… Ｎ-１ Ｎ Ｎ Ｎ ｐ=２ 区間長 Ｎ で１周期 〃 ２ 〃 Ｎ k ） ｅｘｐ（ｊ２π N ｅｘｐ（ｊ２π ｋ ・・・・・ Ｎ ｋ ２k ） N （周期 Ｎ/２の複素指数関数） 周波数２/Ｎ ０ ・・・・・ Ｎ ｋ ｅp(k)の波形（２） ｅp(k) ＝ｅｘｐ｛ｊ2π･(p/N)･k｝ の実部波形 ｐ＝Ｎ/２ ｅｘｐ（ｊ2π k ） ２ （周期 ２の複素指数関数） 周波数 １/２ （上限周波数） ｋ ｐ＞Ｎ/２ Ｎ-ｐ の負の周波数に対応 → 折り返しひずみと同じ理屈 34 系（システム）とは？ ３．線形系（線形システム） 入力 ｘ(k) ① 線形系の概要 出力 ｙ(k) 系 ② インパルス応答 ③ たたみ込み 例） スピーカ 室内音響系、マイクロホン、通信系 自動車、人間、ｅｔｃ 線形系の定義 線形系 ｘ1(k) 系 線形系 非線形系 なぜ重要か？ ① 身近に多く存在 ② 扱いやすい ・数式表現（モデル化） が簡単 ・制御しやすい ｘ2(k) ｙ2(k) の時、 以下が成立 ｃ・ｘ1(k) 系 定数倍 ｘ1(k)＋ｘ2(k) ［代表例］ 受動電気回路 室内音響系 ｙ1(k) 系 和 ｃ・ｙ1(k) 定数倍 系 ｙ1(k)＋ｙ2(k) 和 (p.48) 線形系の重要な性質 正弦波入力に対する線形系の出力 一般の（広帯域）信号を入力した時には、 波形が大きく変化する系であっても、 ｘ (k) 線形系に正弦波信号を入力すると、 波形は変わらない。 同じ周波数で、 振幅と位相のみが変化した 正弦波を出力する。 ｙ (k) 時間 k 線形系 正弦波を入力した時には、 同じ周波数の正弦波を出力する。 ｘ (k) k ｙ (k) k k 注） ｃｏｓ波も正弦波に含める（以下同様） 線形系 k 35 クイズです 正弦波は線形系の固有関数 第７問 Ａｘ＝λｘ 行列 固有ベクトル （線形変換） 線形系 固有値 固有ベクトル 振幅と位相 の変化量 正弦波 代表的な線形操作は、次のうちどれ？ (a) 平方根 (b) 対数操作 (c) 微分･積分 正弦波 正弦波が固有関数（ベクトル） 正弦波の変化が重要な理由 (1) ◇ フーリエ変換 全ての信号 ｘ(k) は、正弦波 ｅ ｊΩk の荷重和として表される。 ｘ(k)＝ΣX(Ω)･ｅ ｊΩk 正弦波の変化が重要な理由 (2) ｅ ｊΩk ◇ 定数 （ X(Ω)) 倍の正弦波 X(Ω)･ｅ ｊΩk を入力 (1) Ω ◇ 一方、正弦波 ｅ ｊΩk を線形系に入力すると、 Ｇ(Ω)倍されて出力されるとする。 ｅ ｊΩk 線形系 Ｇ（Ω） Ｇ(Ω)・ｅ ｊΩk すなわち、 個々の正弦波 ｅ ｊΩk に対する 出力 Ｇ(Ω) ｅ ｊΩk がわかれば、 すべての入力信号 ｘ(k) に対する出力 ｙ(k) がわかる。 ｅ ｊΩk 線形系 Ｇ X(Ω)･ｅ ｊΩk X(Ω)・Ｇ(Ω)･ｅ ｊΩk Ｇ 定数倍 ◇ Ωで総和を取った信号 ｘ(k) を入力 ΣX(Ω)･ｅ ｊΩk ΣX(Ω)・Ｇ(Ω)･ｅ ｊΩk Ｇ Ω Ω ＝ ｘ(k) 正弦波の変化が重要な理由 (3) 正弦波 Ｇ(Ω)･ｅ ｊΩk Ｇ 線形系Ｇの固有値 ◇ Ｇ(Ω) は複素数であり、正弦波 ｅ ｊΩk の、 振幅・位相の変化を表す ◇ Ｇ(Ω) を系の 「周波数特性」 と呼ぶ 和 周波数特性の一般的定義 入力信号スペクトル Ｘ（Ω） 正弦波は、Ｇ(Ω)倍 Ｇ(Ω)・ｅ ｊΩk ＝ ｙ(k) 周波数特性 （定義） 線形系 Ｇ Ｙ(Ω) Ｇ(Ω) = ----Ｘ(Ω) 入力の周波数成分 Ｘ(Ω)が Ｇ(Ω) 倍されて出力される 出力信号スペクトル Ｙ(Ω) 出力スペクトルと 入力スペクトル との比 Ｙ(Ω) = Ｇ(Ω)･Ｘ(Ω) 36 伝達関数 線形系と正弦波（まとめ） 線形系 ｘ(k) ｙ(k) インパルス応答 ｇ(k) 伝達関数 Ｇ(z) Ｘ(z) Ｙ(z) 伝達関数 （定義） Ｇ(z) = ----Ｘ(z) 周波数領域の 入出力関係 Ｙ(z) 出力 ｚ 変換と 入力 ｚ 変換 との比 Ｙ(z) = Ｇ(z)･Ｘ(z) 入力にＧ(z) 倍したのが出力 非線形系 ◇ 正弦波を入力したら正弦波が出力されるのは、 正弦波が線形系（線形変換）の固有関数だから である。 （参考） 微分、積分も線形変換です ◇ その固有値は、周波数特性（伝達関数）である。 ◇ ある線形変形に対する固有値が得られれば、 全ての入力に対する変形が予想できる。 ◇ 正弦波や、 周波数特性の重要性は、以上の理由による 音で感じる非線形 ・ 線形系の性質が満足されない系。 ・ 厳密には、多くの系が非線形系。 非線形の程度が小さい系 （または、入力の範囲など）を 線形系とみなして信号処理を行う。 ・ 非線形系では高調波歪が発生することが多い。 （スピーカ系では、オーバーロードの時など 音が歪む、割れる） ・ ランダム（非周期）信号では、 非線形の発生は 「耳で」 気づきにくい。 ・ 非線形の程度が大きいと信号処理の 効率が低下するので、注意が必要である。 時不変系とは？ 時不変という性質の必要性 時間が経っても特性が変化しない系 ｘ (k) 系 ｙ (k) ｘ (k -τ) 系 ｙ (k -τ) 時変だと困ること ・ 計測：今測った特性が、もう変化している ・ 制御：推定した特性に基づいた制御ができない 時不変性 ＝ （系や信号の） 定常性 とも言う 時間τの後に、同じ入力を入れれば、 同じ出力が出てくる系 ただし、短時間の間に大きく変化しなければ、 ある程度の時変性や非定常性には 対応できる信号処理もある。 → 適応信号処理など 37 線形系の概要 （まとめ） ３．線形系（線形システム） とくに断らない限り、 時不変 な 線形系は、 ① 線形系の概要 信号処理の大前提 ② インパルス応答 ③ たたみ込み ◇ 正弦波は、線形系の固有関数 アナログ系のインパルス応答 （インパルス信号） δ(t) 幅が 0で高さが∞ 積分値が１のパルス インパルス応答の例 部屋 （音響系） δ(ｔ) ｔ 0 インパルス信号を入力した時の系の出力 δ(t) ｇ(ｔ) 0 ｔ ｇ(t) 系 インパルス応答 インパルス応答と周波数特性 250 ms インパルス応答測定と周波数特性の例 ◇ 系のインパルス応答 ｇ(t) のフーリエ変換は、 その系の周波数応答特性 dt 70 60 相対音響出力[dB] Ｇ（ｆ） ＝∫ｇ(t) ｅ - j 2πｆ t 80 測定風景 スピーカ マイク 50 40 30 Ｇ(f) 20 0° 30° 60° 90° 10 δ(t) 周波数応答特性 Ｇ（ｆ） 0 ｇ(t) ◇ 周波数応答特性Ｇ（ｆ）の逆フーリエ変換は、 インパルス応答 ｇ(t) → 証明は付録 スピーカの インパルス応答 100 2 .7 10k スピーカの 周波数特性 （指向性） ｇ(t) 2 .6 5 1k 周波数[Hz] 2 .7 4 x 10 38 ディジタル系のインパルス信号 （インパルス信号） δ(t) アナログ 幅が 0で高さが∞ 積分値が１のパルス ∞ δ(t) ｔ 0 時間 0 で値１、 その他の点では 値 0 の信号 １ δ(k) １ k 0 ◇δ(k) はδ(t) の物まね ではない。 ◇ ディジタルはアナログ ｔ を模倣したり、シミュレー トするもので、模倣に よる誤差（違い）が発生 する、という解釈は誤り。 ◇ δ(k) とδ(t) は等価で k あることを以下に説明 する。 0 （単位サンプル信号） δ(k) ディジタル δ(t) と δ(k) 0 δ(t) をディジタル世界から見ると 理想ローパスフィルタ の インパルス応答 １ 実世界 (Ｒｅａｌ Ｗｏｒｌｄ） アナログの世界 フーリエ 変換 １/ｆｓ δ(t) ‐fｓ/２ 理想 ＬＰＦ D/A 理想 ＬＰＦ A/D fｓ/２ ｆ 理想ＬＰＦの インパルス応答 （sinc 関数） 理想ローパスフィルタ を標本化したもの 標本化周期 Ts=1/fs の整数倍の時刻では 値が零になる 計算機の中の世界 ディジタルの世界 k 証明 は付録 -4Ts -2Ts 0 2Ts 4Ts -3Ts -Ts Ts 3Ts ( インパルス応答 ) （周波数特性） 理想ローパスフィルタ の インパルス応答 を標本化すると ｓin πｆｓ ｔ πｆｓ ｔ ｔ = 1/fs = Ts で、 πｆｓ ｔ =π （分子がゼロ） δ(t) と δ(k) の等価性 実世界 (Ｒｅａｌ Ｗｏｒｌｄ） アナログの世界 １ k -4Ts -2Ts 0 2Ts 4Ts -3Ts -Ts Ts 3Ts δ(k) はδ(t) を δ(t) 理想 ＬＰＦ ディジタル世界 から見たもの 理想 ＬＰＦ sinc 関数 D/A δ(k) 0 A/D 単位サン プル関数 1 0 k 計算機の中の世界 ディジタルの世界 δ(k) 39 ディジタル系におけるインパルス応答 （単位 サンプル応答） δ(k) ｇ(k) 系 k 0 k 単位サンプル信号δ(k) に対する応答 ｇ(k) を ディジタル系の インパルス応答と呼ぶことは妥当 インパルスは全ての周波数成分を含む ｆ0 2ｆ0 3ｆ0 4ｆ0 5ｆ0 6ｆ0 7ｆ0 8ｆ0 9ｆ0 10ｆ0 ∑ -0.025 -0.02 -0.015 -0.01 -0.005 0 0 0.005 0.01 0.015 0.02 インパルス応答の特質 インパルス応答は，線形系の 全ての情報を含んでいる 例えば、 ◇ インパルス応答がわかれば， あらゆる入力に対する出力がわかる ◇ インパルス応答の ｚ 変換は、 その線形系の伝達関数 ◇ インパルス応答の ＤＦＴ は、 その線形系の周波数応答特性 次節で説明 インパルス応答（まとめ） １） ディジタル系の単位サンプル信号δ(k)は アナログ系のインパルス信号δ(t)と等価 ２）ディジタル系の単位サンプル応答を インパルス応答と呼ぶ ３） インパルス応答は、線形系の基本量 0.025 全ての周波数成分を等振幅で 足すと、インパルス信号になる ３．線形系（線形システム） ① 線形系の概要 ② インパルス応答 ③ たたみ込み たたみ込み 線形系の出力 ｙ(k) は、 入力信号 ｘ(k) と、 系のインパルス応答 ｇ(k) との たたみ込み演算＊で表される。 y (k ) = x(k ) ∗ g (k ) = ∞ ∑ x ( m) ⋅ g ( k − m) m = −∞ ＜線形系特有の演算＞ 線形系の性質からたたみ込み演算の 必然性が導出できる 40 単位サンプル信号を用いた信号表現 ｘ(k) ｘ(0) ｘ(0) ｘ(1) 0 ｘ(2) ＝ 2 1 k k 0 ＋ 2 ＋ ｘ(2)δ(k-2)＋・・・ k m ｘ(1)δ(k-1) 入力信号 ｘ(k) は、 単位サンプル信号δ(k) の時間シフトδ(k-m)、 ｘ(2) k ｘ(k) ＝ｘ(0) δ(k) ＋ ｘ(1)δ(k-1) ＝Σｘ(m)δ(k-m) k ＋ m ｘ(0)δ(k) ｘ(1) 1 δ(k-m) 1 0 単位サンプル信号を用いた信号表現（２） ｘ(2)δ(k-2) 定数（ ｘ(m) ）倍、その総和 として表される。 (p.48) 線形系の入出力関係 δ(k) Ｇ ｇ(k) ｘ(k) Ｇ ？ たたみ込み x(k ) : 系に対する入力信号 g (k ) : 系のインパルス応答 δ(k-m) Ｇ ｇ(k-m) ｘ(m)δ(k-m) Ｇ ｘ(m) ｇ(k-m) 定数倍 Σｘ(m)δ(k-m) Ｇ Σｘ(m) ｇ(k-m) 総和 m ＝ ｘ(k) 時不変性 m 時不変線形系の３つの性質を利用 （定理） たたみ込み の 可換性 y (k ) = ∞ y ( k ) = m = −∞ ｙ(k) ＝ ｘ(k)＊ｇ (k) ＝ ｇ(k)＊ｘ(k) インパルス応答が ｇ(k) である系に、 信号 ｘ(k) を入力したときの出力 ＝ インパルス応答が ｘ(k) である系に、 信号 ｇ(k) を入力したときの出力 ∑ x ( m) g ( k − m) m = −∞ インパルス応答 ｇ を知れば、任意の入力 に対する 出力 ｙ が計算で求められる ｘ インパルス応答の因果性を 考慮した表示 ∞ ∑ x ( m ) g ( k − m) = ∑ g ( m ) x ( k − m ) m = −∞ y (k ) : 系の出力信号 ∞ y (k ) = ∞ ∑ x ( m) g ( k − m ) = m = −∞ ∞ ∑ g ( m ) x ( k − m) m = −∞ インパルス応答の因果性 → g (m) = 0 for m < 0 y (k ) = k ∞ m = −∞ m=0 ∑ x ( m) g ( k − m) = ∑ g ( m) x ( k − m ) 41 たたみ込みの解釈 たたみ込み の例 さらに、 x(m) = 0 for m < 0 を仮定すれば、 k m =0 m =0 ∑ x ( m ) g ( k − m) = ∑ g ( m ) x ( k − m) 1 ｇ(2) 大 中 小 0 1 2 は は は ね ね ね だ だ x(2) 0 1 2 ｇ(1) x(1) 2 y (2) = x(0) g (2) + x(1) g (1) + x(2) g (0) ｇ(0) だ 0 x(0) ね （ 入力） ｘ k =2 は y ( k ) = k インパルス 応答 ｇ y(0) インパルス応答の ｚ 変換が 伝達関数となることの証明 y(1) 時間 k だ y(2) ｚ 変換（多項式）の積の例 例） Ｘ(ｚ)＝ ｘ(0)＋ｘ(1)ｚ-1＋ｘ(2)ｚ-2 線形系 ｘ(k) インパルス応答 ｇ(k) 伝達関数 Ｇ(z) Ｘ(z) 時間領域の 入出力関係 ｙ(k) Ｙ(z) ｇ(k)＊ｘ(k) = ｙ(k) ｚ 変換 周波数領域の 入出力関係 フーリエ変換と 同様の性質 証明 → Ｇ(z)･Ｘ(z) = Ｙ(z) ｚ 変換の積とたたみ込み ｚ 変換は多項式 多項式の積の係数は、たたみ込み Ｇ(ｚ)＝ ｇ(0)＋ｇ(1)ｚ-1＋ｇ(2)ｚ-2 Ｇ(ｚ) Ｘ(ｚ) ｙ(1) ｙ(0) ＝ ｇ(0)ｘ(0) ｚ0 ＋（ｇ(0)ｘ(1)＋ｇ(1)ｘ(0)）ｚ-1 ＋（ｇ(0)ｘ(2)＋ｇ(1)ｘ(1)＋ｇ(2)ｘ(0)）ｚ-2＋・・・ ｙ(2) 多項式の積の結果、各項は たたみこみの形になっている ｚ 変換の積は、たたみ込みのｚ 変換 ｇ(k) ＊ ｘ(k) ｙ(k) ｚ 変換 ｚ 変換の積の係数は、たたみ込み ｚ 変換の積は、たたみ込みの ｚ 変換 Ｇ(ｚ) × Ｘ(ｚ) Ｇ(ｚ)Ｘ(ｚ) Ｙ(ｚ) 一致 ｘ(k)、ｇ(k) の ｚ 変換の積 Ｇ(z)Ｘ(z) と、 たたみ込みの結果 ｙ(k) の ｚ 変換 Ｙ(z) とが 一致する 42 ＤＦＴの積 （重要 ！） ｚ 変換の重要な性質、について １） ｙ(k) ＝ｇ(k)＊ｘ(k) Ｙ(ｚ) ＝Ｇ(ｚ) Ｘ(ｚ) DFT の積は、時間領域では円状（巡回）たたみ込み （ｚ 変換の積） （時間領域のたたみ込み） ２）系のインパルス応答 ｇ(k) の ｚ 変換が伝達関数 Ｇ(ｚ) となる。 ∵ 1) より、 Ｙ(z)＝Ｇ(z)･Ｘ(z) であり、 フーリエ変換（アナログ）や ｚ 変換の積は、 直線たたみ込み（通常のたたみ込み）に対応 直線たた み込み 入力がインパルス信号のとき、Ｘ(z)＝１ であるので、 インパルス応答（系の出力）の ｚ変換 Ｙ(z)＝伝達関数 Ｇ(z) ３） 伝達関数Ｇ(ｚ)に、ｚ＝ｅ j Ω を代入すると、 周波数特性が得られる。 ＊ k Ｎ ＝ Ｎ Ｎ 円状たた み込み k k ＊ k ＝ k k Ｎの長さからはみ出た分は回り込む Ω：離散角周波数 ※ 円状たたみ込みに関しては付録を参照 周波数変換とたたみ込みの関係 Ｚ変換の積 直線たたみ込み （実存する物理現象） ＤＦＴの積 円状たたみ込み 系のインパルス応答と周波数特性 畳込み ｙ(k)＝ｘ(k)＊ｇ(k) 直線たたみ込み フーリエ変換 （ＤＦＴ） ｚ変換 （やや人工的な現象） （参考） フーリエ変換 の積 インパルス応答 ｇ(k) 逆ｚ変換 伝達関数 Ｇ(z) 逆ＤＦＴ ｚ＝ｅ ｊΩ 周波数特性 G(ｅ ｊΩ) 積 Ｙ(z)＝Ｘ(z)Ｇ(z) インパルス応答とたたみ込み（まとめ） ◇ 線形系のインパルス応答 （単位サンプル応答） ｇ(k) を知れば 任意の入力信号 ｘ(k) に対する系の出力ｙ(k) が、たたみ込み演算によって計算できる。 ｙ(k) ＝Σｇ(m) ｘ(k-m) m ◇ インパルス応答 ｇ(k) の ｚ 変換 Ｇ(ｚ) が伝達関数、 ＤＦＴが周波数特性 Ｇ(ｅ ｊΩ )となる。 ◇インパルス応答からは、その他にも さまざまな特徴量（残響時間、ＭＴＦなど）が得られる。 43 ① ディジタルフィルタとは ４．ディジタルフィルタ ディジタルフィルタ （広義） ① ディジタルフィルタとは アナログ 信号 ② フィルタの種類 と 伝達関数 ③ 極と零によるフィルタの 特性解析 Ａ/Ｄ 数列 ディジタル フィルタ （狭義） アナログ 信号 Ｄ/Ａ 数列 ④ フィルタの実行 ⑤ フィルタの設計 （Ａ/Ｄ、Ｄ/Ａは、低域フィルタを含んでいる） ディジタルフィルタの構成要素 遅延器 最も簡単なディジタルフィルタの例 入力 出力 ｘ(k) ｙ(k) ＋ ｘ(k) ｚ-1 ｘ(k-1) ｃ･ｘ(k-1) × ｃ (フィルタ係数) 構成要素： × 乗算器 ＋ 加算器 ｘ(k) ｚ-1 ｘ(k) ｘ(k-1) ｘ(k-1) 遅延器とは， 現在の信号ｘ(k) を 入力し， １時刻前の信号ｘ(k-1) を出力する。 (単位)遅延器 遅延器 の伝達関数 アナログフィルタ の特性表現 L ｘ(k) 遅延器 ｚ変換 Ｘ(ｚ) ｘ(k-1) 構成要素は、 ＬＣＲ ｅi R ｚ変換 遅延器の 伝達関数 ｚ-1 ｚ-1Ｘ(ｚ) C 入出力関係の 表現は 微分方程式 解析は ラプラス変換 ｖ＝L ｅo ｄ2 ｑ ｄｑ ＋ R ｄｔ ＋ 1 ｑ C ｄｔ2 Ｖ ＝ ｓ2 L Ｑ ＋ｓ R Ｑ ＋ 1 Ｑ C 44 フィルタの伝達特性 ディジタルフィルタの特性を計算 入力 出力 ｘ(k) ｙ(k) ＋ ｚ 変換 Ｙ（ｚ）＝ 伝達関数 ｚ-1 ｘ(k-1) ｘ(k) 入出力関係は 差分方程式 ｃ･ｘ(k-1) × ｃ Ｇ（ｚ）＝ (フィルタ係数) 周波数特性 y(k) = Ｇ(Ω)＝ 解析は Z変換 ｃ (p.42) フィルタの周波数特性 フィルタの周波数特性図 ＝１ ｃ Ｇ(Ω)＝１＋ｅ -ｊΩ ＝ ｅ -jΩ/２（ ＝ ｅ -jΩ/２ ・２・（ ） ） ＝１ | Ｇ(Ω) | ＝２・|cos(Ω/２)| Ｇ(Ω) ｙ(k) ＝ｘ(k) ＋ｘ(k-1) 6 dB 直流 振幅特性 |Ｇ(Ω)| ＝ ２・| | Ω -π ～ ０ ～ π ｆ -ｆｓ/2 ～ ０ ～ ｆｓ/2 フィルタの周波数特性図 (2) ｃ ＝ -１ | Ｇ(Ω) | ＝２・|ｓｉｎ(Ω/２)| Ｇ(Ω) 直流 ｆｓ/2 0 π （高域通過フィルタ） （ｆｓ/2） （低域除去フィルタ） 0 π （低域通過フィルタ） （ｆｓ/2） Ω ｘ(k) ｙ(k) ＋ ｘ(k) ｚ-1 ｘ(k-1) ｙ(k) ＝ｘ(k) ＋ｃ･ｘ(k-1) 差分方程式を計算する プログラム例 特性はフィルタ係数に依存 ｆｓ/2 Ω ディジタルフィルタの実体はプログラム ｙ(k) ＝ｘ(k) －ｘ(k-1) 6 dB -40 -40 × ｃ ｃ･ｘ(k-1) 入力 出力 Function Filt(xk, yk) double xk, xk1, c, yk yk=xk+c*xk1 xk1=xk return 45 ディジタルフィルタとは、（まとめ） ４．ディジタルフィルタ ◇ 乗算・加算・遅延・フィルタ係数 ◇ フィルタ特性の求め方 ・ 差分方程式 → ｚ 変換 ・ 伝達関数 Ｇ（ｚ）＝Ｙ（ｚ）/Ｘ（ｚ） 周波数特性 Ｇ(Ω) ・ ｚ＝ｅ ｊΩ ① ディジタルフィルタとは ② フィルタの種類 と 伝達関数 ③ 極と零によるフィルタの 特性解析 ④ フィルタの実行 ◇ 周波数特性は フィルタ係数に依存 ◇ ディジタルフィルタはプログラム ⑤ フィルタの設計 Ｆ Ｉ Ｒ フィルタ ２種類のフィルタ δ(k) Ｆ Ｉ Ｒ（Finite Impulse Response ） フィルタ 有限 ｚ-1 インパルス応答 Ｉ Ｉ Ｒ（Infinite Impulse Response ） フィルタ 無限 k 0 インパルス応答 インパ ルス 応答 ｙ(k) ＋ 入力を遅ら せて荷重和 × ｃ k ・・・ -2 -1 0 1 2 3 ・・・ 入力 δ(k) ・・・ 0 0 1 0 0 0 ・・・ 出力 y(k) ・・・ 0 0 1 c 0 0 ・・・ 時刻 有限長 (p.44) Ｆ Ｉ Ｒ（Finite Impulse Response ） フィルタ （または、非巡回型フィルタ） Ｉ Ｉ Ｒフィルタ δ(k) ｙ(k) ＋ 出力を遅ら せて荷重和 インパ ルス 応答 不安定なフィルタ × ｃ ｕ(k) ｙ(k) I I R フィルタ ｚ-1 ｙ(k-1) ｃｙ(k-1) 時刻 k ・・・ -1 0 1 2 3 ・・・ 入力 δ(k) ・・・ 0 1 0 0 0 ・・・ 出力 y(k) ・・・ 0 1 c c２ c３ ・・・ ｚ-1 ｃ インパルス応答は ｛1, c, c2, c3, …｝ 例 |c | ≧１ の時、不安定（発散） 無限長 Ｉ Ｉ Ｒ（Infinite Impulse Response ） フィルタ （または、巡回型フィルタ） |c | ≧１ の時、不安定（発散） 46 一般形 （参考） 不安定な 「系」 の例 ハウリング 不安定なアナログ系 b0 x(k) ｘ(k) ｙ(k) ｚ-1 b1 a1 ｚ-1 x(k-1) マイクロホン アンプ ピ～～ b2 a2 ｚ-1 x(k-2) y(k-2) ｚ-1 bQ x(k-Q) 不安定系は使い物にならないことがわかる aP ｚ-1 y(k-P) P Q ｉ =1 Q ｉ=0 ｙ(k) ＝Σａi ｙ(k-i) ＋Σｂiｘ(k-ｉ) IIR ｙ(k) ＝Σｂi ｘ(k-i) ＦIR （ａ1,･･･ ａP＝0） ｉ=0 （用語の補足）Ｆ Ｉ Ｒフィルタの次数とタップ数 Q ｘ(k) b0 x(k) ｙ(k) ｉ =0 ｚ-1 b1 「次数」 伝達関数の 多項式の次数 遅延器の数 ｚ-1 b2 x(k-2) 「タップ数」 係数 ｂi の数 ｚ-1 bQ x(k-Q) P Q ｉ =1 ｉ =0 Q Ｇ(z) ＝ ｂ0Π(1- ｑi ｚ-1) ｉ =1 x(k-1) フィルタ（一般形）の伝達関数 ｙ(k) ＝Σａi ｙ(k-i) ＋Σｂiｘ(k-ｉ) ｙ(k) ＝Σｂi ｘ(k-ｉ) P Q ｉ =1 ｉ =0 P (1- Q B( z ) = ∑ bi z = b0 + b1 z + b2 z + L + bQ z −2 −Q i =0 に対して、 B( z ) = 0 を満たす ｚ の値 q を 「根」 と呼ぶ。 ② 多項式の根の数は、複素根、重根も含めると、 多項式の次数 （Ｑ次多項式はＱ個） の根を持つ 移項 Q Ｇ(z) ＝ Y(ｚ) ＝ Ｘ(ｚ) Σ ｂ ｚ-ｉ ｉ =0 i P 1- Σ ａ ｚ-ｉ ｉ =1 i 伝達関数 多項式に関する性質 (2) ① 多項式 −1 ｚ 変換 Σ ａ ｚ-ｉ) Y(z) ＝ （Σｂi ｚ-ｉ )Ｘ(ｚ) ｉ =1 i ｉ =0 多項式に関する性質 (1) −i 差分方程式 Y(z) ＝Σａiｚ-ｉ Y(ｚ) ＋Σｂiｚ-ｉ Ｘ(ｚ) この場合は、 Ｑ次、 Ｑ＋１タップ Q すべての ディジタル フィルタが この差分 方程式で 書ける y(k-1) ｚ-1 ③ 多項式は、Ｑ個の根 ｑi （ ｉ=1,2,･･･,Q ）を用いて 因数分解した形で表すことができる。 B( z ) = b0 (1 − q1 z −1 )(1 − q2 z −1 )(1 − q3 z −1 ) L (1 − qQ z −1 ) Q = ∏ (1 − qi z −1 ) i =0 （ ｑi が根であることの確認 ） q 1 − qi z −1 = 1 − i であるので、 z = qi で、 z B（ｚ） はゼロになる。 47 伝達関数の因数分解 極 と 零 Q Ｇ(z) ＝ Y(ｚ) ＝ Ｘ(ｚ) Σ ｂ ｉ =0 i ｚ-ｉ Q Σ ａ ｚ-ｉ ｉ =1 i Ｇ(z) ＝ ＝ ｂ0 ＋ ｂ1ｚ-1 ＋ ｂ2ｚ-2 ＋ｂ3ｚ-3 ＋ ・・・ ＋ｂQｚ-Q 1- ａ1ｚ-1 - ａ2ｚ-2 -ａ3ｚ-3 - ・・・ -ａPｚ-P ＝ ｂ0 ･(1- ｑ1 ｚ-1)･(1- ｑ2 ｚ-1)・・・(1- ｑQ ｚ-1) (1- ｐ1 ｚ-1)･(1- ｐ2 ｚ-1)・・・(1- ｐP ｚ-1) ｚ = ｑi で 0 P Π (1- ｐi ｚ-1) ｉ =1 ｐi (1- － ｚ ) Ｇ(z)が 0 ｚ = ｐi で 0 Ｇ(z)が∞ 極： Ｇ(z)の値を ∞とする ｚ の値 Q ＝ ｂ0Π(1- ｑi ｚ-1) ｉ =1 P 1- ｂ0Π(1- ｑi ｚ-1) ｉ =1 ｑi ： 分子多項式の根 ｐi ： 分母多項式の根 P Π (1- ｐi ｚ-1) ｉ =1 pi 零： Ｇ(z)の値を ０ とする ｚ の値 ｑi 周波数の p と混乱しない ４．ディジタルフィルタ ① ディジタルフィルタとは ③ 極と零に基づくフィルタの特性解析 (1) フィルタの安定条件 伝達関数の部分分数展開 ② フィルタの種類 と 伝達関数 (2) 周波数特性の定性的理解 ③ 極と零によるフィルタの 特性解析 伝達関数の ｚ 平面上表現 (3) フィルタの分類 ④ フィルタの実行 ⑤ フィルタの設計 フィルタの安定条件 （ I I R フィルタ） Q Ｇ(z) ＝ ｂ0Π(1- ｑi ｉ =1 ｚ-1) 伝達関数の部分分数展開 ｂ0 ･(1- ｑ1 ｚ-1)･(1- ｑ2 ｚ-1)・・・(1- ｑQ ｚ-1) (1- ｐ1 ｚ-1)･(1- ｐ2 ｚ-1)・・・(1- ｐP ｚ-1) Ｇ(z) ＝ P Π (1- ｐi ｚ-1) ｉ =1 N ＝ Σ C ｚ-ｉ ｉ =0 i ＝ Σ C ｚ-ｉ ｉ =0 i Ｐ ＋ N フィルタが安定である（発散しない）条件は、 その伝達関数 Ｇ(z) の すべての極 pi （i = 1,2,･･･,P）の絶対値が ｜pi｜＜１ FIR 部分 ＋ Σ ｉ =1 Di 1- ｐi ｚ-1 Ci，Di ： 定数 D1 D2 ＋ ＋ 1- ｐ1 ｚ-1 1- ｐ2 ｚ-1 DP ・・・ ＋ 1- ｐP ｚ-1 一次極特性 （ IIR 部分 ） 48 伝達関数の和は並列接続 部分分数展開のブロック図 Ｇ（ｚ） Ｇ（ｚ） N Ｇ2（ｚ） Ｙ（ｚ） G2X Ｘ（ｚ） Y=（G1X＋G2X） G= Y/X Ｇ（ｚ）＝Ｇ1（ｚ）＋Ｇ2（ｚ） (1- ｐi ｚ-1) −1 (1 − pi z ) ⋅ V = U ( z ) V ( z ) = U ( z ) + pi z −1V ( z ) v ( k ) = u ( k ) + pi v ( k − 1) 逆 ｚ 変換 １個の FIR部と P 個の 一次極フィルタ の並列接続 v ( k ) = u ( k ) + pi v ( k − 1) ｕ(k) 1 V ( z) = U ( z) 1 − pi z −1 1 1- ｐ2 ｚ-1 極はフィードバックゲインを表す 出力 V(z) 1 Y（ｚ） 1 1- ｐP ｚ-1 一次極フィルタ 入力 Ｕ(z) 1 1- ｐ1 ｚ-1 ・・・ Ｇ1（ｚ） ・・・ Ｘ（ｚ） Σ C ｚ-ｉ ｉ =0 i G1X 一次の I I R フィルタ フィルタの安定条件 インパルス応答は ｛1, p, p2, p3, …｝ ｖ(k) ＋ × ｚ-1 極は等比 数列の公比 ｐi フィードバック信号（巡回信号）は ｐi 倍 （極倍）される |ｐi | ≧ 1 だと信号が発散（不安定） ③ 極と零に基づくフィルタの特性解析 Ｇ（ｚ） (1) フィルタの安定条件 N Σ C ｚ-ｉ ｉ =0 i Ｘ（ｚ） 1 1- ｐ1 ｚ-1 ・・・ ・・・ 1 1- ｐ2 ｚ-1 1 1- ｐP ｚ-1 伝達関数の部分分数展開 Y（ｚ） どれか一つでも 不安定だと フィルタは 不安定 (2) 周波数特性の定性的理解 伝達関数の ｚ 平面上表現 (3) フィルタの分類 安定条件は、 すべての極 pi が、 ｜pi｜＜１ 49 極と零に関する一つの性質 Ｇ(z) ： 複素数 ｚ に対して、値を与える複素関数 ｚ-平面： 複素数 ｚ に対応した平面 |Ｇ(z)| は、ｚ-平面上に等高線表示できる。 Q ｂ0Π(1- ｑi ｚ-1) ｉ =1 Ｇ(z) ＝ 複素平面（ｚ-平面）で見る伝達関数 P Π (1- ｐi ｚ-1) ｉ =1 ◇ Ｇ(z) の分子分母は実係数多項式なので、 共役数はともに根。 Imag(z) ｚ-平面 ｘ(t) |Ｇ(z)| 共役 t 例） ａ＋ｊｂ が極（零）なら、ａ-ｊｂ も 極（零） Real(z) ：極 （山） (参考) 実関数 ：零 （谷） 複素指数関数 ｅ ｊΩ ｅ ｊΩ＝cosΩ＋ｊ ｓｉｎΩ ｅ ｊΩ は単位円の上を動く Ω：実数 単位円： 半径が１の円 実数成分 cosΩと、虚数成分 ｓｉｎΩを持つ複素数 ｅ ｊΩ Imag 角度がΩで 原点からの 距離が１ ｊ ｓｉｎΩ Ω 0 フィルタの周波数特性 Ｇ(z)に ｚ＝ｅ ｊΩを 代入し、Ωの値を -π～ 0 ～πと 変化させたものが 周波数特性 Imag(z) ｊ Ｇ(ｅ ｊπ) z=ｅ ｊΩ Ω＝0 Ω -1 1 Real(z) Ｇ(ｅ-ｊπ) Ω＝-π/2 伝達関数と周波数特性の例 例) ２つの極と２つの零とを持った伝達関数 Ｇ(z) ＝ (1- ｑ1 ｚ-1) (1- ｑ2 ｚ-1) (1- ｐ1 ｚ-1) (1- ｐ2 ｚ-1) ｚ＝ｑ1, ｑ2 で Ｇ(z) → 零 単位円上の 角度が角周波数 -ｊ ｚ＝ｐ1, ｐ2 で Ｇ(z) → 無限大 1 Real(z) -ｊ ｅ ｊΩ Ω＝π を考える。 Ｇ(ｅ ｊ0 ) -1 Ω＝π/2 Ω＝-π ◇ Ωの値によらず、|ｅ ｊΩ| は常に 1 Ｇ(ｅ ｊπ/2) ｊ Real cosΩ -π～ 0 ～π と変化させると 単位円上を一周 Imag(z) ｅ ｊΩは、 1 ｅ ｊΩ のΩの値を Ｇ(z)の値が 極・零 を 図上に表すと、 周波数特性 50 極・零と周波数特性 単位円 との関連 ｚ-平面 |Ｇ(z)| の等高線 ｚ-平面 Ω＝π/2 (1) 安定条件 極 ｜pi｜＜１ ：極 ：零 Imag Ω2 pi が単位円の 内部に含まれる。 Ω＝0 ｐ1 ｑ1 Real ｑ2 |Ｇ(Ω)| 周波数特性 Ω1 ｐ1 Ω＝π Ω＝ 0 ｑ1 Ω Real ｑ2 (2) Ωは角周波数 アナログ周波数と は、π → ｆs/2 ｐ2 Ω＝π ：極 ：零 Imag ｐ2 0 Ω1 Ω2 π 極の近くで山、零の近くで谷 周波数特性は極と零の分布で決定 |Ｇ(z)| の等高線を 図上に重ねると、 伝達関数の３次元表示 ③ 極と零に基づくフィルタの特性解析 |Ｇ(Ω)| (1) フィルタの安定条件 2 1.5 伝達関数の部分分数展開 1 0.5 (2) 周波数特性の定性的理解 0 -0.5 伝達関数の ｚ 平面上表現 -1 -1.5 (3) フィルタの分類 -2 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2 -1.5 -0.5 -1 0 0.5 1 1.5 2 極と零 によるフィルタの分類 極・零フィルタ Q １） Ｇ(z) ＝ ｂ0Π(1- ｑi ｚ-1) ｉ =1 P 1- ２） Π (1- ｐi ｉ =1 極・零フィルタ ｚ-1) （Ｉ Ｉ Ｒ） Π (1- ｐi ｚ-1) ｉ =1 全極フィルタ （Ｉ Ｉ Ｒ） ｂ0 Ｇ(z) ＝ P 1- Q Ｇ(z) ＝ ｂ0Π(1- ｑi ｚ-1) ｉ =1 P Π (1- ｐi ｚ-1) ｉ =1 |Ｇ(Ω)| Q ３） Ｇ(z) ＝ ｂ0Π(1- ｑi ｚ-1) ｉ =1 全零フィルタ （Ｆ Ｉ Ｒ） Ω 0 Ω1 Ω2 π 例） 極の周波数 Ω1と零の周波数 Ω2とを近づけるこ とで急峻な（遮断 特性の良い）フィ ルタを作ることが できる 51 全極フィルタ 全零フィルタ （＝ FIR フィルタ） |Ｇ(Ω)| |Ｇ(Ω)| ｐ1 Ｇ(z) ＝ ｂ0 Q Ｇ(z) ＝ ｂ0Π(1- ｑi ｚ-1) ｉ =1 ｐ3 P Π (1- ｐi ｚ-1) ｉ =1 Ω 0 π ・ 極によるスペクトルの表現 ・ ピークを表現 → 共振系の表現に適合、知覚に適合 ・ 音声生成のモデル：ＡＲ（Auto Regressive)モデル ③ 極と零に基づくフィルタの特性解析 （まとめ） (1) フィルタの安定条件 すべての極が ｚ 平面の単位円内 (2) 周波数特性 フィルタの周波数特性は、ｚ 平面 においてＧ(z) の 単位円上の値 (3) フィルタの分類 極・零フィルタ、全極フィルタ、全零フィルタ ４．ディジタルフィルタ ① ディジタルフィルタとは ② フィルタの種類 と 伝達関数 ③ 極と零によるフィルタの 特性解析 ④ フィルタの実行 ⑤ フィルタの設計 q1 q3 π 0 ・ 零によるスペクトルの表現 ・ F I Rフィルタであるので、安定性が保証！ （係数を、どのように定めても発振の危険性が無い) ・ 次数を高くすれば急峻なフィルタも作れる ・ ＭＡ（Moving Average)モデル Ω （補足） ＩＩＲ と ＦＩＲ の比較 ◇ IIR （長所） 少ない次数で急峻なフィルタ （短所） 安定性に注意を払う必要 プログラムが複雑 （用途） 固定的・専用 ◇ ＦＩＲ ＩＩＲ の逆 （長所） 安定性が保障 （短所） 急峻なフィルタは大きな次数 （用途） 実験室、可変フィルタ、適応フィルタ フィルタの実行 基本的には、差分方程式 P Q ｉ =1 ｉ =0 ｙ(k) ＝Σａi ｙ(k-i) ＋Σｂiｘ(k-ｉ) を、ＤＳＰやＣＰＵのプログラミングで実行 （ただし、ＦＩＲ フィルタには いくつかの 便利な性質がある） → 関連事項も交えて説明する 52 Ｆ Ｉ Ｒフィルタとインパルス応答 ｘ(k) b0 x(k) ｚ-1 ｙ(k) b1 ◇ フィルタ係数 b0， b1，b2，… がインパルス応答 ◇ 出力ｙ(k) は ｘ(k) と ｂi とのたたみ込み x(k-1) ｚ-1 b2 x(k-2) Q ｙ(k) ＝Σｂi ｘ(k-ｉ) ｉ =0 ｚ-1 bQ x(k-Q) Ｆ Ｉ Ｒフィルタのベクトル表現 Q ｙ(k) = ｂTｘ(k) ｙ(k) ＝Σｂi ｘ(k-ｉ) ｉ =0 ｂ= ｂ０ ｂ１ ｂ２ ･ ･ ･ ｂＱ ｘ(k) = ｘ(k) ｘ(k-1) ｘ(k-2) ･ ･ ･ ｘ(k-Q) T：転置 たたみ込みは、ｂ と ｘ(k) の 内積 演算 Ｆ Ｉ Ｒフィルタのたたみこみ計算 （フィルタの実行の標準形) たたみ込みは計算量が多い ｂ ｘ(k) ｘ(k) ｘ(k) ｘ(k-1) ｘ(k-2) ｘ(k) ｘ(k-1) ｘ(k) ｘ(k-1) ｘ(k-2) ･ ｙ(k) = ｂTｘ(k) ･ ｘ(k-Q) k=ｋ＋１ ｙ(k) たたみ込み と Ｆ Ｆ Ｔ （＝ＤＦＴ） ･ ･ ｘ(k-Q) ｘ(k-Q+1) Real Time 処理 関連事項： 信号と信号の乗算 信号どうしの乗算は、 信号の周波数範囲を拡大する！ ◇ 正弦波の積 ｓｉｎ 2πｆ1ｔ ×ｓｉｎ 2πｆ2ｔ ＝ ｛cos 2π(ｆ1-ｆ2)ｔ－cos 2π(ｆ1+ｆ2)ｔ｝/２ ◇ 一般に、 ａ(t) に含まれる 上限周波数 ｆma ｂ(t) に含まれる 上限周波数 ｆmb ａ(t)×ｂ(t) に含まれる上限周波数 ｆma＋ｆmb （時間領域） たたみ込み （周波数領域) 積 ＦＦＴを使って周波数領域で積として計算すれば、 計算量は Ｑ → ｌｏｇ2 Ｑ に比例して低減 例） 1000 → 10 ただし、信号に FFT データ長の遅延が発生する。 また、ＤＦＴスペクトル同士の積、および ディジタル信号同士の積には注意が必要である。 アナログ信号どうしの乗算結果と ディジタル信号どうしの乗算結果は必ずしも一致しない アナログの世界 ｘ（ｔ） × ｕ(t) ＝ ｙ(t) ＬＰＦ ｆｓ/2 ＬＰＦ ｆｓ/2 ＬＰＦ ｆｓ/2 A/D：ｆｓ A/D：ｆｓ A/D：ｆｓ ｘ（k） × ｕ(k) ＝ ｙ(k) ディジタルの世界 53 乗算結果不一致の具体例 何と一致するか？ アナログの世界 2kHz の 2kHz の 4kHz の 正弦波 × 正弦波 ＝ 正弦波 ＋直流 アナログの世界 2kHz の 2kHz の 4kHz の 正弦波 × 正弦波 ＝ 正弦波 ＋直流 3kHz ＬＰＦ LPFなし 3kHz ＬＰＦ 折り返し発生 直流 ｆｓ＝6kHz のA/D と一致 ｆｓ＝6kHz のA/D 2 ＝ 直流 ＝ ディジ タル内で 折り返しが 発生 2 ＝ ディジタルの世界 ディジタルの世界 折り返し発生の例 乗算による折り返し誤差の発生 120 ２乗信号 ｓ(k)×ｓ(k) 100 80 周波数 周波数 折り返し誤差は ＡＤするときに ｆs/2 以上の成分が 含まれていると発生する ことは知っている。 スイープ 正弦波 ｓ(k) 60 40 折り返し しかしそれだけでなく、 20 10 20 30 40 50 60 時間 時間 ディジタル信号どうしを乗算した場合にも、 折り返し誤差が発生する。 （非線形処理も） コンピュータ内で信号の乗算を行なうと 折り返しが発生するという例 対策は ？ ディジタル信号で 信号どうしの乗算を行う前には補間を！ 信号どうしの乗算 信号周波数の増加（最大２倍) 乗算と補間の効果の例 サンプ リング 周波数 ｆｓ サンプリング周波数の見なおし 補間 ２倍の補間を行って （サンプリング周波数を上げて） から乗算 2ｆｓ ｘ(k) の２乗信号 = 直流 ＋ ( ｆｓ のｃｏｓ波) ｘ(k) = (ｆｓ/2) の ｃｏｓ波 信号 1 0 1 k -1 0 1 1 0 -1 k k 0 補間しておけば、乗算しても大丈夫 k 54 アナログ信号のたたみ込み結果と ディジタルの信号のたたみ込み結果は一致する 周波数領域（DFT結果) の乗算にも同様の注意 例）パワースペクトルはスペクトルの２乗 アナログの世界 ｘ（ｔ） ＊ ｕ(t) ＝ ｙ(t) ＬＰＦ ｆｓ/2 ＬＰＦ ｆｓ/2 ＬＰＦ ｆｓ/2 A/D：ｆｓ A/D：ｆｓ A/D：ｆｓ ｘ（k） ＊ ｕ(k) ＝ ｙ(k) |Ｘ(p)|２ 補間なし スペクトル形状の 凸凹の頻度 は倍増 実際の パワースペクトル ２倍補間 ｆ ディジタルの世界 たたみ込みは安心して実行できる。 乗算は注意。 ◇最低でも２倍の補間は必要 （周波数軸上のサンプリング定理） ◇視覚的に谷などを明確にするためには多数倍の補間 周波数領域： ｚ変換とDFT 周波数領域（DFT結果) の積の解釈 周波数軸上の積 Ｚ 変換 ［時間］ ： 無限長 ［周波数］： 連続 周波数スペクトルの曲線に 高い周波数成分が発生する。 DFT ［時間］ ： 有限（周期） ［周波数］： 離散 補間の必要 周波数領域の乗算 周波数軸上の標本化周波数を増加 （＝ 補間）する必要性 アナログ 時間領域の直線たたみ込みに相当 Ｚ 変換 時間領域では、円状たたみ込み DFT アナログに合わせるには対策が必要 補間効果の例： 1025Hz 付近拡大 補間効果の例 1025Hz 正弦波の分析結果 （パワースペクトル） 0 補間なし -10 -20 -20 -30 -30 -40 -40 -50 -50 -60 -60 0 1000 0 2000 3000 周波数 [Hz] 4000 8倍補間 -10 -30 -30 -40 -40 -60 1000 2000 3000 周波数 [Hz] 4000 -60 2倍補間 16倍補間 -10 -15 0 1000 2000 3000 周波数 [Hz] 4000 -20 -25 16倍補間 -30 -35 700 -50 0 0 -5 0 -20 -50 2倍補間 -10 -20 補間なし 5 0 -10 0 1000 2000 3000 周波数 [Hz] 800 900 1000 周波数 [Hz] 1100 1200 4000 55 補間効果の例を見て 周波数軸上での補間の一方法 ・ 周波数領域のサンプリング定理を満たす ためには、２倍の補間は必須 周波数軸上での補間 ＝ 時間軸上のゼロ詰め ・ 視覚的に十分な理解を得るためには、 それ以上の補間が必要な場合もある ｘ0, ｘ１,…,ｘN-1 0, 0, …, 0 ・ 何倍の補間が必要か、の一般論は無いが、 目的とする（視覚的な）結果が得られるためには、 補間が必要であることは理解できる。 N (または、3N,7N,...) N 必要に応じて窓関数をかける 信号やスペクトルの乗算の際の注意点 （まとめ） あらかじめ 乗算による周波数成分の増加の 準備をしておく たたみ込み ( FIRフィルタリング）をＦＦＴで行う手順 N点の信号 ｘk ,･･･, ｘk+N-1 と N点のインパルス応答 ｂ0 ,･･･,ｂN-1 とのたたみ込み ２Ｎ点 ２Ｎ点 ｛ｘk ,･･･, ｘk+N-1,0,..0｝ → 補間、または ＦＦＴ Ｘ(p, k) N N ｛ｂ0 ,･･･,ｂN-1 ,0,..0 ｝ ＦＦＴ B(p) ２Ｎ点 DFTの結果（＝FFTの結果）の積は 円状たたみ込みだが、 ゼロを付加することで 直線状たたみ込みと同等の結果が 得られる ２Ｎ点 逆ＦＦＴ ２Ｎ点 ｙ(k) Y(p, k) 乗算を行う逆領域での零づめ 直線状たたみ込みをＦＦＴでできる トータル長を， 2N，4N, 8N,... などとして、ＦＦＴ ２Ｎ点 たたみ込み演算の 結果と一致 ＦＦＴを用いたフィルタの実行（１） ｂi ｘ(k) 切り出し（方形窓） Ｎ B(p) 0 ＦＦＴ 2Ｎ 零づめ Ｘ(p, k) 0 2Ｎ Y(p, k) ＝B(p)Ｘ(p, k) ＦＦＴ 逆ＦＦＴ 2Ｎ ｙ(k) 56 ＦＦＴを用いたフィルタの実行（２） フィルタの実行（まとめ） （overlap & add 法） ｘ(k) ◇ 差分方程式のプログラミング実行 Ｎ ◇ ＦＩＲ は、入力と係数のたたみ込みで、 F-1(B･Ｘ) ＦＦＴ を利用すれば、演算量が低減 2Ｎ ＋ ◇ 時間波形どうし、周波数スペクトル（特性） どうしの乗算は、補間をしてから ＋ ｙ(k) フィルタの設計 ４．ディジタルフィルタ ① ディジタルフィルタとは １） アナログフィルタの設計式を利用 （ＩＩＲフィルタ：双一次変換、など） ② フィルタの種類 と 伝達関数 ２） 設計ソフト、プログラムを利用 ③ 極と零によるフィルタの 特性解析 ３） ＤＦＴを使った簡易法（ＦＩＲ） ４） その他 （特別な注意が必要な、 固定小数点ＤＳＰ用の ＩＩＲ フィルタ設計 などは、専門書参照） ④ フィルタの実行 ⑤ フィルタの設計 簡易法 例）Ｎ＝８ Ｇ(p) G(0) G(1) G(2) G(3) Ｇ(5)～Ｇ(8) ＝ Ｇ(3)～Ｇ(1) G(4) ◇ 所望するフィルタの周波数特性を、ＤＦＴ の 値 Ｇ(p) として与える。これを 逆ＤＦＴして、 F -１ Ｇ(p) ｇ(k) 得られたインパルス応答 ｇ(k) を係数として ＦＩＲ フィルタリングを行う。 p Ｇ(p) ゲイン メリット： 自由な特性を 設計できる 勘違い例 (その１) 特性を設定 周波数 p この特性の逆ＤＦＴを用いれば， 理想的な低域フィルタになるか？ しかし、このままでは不適切な場合も・・ 57 勘違い例 (その２) 誤ったフィルタ動作 「すきま」の向こうに注意 サンプルのすきまに リップルが存在 k ＤＦＴ 阻止帯域の 減衰量も 実は小さい Ｇ(p) p 見えていないだけ 不要な成分を ゼロにする k p 補間すれば 見える p 不適切な理由 Ｎ k ＤＦＴ p 逆ＤＦＴ コメント この操作は、 周波数スペクトルの乗算 （＝ 時間領域のたたみ込み） なので、結果の時間波形は、 長さが２Ｎになるべきなのに、 なっていない → 時間領域で折り返し誤差 ＝回りこみが発生している。 Ｎ 短時間スペクトル分析に基づいた信号処理 例えば、雑音除去などの時変フィルタリングを 周波数領域の乗算でおこなう場合には、 この点に注意する必要がある。 必要なら、時間信号の適切な窓掛けや、 フィルタ特性は一度時間に戻して窓をかけてから 利用、などの対策を行う。 k p 逆ＤＦＴ その他の設計（実現？）法 １） 所望のインパルス応答を実現するフィルタ →ＦＩＲフィルタ の場合は、 所望のインパルス応答をフィルタ係数とすればよい 頭部伝達関数（ＨＲＴＦ）、 室内音響シミュレータ、など ２） その他の演算によって得られるフィルタ 逆フィルタなど （後述） フィルタの設計（まとめ） ◇ 通常は各種設計プログラムの利用が 便利 ◇ 簡易法を用いる場合は、いくつかの 注意が必要 58 信号処理対象のモデル化 講習概要 （発展編） ５．伝達関数による対象系のモデル化 ６．インパルス応答の測定法 信号処理を行うにあたって、 対象とする系をモデル化し 数式で表現することが必要 ７．逆フィルタ ８．適応フィルタとその応用 音響系を例に説明する ９．受音系における信号処理の応用 音の伝播の表現方法 波動方程式 （古典的表現） スピーカから出た音がマイクロホンにどう伝わるかを、 数式で表現する（モデル化） p ( x, y , z , t ) ： x,y,z 点における時刻 ｔ の音圧 ∂2 p ∂2 p ∂2 p ∂2 p = c2 ( 2 + 2 + 2 ) 2 ∂t ∂x ∂y ∂z ｃ： 音速 長所：任意の場所、任意の時刻における音圧を表現 短所：方程式が解けない （解は存在するが、簡単な関数で表現できない） → 有限要素法 マイクロホン スピーカ （ p.118） 伝達関数 （信号処理に適した表現） Ｘ(z) 豆知識：マイクロホン記号 Ｙ(z) Ｇ(z) 「点から点まで」 の伝達特性 Ｇ(z) で伝搬を表現 コンデンサ スピーカ ○ × Y(z)＝G(z)Ｘ(z) 長所： 複雑な系であっても、測定によって 求められる場合が多い マイクロホン 59 伝達関数を用いたモデル化の例 (1) 音源 ⇒ 発振器 空間 ⇒線形系 ⇒複数信号 は加算 対象となる系の入力と出力に注目して 伝達関数を用いた数式モデル（等価回路） として記述する。 入力 出力 ｘ(k) ｙ(k) 伝達関数を用いたモデル化の例 (2) ｙ(k) 楽器 （室内音響系） ～ 人 Ｇ(z) ～ Ｘ(z) Ｙ(z) （ p.121） 伝達関数を用いたモデル化の例 (3) S(z) ～ 声帯 声道 （周期（パルス）信号） （全極モデル） Ｘ2(z) （歌声） Ｇ2(z) Ｙ(z) マイク －信号分離問題－ Ｘ (z) 雑音が加わった音声 マイク 音声 Ｙ(z) Ｇ１(z) 伝達関数を用いたモデル化の利点の例 （音声生成系） Ｇ (z) Ｘ1(z) （音楽） 音声 雑音 N(z) 対象や目的に適合したモデルの選定が必要 音声に加わった雑音を除去し、 クリーンな音声を再現したい － 問題の定式化 － 伝達関数を用いたモデル化 問題の定式化 S(z) 音声 Ｙ1(z) ＧS1(z) ＧN1(z) マイク ＧS2(z) 雑音 N(z) ＧN2(z) Ｙ1(z)＝ＧS1(z)S(z)＋ＧN1(z)N(z) Ｙ2(z) Y1 ( z ) = GS1 ( z ) S ( z ) + G N 1 ( z ) N ( z ) Y2 ( z ) = GS 2 ( z ) S ( z ) + G N 2 ( z ) N ( z ) 行列表現 ⎡ Y1 ( z ) ⎤ ⎡ GS1 ( z ) G N 1 ( z ) ⎤ ⎡ S ( z ) ⎤ ⎢Y ( z ) ⎥ = ⎢G ( z ) G ( z ) ⎥ ⎢ N ( z ) ⎥ ⎦ ⎣ 2 ⎦ ⎣ S2 ⎦⎣ N2 Y ( z ) = G ( z ) X( z ) Ｙ1(z)、Ｙ2(z) は受音信号であるので既知。 G(z) が既知であれば、S(z), Ｎ(z) を未知数とした連立方程式 60 －信号分離問題の解－ G(z) の逆行列、G－１(z)を用いて解ける X ( z ) = G −1 ( z ) Y ( z ) −1 ⎡ S ( z ) ⎤ ⎡ GS 1 ( z ) G N 1 ( z ) ⎤ ⎡Y1 ( z ) ⎤ ⎢ N ( z ) ⎥ = ⎢G ( z ) G ( z ) ⎥ ⎢Y ( z ) ⎥ ⎣ ⎦ ⎣ S2 N2 ⎦ ⎣ 2 ⎦ ◇ モデル化を行うことで、複数のマイクロホンを 用いて （形式的には）複数の音源が分離できる ことがわかる。 ◇ ただし、実用的には、次の課題がある。 １) G(z) をどのようにして得るか（測るか）？ ２）G(z) の逆行列は実現可能か（安定か）？ 実用的な ｚ 変換多項式の計算 理論式では ｚ 多項式の計算であるが、実用的には コンピュータ内部のデータがインパルス応答数列 であるので、それに応じた計算を行う ◇ 乗算 Ｇ (z)･Ｘ(z) 多項式の積の実体は、多項式係数のたたみ込み なので、 時間領域と同様に、２つのインパルス 応答数列 （＝多項式の係数）のたたみ込みを 行う Ｇ (z)Ｘ(z) ＝ Z［ｇ(k)＊ｘ(k) ］ 伝達関数の代表的な求め方 インパルス応答 ｇ(k)＝｛ｇ(0), ｇ(1), ｇ(2), ｇ(3) ．．．｝ を測定できる場合は、これを ｚ変換して求める。 L −1 G ( z ) = ∑ g (k ) z − k k =0 = g (0) + g (1) z −1 + g (2) z −2 + g (3) z −3 + L ・ 理論的には上記の Ｇ(z) を用いる。 ・ 実用的には、コンピュータ内部で 伝達関数と等価な、インパルス応答の数列で保存 実用的な ｚ 変換領域（多項式）の計算 (2) ◇ 除算 Ｘ(z)/Ｇ (z) ・ 安定かどうかの判定が必要 ・ 安定なら、多項式の除算をして ｚ－１の級数として 展開し、有限で打ち切った（近似）展開係数 （多項式係数）をインパルス応答として保存 ・ 必要がなければ割り算を実行せず、分数のまま 保存しておくのが良いかもしれない （安定判別不要） ◇ 加算 インパルス応答（＝多項式係数） の和 → Ｇ (z)＋Ｘ(z) ＝ Z［ｇ(k)＋ｘ(k) ］ モデル化は必須か？ ◇ 近年、Blind （ブラインド）処理、ニューラルネット 適応フィルタ、など、必ずしもモデル化を行わな くても処理効果が得られる方法が提案されている。 ◇ しかし、ブラックボックス処理では、性能が 十分でないとき、理由がわからず、改善方針が たてられない。 （ 後述 ） 信号処理対象のモデル化 （まとめ） ◇ 伝達関数を用いて、対象（音響系）を モデル化 （等価回路表現） することがディジタル信号処理の第一歩 ◇ 伝達関数はインパルス応答より得られる 実用的演算は、インパルス応答数列 （＝多項式係数）に対して計算する ⇒ 常にモデル化を意識することが大切 ◇ 相対的な伝達関数で良いことも多い （マイクロホンアレー処理など） 61 ６．インパルス応答の測定法 講習概要 （発展編） ５．伝達関数による対象系のモデル化 ６．インパルス応答の測定法 ７．逆フィルタ ８．適応フィルタとその応用 ９．受音系における信号処理の応用 ◇ 測定法の概要 ◇ 測定原理 （測定信号を選べる場合） ◇ 測定信号 ・ TSP 信号： 信号と逆関数、演算量削減法など ・ M系列信号 ◇ 測定上の注意点（過大入力による非線形誤差） ◇ その他： 同期加算、など パルス法 インパルス応答測定法の概要 入力 ｓ(ｋ) 被測定系 ｇ(k) 出力 ｙ(ｋ) 被測定系 δ(k) ？ ｇ(k) 被測定系の入出力からインパルス応答 ｇ(k) を測定 １）入力 ｓ(k) を自由に設定できる場合 パルス法、ＴＳＰ法、掃引正弦波法、Ｍ系列法 ・ インパルスを入力とし、定義どおりに インパルス応答を測定する方法 （欠点）入力信号エネルギーが大きくできない ので、ＳＮ比が悪い ２）入力 ｓ(k) は制御できない場合 クロススペクトル法、適応フィルタの利用 （ p.158） （注意） 計算機から見た 測定結果 δ(k) D/A & LPF 被測定系 ？ LPF & A/D ｇ(k) 望ましい 測定入力信号の条件 １） 大きなエネルギを持つ信号 → ＳＮ比向上 ２） ただし、ある特定の時間にエネルギが集中 すると、系の非線形が発生するので、 ほぼ一定の振幅で持続する信号 計算機 ｇ(k) ３） 測定対象となる周波数成分を、 欠落無く含んでいる信号 ４） 扱いやすく、性質の良い信号 測定結果には、AD，DA の LPF の特性も含まれる 62 望ましい 測定入力信号の例 １） 掃引正弦波 （Swept Sine Signal） ・ 直線掃引 （周波数の遅れ時間が周波数に比例） TSP (Time Streched Pulse） ［鈴木、浅野］ ・ 対数掃引 （周波数の遅れ時間が周波数の対数に比例） ピンクTSP ［藤本］（Log-TSP、Log-SS) ・ 雑音スペクトルに応じた掃引 ［守谷、金田］ ２） Ｍ系列 （擬似ランダム雑音） ６．インパルス応答の測定法 ◇ 測定法の概要 ◇ 測定原理 （測定信号を選べる場合） ◇ 測定信号 ・ TSP 信号： 信号と逆関数、演算量削減法など ・ M系列信号 ◇ 測定上の注意点（過大入力による非線形誤差） ◇ その他： 同期加算、など 測定原理は共通 測定原理 円状たたみこみを前提とした逆関数 － 逆関数と円状たたみ込みの利用 － DFT ｓ(k) ① 測定信号 ｓ(k) と、その 逆関数 を ｓ-1(k) を作成 する。 ただし、ｓ(k)，ｓ-1(k) は、ともに 長さがＮ の信号であり、次の関係を満たす。 Ｓ-1(p) ＝ ｓ-1 (k) ＠ ｓ (k) ＝ δ(k) 1 Ｓ(p) 逆 DFT Ｓ-1(p) ただし、＠ は円状（巡回）たたみ込みを表す。 Ｓ(p) ｓ-1(k) この方法は、直線たたみ込みを前提とした 逆フィルタの設計方法としては、不適当である（後述） 測定原理 （つづき） 円状（巡回）たたみこみ （再） ② 測定信号 ｓ(k) を系に入力して、出力 ｙ(k) を得る。 ただし、ｓ(k) と 系のインパルス応答 ｇ(k) とが、 円状たたみ込み となるようにする。 ｙ(k)＝ ｓ(k) ＠ ｇ(k) ｓ(k) ｙ(k) ｇ(k) N個とN個のデータを畳み込んで、N個の結果を得る ｓ(k) を ２周期発生させ、 ２周期目に対応する出力 を ｙ(k) とする。 直線たたみ込みでは ２N-1 個となる Ｎ 直線たた み込み ＊ k Ｎ 円状たた み込み 2Ｎ-1 Ｎ ＝ Ｎ Ｎ ＠ k k k ＝ k k Ｎの長さからはみ出た分は回り込む 63 アナログ世界 （現実の系）で 円状たたみ込みを実行する方法 ｓ(k) ｓ(k) ・ ｓ(k) の長さＮ を インパルス応答長 より長く定める。 ・ ｓ(k) を ２周期発生 させ、２周期目に対 応する出力を ｙ(k) とする。 ＩＮ ｇ(k) OUT Ｎ Ｎ ｙ(k) は、ｓ(k) と ｇ(k) との 円状たたみ 込み ｙ(k) Ｎの長さからはみだした応答が、 回り込んでいるように見える ６．インパルス応答の測定法 ◇ 測定法の概要 ◇ 測定原理 （測定信号を選べる場合） ◇ 測定信号 ・ TSP 信号： 信号と逆関数、演算量削減法など ・ M系列信号 ◇ 測定上の注意点（過大入力による非線形誤差） ◇ その他： 同期加算、など 波形の引き伸ばしのイメージ 5 0 -5 0 500 1000 1500 2000 2500 3000 3500 4000 4500 4 3 2 1 0 -1 -2 -3 -4 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 ② 測定信号 ｓ(k) を系に入力して、出力 ｙ(k) を得る。 ただし、ｓ(k) と 系のインパルス応答 ｇ(k) とが、 円状たたみ込み となるようにする。 ｙ(k)＝ ｓ(k) ＠ ｇ(k) ｓ(k) ｙ(k) ｇ(k) ③ 得られた ｙ(k) と 逆関数 ｓ-1(k) と を計算機上で 円状にたたみこむ。 ｓ-1(k) ＠ ｙ(k)＝ｓ-1(k) ＠ ｓ(k) ＠ ｇ(k) ＝δ(k) ＠ ｇ(k) ＝ ｇ(k) となって、インパルス応答 ｇ(k) が得られる。 ＴＳＰ信号 (Time Stretched Pulse) ◇ インパルス信号δ(k)のＤＦＴ（フーリエ変換）Δ ( p )は、 Δ( p) = 1 for all p （定数倍は無視） 全ての周波数 p に対して等振幅・零位相 → エネルギーが時刻 0 に集中 ◇ ＴＳＰ信号は、 周波数に比例した遅延（または進み）を与えて、 信号を引き伸ばす （エネルギーを時間軸上に分散する） → 周波数 p の二乗に比例した位相成分 遅延の周波数表現 ∫ｘ(t) ｅ - j 2πｆ t dt ＝ X（ｆ） 10 -10 測定原理 （つづき） 時刻 0 に 集中していた エネルギーを 時間軸上に 引き伸ばす （分散させる） ∫ｘ(t-τ) ｅ - j 2πｆ t dt ＝ｅ - j 2πｆτ X（ｆ） τの遅延操作 →ｅ - j 2πｆτ ◇ 周波数 ｆ に比例した遅延 τ(f)＝ａ ｆ （ａ：定数） を与える周波数操作は、 2 → ｅ - j 2πａｆ ◇ インパルスの周波数スペクトル（＝１）に対して この周波数操作をした信号のスペクトルは、 2 → ｅ - j 2πａｆ 64 ＴＳＰ信号 ＴＳＰ信号 の 定義式 ◇ 周波数 p の二乗に比例した位相成分 ◇ Ｓ(p) を逆ＤＦＴ した信号 ｓ(k) がＴＳＰ信号 （TSP-up） ⎧ exp(− jap 2 ) Sup ( p) = ⎨ 2 ⎩exp{ ja ( N − p) } p = 0 ,1,L N/ 2 p= a = mπ(2 / Ｎ) N/ 2 + 1,L ,N − 1 6000 m : 整数 2 5000 ⎧ exp( jap 2 ) S dwn ( p ) = ⎨ 2 ⎩exp{− ja ( N − p ) } 周波数[Hz] （TSP-dwn） p = 0,1, L N / 2 p= N / 2 + 1,L , N − 1 （TSP-up） 4000 3000 2000 1000 0 0 0.2 0.4 0.6 0.8 1 1.2 時間[s] （ p.160～） ＴＳＰ信号 の 逆関数 Up & Down TSP （TSP-up） Ｓup（p）＝ ｅｘｐ（-ｊａp2） と （TSP-dwn） Ｓdwn（p）＝ ｅｘｐ（ｊａp2） を乗算すると、全ての周波数 p で １ となる。 （TSP-up） ｔ （TSP-dwn） δ(k) ｔ ※ （up が おすすめ ） （TSP-up） と （TSP-dwn） とは、 お互いに逆関数の関係 演算量の削減法 円状たたみ込みを利用したTSP測定 ＩＮ ＴＳＰ信号の長さＮをインパルス応答より長く定めて、 ＴＳＰ信号 ｓ(k) を２周期入力した時の出力 ｙ(k) の 第２周期目と、ｓ(k)の 逆関数とを、円状たたみ込み すれば、インパルス応答 ｇ(k) が得られる。 OUT ｓ(k) ｓ(k) ｇ(k) ｙ‘(k) ｙ(k) cut FFT 逆FFT Y(p) × たたみ込み （演算量） Ｎ×Ｎ 例）Ｎ=２１0 ：１００万回 ＦＦＴ ４Ｎ･log2Ｎ Ｓ-1(p) ４万回 65 たたみ込みと ＦＦＴ の補足 ◇ 円状たたみ込みを周波数（ＤＦＴ）領域で 行う場合には、零づめは不要 （長く延びた時間信号が回り込むのが 円状たたみ込み） 直線たたみ込みを周波数領域で行う 場合には、零づめ必要 （ 参照： ４章ディジタルフィルタ ④フィルタの実行） ＴＳＰ は逆関数の 直線たたみ込みが可能 鈴木、浅野らの設計による ＴＳＰ信号は、 波形の収束性が良いので、円状たたみ込みの 代わりに、直線たたみ込みを行っても、 誤差が小さい。 → ＴＳＰ 信号を２回鳴らす必要がない → あらかじめインパルス応答長を知る必要がない 欠点： FFT計算が使えない 参考文献： http://tosa.mri.co.jp/sounddb/tsp/index.htm Ｍ系列信号 Ｍ系列 （例：｛0, 1, 0, 1, 1, 0, … ｝） の 0 と 1 をそれぞれ +1 と -1 に対応させた信号。 ＤＡ後 ランダムな 1，-1 信号 ±1 しかとらない特 殊信号のようだが Ｍ系列 の作り方 (Maximum length sequence) ｎ 次のＭ系列とは、 ２ｎ－１の周期を持つ ０ と 1 のランダム系列 ◇作り方 ０か１の値を持つ ｎ段のシフトレジスタの 最終段を含む適切な段からの出力を 排他的論理和して、入力にフィードバック。 0 または 1 ｎ=4 Ｄ１ Ｄ２ Ｄ３ Ｄ４ ＬＰＦ後 ほぼ 白色雑音 Ｍ系列信号による測定 Ｍ系列信号 ｍ(k)の長さＮを インパルス応答より長く定めて、 Ｍ系列信号 ｍ(k) を２周期入力しした時の出力 ｙ(k) の第２周期目と、ｍ(k)の 逆関数ｍ(-k)とを、 円状たたみ込みすれば、インパルス応答 ｇ(k) が 得られる。 円状たたみ込み → アダマール変換を利用 排他的論理和 (mod２の和） ６．インパルス応答の測定法 ◇ 測定法の概要 ◇ 測定原理 （測定信号を選べる場合） ◇ 測定信号 ・ TSP 信号： 信号と逆関数、演算量削減法など ・ M系列信号 ◇ 測定上の注意点（過大入力による非線形誤差） ◇ その他： 同期加算、など 参考文献のリストが掲載されている文献： 金田、“Ｍ系列を用いたインパルス応答測定．．．” 日本音響学会誌、52巻10号、(1996.10) 66 測定上の注意点 非線形誤差の例（Ｍ系列法） ー 過大入力による非線形誤差 － ＳＮ比を稼ごうとして、入力が過大になると、 非線形誤差が発生する。（特にスピーカ） ・ Ｍ系列法では、パルス性雑音が発生し 特に時間後半で目立つ ・ TSP-upでは、非因果（主応答の前）の 波形として現れる 対策： スピーカレベルの低減、 測定系に含まれる非線型要素の点検 ＴＳＰ-up に現れる非線形誤差 TSP-up 4000 3500 3000 2500 Frequency 周波数 4000 3500 3000 Frequency 2500 3000 2000 1500 2000 1000 1000 0 0 500 0 0.1 0.2 0.3 Time 0.4 0.5 0.6 0 0.1 0.2 0.3 時間 0.4 0.5 Time 0.6 0.7 0.8 0.9 時間 TSP-up を使って、スペクトログラムで歪の有無を 検証することがおすすめ 測定信号レベルと測定誤差 雑音による誤差 測定誤差 周波数 TSP-dwn 5000 大 周波数 ＴＳＰ-up とＴＳＰ-down に現れる非線形誤差比較 6000 4000 Frequency ・ ＴＳＰ-dwn では、 インパルス応答の中に非線形誤差が含まれ、 目立たない。用途や誤差の大きさによっては、 良い場合もある。 ・ TSP-up では、 スペクトログラムで非線形誤差の存在がよく わかる。インパルス応答の負の時間方向に 非線形誤差が出現するので、信号レベルを 小さくして非線形誤差を小さくしないと不自然。 周波数 ＴＳＰと非線形誤差 2000 小 1500 非線形による誤差 最適レベル 1000 500 0 時間 0 0.1 0.2 0.3 0.4 0.5 Time 0.6 0.7 時間 0.8 0.9 小 測定信号レベル [dB] 大 TSP-dwn では、歪が見えないことがある。 67 ６．インパルス応答の測定法 ◇ 測定法の概要 ◇ 測定原理 （測定信号を選べる場合） ◇ 測定信号 ・ TSP 信号： 信号と逆関数、演算量削減法など ・ M系列信号 ◇ 測定上の注意点（過大入力による非線形誤差） ◇ その他： 同期加算、など 同期加算 ＳＮ比の改善に良く利用される手法 測定したインパルス応答の 時間軸を合わせて平均化 ｇ(k)＋ノイズ ランダムな雑音成分が減衰し Ｎ回の平均で、ＳＮ比を 10log10(Ｎ) [dB] 改善 時不変非線形による誤差には 効果がない 入力が制御できない場合の測定法 インパルス応答の測定（まとめ） ◇ クロススペクトル法（省略）は 経験とカンが必要 ・ 逆関数と円状たたみ込みを利用して インパルス応答を計測する。 ◇ 多少誤差があっても良いなら、 ・ 性質の良い測定信号として、 TSP や Ｍ系列が用いられる。 適応フィルタ （後述） の利用が簡便 （ p.159） ・ 過大入力による非線形誤差には 注意が必要である。 （同期加算でも低減できない） 今後の課題 ・ 各測定信号の長・短所の比較検討 ・ 目的に応じた測定信号、測定手順の 選択指針作成 ・ ノートパソコンなどで質のよいインパルス 応答測定ができるような環境作り、 ノウハウの蓄積、開示。 68 講習概要 （発展編） ７．逆フィルタ ５．伝達関数による対象系のモデル化 ◇ 逆フィルタとは ６．インパルス応答の測定法 ◇ 安定性 ７．逆フィルタ ◇ 近似的実現方法 ◇ 具体的計算方法 ８．適応フィルタとその応用 ◇ 多入出力形逆フィルタ ９．受音系における信号処理の応用 逆フィルタの例 (理想的残響除去) 逆フィルタ Ｘ(z) Ｇ(z) Ｘ(z) Ｙ(z) Ｙ(z) 逆フィルタ Ｈ(z) Ｘ(z) Ｙ(z)＝Ｇ(z)Ｘ(z) Ｇ(z) Ｙ(z)＝Ｇ(z)Ｘ(z) 伝達系で変形した信号 Y を元の信号 X に戻す Ｘ(z) 出力Ｙ(z) から 入力Ｘ(z) を復元するフィルタ 逆フィルタ Ｘ(z) Ｙ(z) Ｈ(z) Ｇ(z) H(z)＝１／Ｇ(z) Ｈ(z)Ｙ(z) ＝ Ｈ(z) Ｇ(z) Ｘ(z) ＝ Ｘ(z) Ｈ(z) Ｙ(z) ＝ Ｘ(z) （ p.124～） ｚ 変換 ・ DFT と逆フィルタ 振幅の逆フィルタ ・ （いわゆる）イコライザ ・ グラフィックイコライザ ｚ 変換 → 無限の信号長 を扱える Ｘ(z) 残響などの波形の変形は 復元できない Ｇ(z) Ｙ(z) DFT → 有限の信号長 インパルス応答（有限長）を DFTして求めた周波数特性 Ｇ(p) H(z)＝１／Ｇ(z) H(p)＝１／Ｇ(p) 正しい逆フィルタ： 問題は H(z) の安定性 誤った逆フィルタ： 有限長の逆数は無限長 であり、DFT では表せない 69 再度、ディジタル信号と ＤＦＴスペクトルの乗・除算に注意 DFT 周波数特性の逆数 １） 伝達系のインパルス応答 ｇ(k) を測定 ２） DFT により周波数特性 Ｇ(p)を計算 p：離散周波数 ３） Ｇ(p) の逆特性 Ｈ(p) ＝ １／Ｇ（p）を計算 ４） Ｈ(p) を逆ＤＦＴして、逆フィルタの係数 （＝インパルス応答） ｈ(k)を計算 ・ Ｇ(p) Ｈ(p) ＝ １ とはなる。 ・ しかし、ステップ ３） において、周波数領域の サンプリング定理が満たされなくなる。 ◇ 周波数成分の乗算・除算をすると 周波数特性の凹凸の頻度 （ｆ 軸を時間軸と考えた場合の周波数成分） が増大し、周波数軸の 標本化定理を満たさなくなる。 ◇ 乗算は、標本周波数を２倍にすれば良い。 ◇ 割り算はたちが悪い。 （かなりの高周波も発生するので、 ２倍程度ではだめ）。 DFT 周波数特性の逆数 １／Ｇ（p）の意味 Ｈ(p) ＝ ７．逆フィルタ 1 Ｇ（p） ・ ＤＦＴとは、Ｎ点の周波数特性と Ｎ点の時間波形の関係づけ ◇ 逆フィルタとは ・ 本来は無限長である逆フィルタのインパルス応答 を表すことができない。 ◇ 近似的実現方法 ・ H（p) を逆ＤＦＴしたＮ点の時間信号 ｈ(k) は、 無限長の応答をＮ個に圧縮したもの（時間軸での 折り返し誤差）。 → 誤った逆フィルタ ◇ 多入出力形逆フィルタ ◇ 安定性 ◇ 具体的計算方法 ・ 円状たたみこみに対しては逆フィルタとなる 逆フィルタの安定性 逆フィルタの簡単な例（伝達関数） ｘ(k) ｕ(k) Ｇ(z) の 逆フィルタ ：ＩＩＲ ｚ-1 は、不安定となる場合も多い。 u (k ) = x(k ) − cx(k − 1) 1 H ( z) = G( z) ｘ(k-1) -c･ｘ(k-1) -ｃ G( z) = 簡単な例を用いて、安定条件を 調べる。 U ( z) = 1 − cz −1 X ( z) この回路の伝達関数 （ p.125～） Z U ( z ) = (1 − cz −1 ) X ( z ) H ( z) = 1 1 = G ( z ) 1 − cz −1 この回路の逆フィルタの 伝達関数 70 練習： 逆 ｚ 変換 逆フィルタの簡単な例（回路導出） 1 1 H ( z) = = G ( z ) 1 − cz −1 逆フィルタ H(z) Y ( z) = のインパルス応答を計算してみる （ 時刻 k= 0,1,2,3 まで） (1 − cz −1 )Y ( z ) = U ( z ) H (z ) y ( k ) = u ( k ) + cy ( k − 1) Y ( z ) = U ( z ) + cz −1Y ( z ) ｚ 変換 インパルス応答 Z-1 伝達関数 ｙ(k) ｕ(k) 逆 ｚ 変換 逆 ｚ 変換 1 U ( z) 1 − cz −1 ｙ(k-1) ｃｙ(k-1) ・ 複素積分 ・ ｚ-1 で級数展開 ｚ-1 ｃ 逆フィルタの回路 FIR フィルタの逆フィルタ 逆フィルタの簡単な例（インパルス応答） F I R フィルタの ｙ(k) ｕ(k) ｙ(k-1) ｃｙ(k-1) ｚ-1 ｘ(k) ｕ(k) ｚ-1 逆フィルタは ｘ(k-1) -ｃ ｃ インパルス応答は ｕ(k) I I R フィルタ （無限長応答） ｃｙ(k-1) ｛1, c, c2, c3, …｝ インパルス応答は ｙ(k-1) ｚ-1 ｛1, c, c2, c3, …｝ 例 逆フィルタは不安定となる場合がある。 逆フィルタが安定な条件 最小位相だと安定な理由 Imag Q 線形系 Ｇ(z) が、最小位相系であること Ｇ(z) ＝ 最小位相系とは （定義）： G(z) の全ての極 ｐｉ と 零 ｑｉ の絶対値が１未満 | ｐｉ |＜１， | ｑｉ |＜１ i ＝ 1,2,3,... 複素平面上 単位円内 ｂ0Π(1- ｑi ｚ-1) ｉ =1 P 逆フィルタでは、 零と極が入れ替わる Π (1- ｐi ｚ-1) ｉ =1 Real −1 先ほどの例： G ( z ) = 1 − cz 零点は ｚ＝ｃ ｙ(k) ｃ ・ この逆フィルタは ＩＩＲフィルタ。 ・ |ｃ|≧１ の時、逆フィルタは不安定となる。 （発散する）。 ：極 ：零 -c･ｘ(k-1) 逆フィルタ P H(z) ＝ Π (1- ｐi ｚ-1) ｉ =1 よって、G(z) の零 qi が 単位円内に無いと 逆フィルタは不安定 Q ｂ0Π(1- ｑi ｚ-1) ｉ =1 G（ｚ）が最小位相系、が 逆フィルタの安定条件 （極 ｐｉ が | ｐｉ |＜１ は、線形系 G(z) の安定存在条件） 71 最小位相系でない、とは (1) 最小位相系でない、とは (2) 最小位相系でないのは、どのような時？ ① G(z) の零点 ｑｉ が | ｑｉ |＝１ の時。単位円上の零。 → 周波数特性上の零がある 逆特性は 無限大に発散 ある周波数でゲインが零 Ω c 1 0 １ |Ｇ(ｅｊΩ)| |G(ｅｊΩ)| ② G(z)の零点 ｑｉ が | ｑｉ |＞１ の 時、 信号エネルギーの遅れがある。 1 k 左の図の応答で、 |ｃ|＞１ | ｑｉ | ＞１ 信号エネルギーの 主要部が遅れる。 遅れを回復する事（＝時間を進める事）はできない。 Ω （ p.129） ７．逆フィルタ ◇ 逆フィルタとは ◇ 安定性 ◇ 近似的実現方法 ◇ 具体的計算方法 逆フィルタの近似的実現 ・ 最小位相系でない伝達系も多い ・ 非最小位相系の信号回復 の要求も多い → 安定性の保証されたＦＩＲフィルタで 逆フィルタを近似的に実現 ◇ 多入出力形逆フィルタ （ p.129） 代表的な近似方法 ◇ 最小２乗原理に基づく方法 （時間領域での設計法） ・ 行列計算 （測定したインパルス応答を用いて） ・ 適応フィルタの利用 （実環境での計算） （インパルス応答は未知でよい。 計算量が少ない） ◇ 直線たたみ込みを円状畳み込みで近似 （周波数領域の設計法） 最小２乗原理の考え方 原信号 δ(k) ＋ 回復 変形 Ｇ ｇ(k) Ｈ ｙ(k) 近似逆フィルタ ｅ(k) 誤差 原信号を良好に回復 → 逆フィルタの出力 ｙ(k) と原信号δ(k) との 誤差ｅ(k)のエネルギー J を最小にするような ＦＩＲフィルタ を求める。 J = ∑{δ(k) − y(k)}2 = ∑{δ(k) −∑g(k −i)h i }2 k k i g(k), h i ： G，H の インパルス応答 72 近似的逆フィルタの導出式 J = ∑ {δ(k ) − y (k )} 2 = ∑ {δ(k ) − ∑ g (k − i ) h i } 2 k k i この Ｊ を最小とするような ＦＩＲ型近似逆フィルタ Ｈ のフィルタ係数ｈｉ を求める。 Ｊ は、ｈｉ の２次関数であるので、 （解法１） 次の連立方程式を解いてｈｉ が求まる。 ∂J ∂J = 0 = 0, L , ∂ h Lh ∂ h1 ∂J = 0, ∂h0 近似的逆フィルタの問題点 Ｇ が非最小位相系である場合、 不安定な逆フィルタをそのまま近似しても 良い結果は得られない。 対策 (1) 雑音の付加 (2) 遅延の付加 ◇ ただし、実用的には、後で述べる行列算法や、 適応フィルタの利用（低演算量）で計算される 雑音の付加 遅延の付加 （Delayed Inverse) 遅延 δ(k) ＋ ＋ Ｇ ｇ(k) ｙ(k) Ｈ ＋ 誤差 δ(k-m) ｚ-ｍ ｅ(k) δ(k) Ｇ ｇ(k) Ｈ 回復 雑音 ｎ(k) ＋ ｙ(k) ｅ(k) 誤差 ・ Ｇ による遅れを補償する 不安定な部分を近似すれば、 雑音が増幅されて誤差が増大。 → 不安定部の近似を回避 （性能は雑音の大きさで調整） ・ 非最小位相系の逆フィルタは 因果だが不安定、または、 安定だが非因果な応答になる。 そのうち、安定だが非因果な成分を反映。 ・ 遅れても良いなら 大変有効 良好な近似逆フィルタを得るための ブロック図 遅延 ｇf(k-m) ｚ-ｍ δ(k) ＋ ＋ Ｇ ｇ(k) ｙ(k) Ｈ ＋ ｅ(k) 誤差 雑音 ｎ(k) ７．逆フィルタ ◇ 逆フィルタとは ◇ 安定性 ◇ 近似的実現方法 ◇ 具体的計算方法 ◇ 多入出力形逆フィルタ 上図において J = ∑ e (k ) 2 を最小化 k 73 行列算法 (たたみ込み行列 Ｇ) たたみ込み行列演算 Ｇ 0 ⎤ 0 0 L ⎡ g (0) ⎢ g (1) g ( 0 ) 0 M ⎥⎥ ⎢ ⎢ g (2) ⎥行で見ると g (1) g (0) 0 ⎢ ⎥インパルス g (2) g (1) g (0) ⎢ g (3) ⎥ 応答の逆順 ⎢ M ⎥ M M O ⎢ ⎥ − − ( 0 ) g ( N ) g ( N 1 ) g ( N 2 ) g ⎢ ⎥ ⎢ 0 g(N ) g ( N − 1) M ⎥ ⎢ ⎥ 0 g(N ) M ⎥ ⎢ 0 ⎢ 0 ⎥ 0 0 O ⎢ ⎥ ⎢⎣ g ( N )⎥⎦ たたみ込みを表す行列方程式 G と ｙ が与えられた 時、方程式を満たす ような未知数 ｈ を求 G ｈ める。 0 0 0 ⎤ ⎥ 0 g (0) ⎥ ⎡ h(0) ⎤ ⎥⎢ g (1) g (0) ⎥ ⎥ ⎢ h(1) ⎥ M O ⎥ ⎢ h(2) ⎥ g (0) ⎥ ⎢ ⎥ ⎥ M ⎥ g(N ) M ⎥⎢ ⎢⎣h( L)⎥⎦ ⎥ 0 g(N ) ⎥ 0 O g ( N )⎥⎦ ｙ(k)＝ｇ(k)＊ｈ(k) ｙ ＝ ⎡ y (0) ⎤ ⎡ g (0) ⎢ y (1) ⎥ ⎢ g (1) ⎢ ⎥ ⎢ ⎢ y (2) ⎥ ⎢ g (2) ⎢ ⎥ ⎢ ⎢ ⎥=⎢ M ⎢ ⎥ ⎢g(N ) M ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ 0 ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢⎣ y ( N + L)⎥⎦ ⎢⎣ 信号ベクトル ｈ = ［ｈ(0),ｈ(1),ｈ(2), ･･･, ｈ(L) ］ との積はたた み込みになっている Ｌ＋1 インパルス応答 ［ｇ(0), ｇ(1), ｇ(2), ｇ(3), ･･･, ｇ(N-1), ｇ(N)］ を 縦ベクトルとして、1段ずつずらして並べた行列 0 ⎡ g (0) ⎢ g (1) g (0) ⎢ ⎢ g (2) g (1) ⎢ M ⎢ M ⎢g(N ) ⎢ g(N ) ⎢ 0 ⎢ 0 ⎢ ⎢⎣ 0 ⎤ ⎥ Ｇ ⎥ ⎡ h(0) ⎤ ⎥⎢ ⎥ ⎥ ⎢ h(1) ⎥ O Ｎ＋1 ⎥ ⎢ h(2) ⎥ Ｌ＋1 ＋Ｌ g (0) ⎥ ⎢ ⎥ ⎥⎢ M ⎥ M ⎥ ⎢ ⎥ ⎥ ⎣ h( L ) ⎦ g(N ) ⎥ O g ( N )⎥⎦ 例えば 0 3行目＝ｇ(2)ｈ(0)＋ ｇ(1)ｈ(1)＋ｇ(0)ｈ(2) 行列は、縦長行列 0 0 g (0) 最小２乗原理に基づく計算法 δ ＝ G ｈ ただし、δ はδ関数 ［ 1, 0, 0, 0, ･･･ ］T T：転置 理想解として、 ｙ＝δ とおく。 しかし、 Ｇ が縦長行列なので、この方程式を満たす 解 ｈ は存在しないが、（δ－Ｇｈ）の２乗誤差 を最小にする ｈ は、 次式で求められる ( h = GT G ) −1 GT δ （ p.155） 不安定対策 １） 雑音の付加 ｇ(k) に 微小な（-40～-60ｄB程度、用途に依存） 雑音を付加して計算をする。 ２） 遅延の付加 目標出力 δ （または δ’）の時間をずらす （全長の 1/2 程度） ［ 1, 0, 0, ･･･］ → ［ 0, 0, ･･･, 0, 1, 0, ･･･, 0, 0 ］ 周波数領域の計算法 ◇ 基本的には、インパルス応答 ｇ の DFT Ｇ の逆数 1/Ｇ の 逆DFT で求めるが、 そのままでは時間軸での折り返し現象が σ おきて誤差が増大する。 ◇ 最も簡易な対策はＧの谷を埋めることである。 例えば、Ｇ(f) → Ｇ(f) ＋ σ σは、定数で、Ｇの最大値の-40～-60dB程度の値 Ω 詳しくは、例えば、 鈴木、浅野ほか：“音響系の伝達関数の模擬をめぐって（その１）” 音響学会誌、44巻12号、pp.936-942 （1988） 74 多入出力形逆フィルタ（ＭＩＮＴ） ７．逆フィルタ ｓ ◇ 逆フィルタとは ◇ 安定性 G1 H1 G2 H2 ｓ ＋ ◇ 近似的実現方法 G1やG2 が非最小位相であるため、個々の経路の 逆フィルタ Ｈｉ(z)＝１／Ｇｉ(z) （ ｉ=1,2 ） が不安定 であっても、経路を複数設ければ、 Ｈ1(z)Ｇ1(z)＋Ｈ2(z)Ｇ2(z)＝１ となる Ｈ1(z)、Ｈ2(z) は安定なフィルタとなる。 ◇ 具体的計算方法 ◇ 多入出力形逆フィルタ （ p.133） 多入出力形逆フィルタの 簡単な解釈 ｓ G1 H1 G2 H2 ＋ 音場制御における多入力逆フィルタ ｓ ｓ sss H G 受聴点 ｓ Ｇ1 に 非最小位相零があって、 逆特性が得られない周波数帯域は、 Ｇ2 の逆特性を利用する。また逆も同じ。 自由度を増やせば、無理なく逆フィルタ が得られる。 H1 G1 H2 G2 ｓ Ｈ1Ｇ1＋Ｈ2Ｇ2＝１ 部屋やスピーカの特性を キャンセル（逆）しないと 所望の音にはならない。 が、１チャンネルで無理 をするとうまく行かない。 ２チャンネルなら うまくいく。 MINT は チャンネル数を増すと 音場制御効果が上がる 理論的根拠を与える 逆フィルタのまとめ ◇ Ｇ(z) の特性の逆フィルタ 1/Ｇ(z)は、 Ｇ(z) が最小位相系でないと不安定になる ◇ 逆フィルタが不安定な時は、 以下の点を考慮してＦＩＲフィルタで近似。 １）系Ｇの谷をつぶしてから逆フィルタを計算 （雑音付加が相当） ２）回復の目標信号に遅延を許す （Delayed Inverse） ◇ 構成することが可能であれば 多入出力逆フィルタは有効 75 適応処理とは 講習概要 （発展編） ５．伝達関数による対象系のモデル化 状況に応じた適切な処理 ６．インパルス応答の測定法 非適応な処理 （旧型の処理形態） 状況によらず一定の処理 状況が変化しても同一の処理 ７．逆フィルタ ８．適応フィルタとその応用 ９．受音系における信号処理の応用 一次元音場*の騒音制御 適応処理の例 （* 音の進行方向が一次元） 簡単な適応処理 ダクト 掃除機や洗濯機の （強 ・ 中 ・ 弱） 高度な適応処理 （波形の制御） ＧＸ 消去点 Ｘ Ｎ マイク Ｇ ＷＸ スピーカ 騒音源 Ｗ 騒音制御、ほか ・ Ｗ＝－Ｇ とすれば、騒音は消去される ・ ＧＸ と ＷＸ の波形の一致が必要 ・ Ｇの変化に対応する必要がある 高度な適応処理 適応フィルタとは ｄ(k) 目標信号 適応フィルタ 最適なフィルタ特性を自動的に計算し、 実行するディジタルフィルタ ・ ＤＳＰ （Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ） の進歩に伴って実用化が進む ・ さまざまな適応信号処理 発展の原動力 入力 ｘ(k) 適応フィルタ 目的信号 出力 ＋ 希望信号 ｙ(k) － ｗ(k) 誤差 ｅ(k) ｘ と ｄ が与えられた時、ｄ に似た ｙ を合成して 誤差 ｅ のパワーを最小とする機能を持つフィルタ 自動的に誤差パワーを最小化する （ p.136） 76 一次元騒音制御における 適応フィルタの適用例 適応フィルタの代表的応用例 ダクト ｄ Ｎ １） 未知系の同定 （未知系の持つ入出力特性の推定） 消去点 Ｇ 騒音源 Ｘ ２） 逆フィルタの設計 適応フィ ルタ Ｗ ３） 予測 （未来の信号の予測） 誤差 ｅ 自動的に －Ｇ の特性を持つ 未知系の同定 入力 ｘ(k) 未知系 ？ ｄ(k) Room-A Room-B エコーキャンセラ 適応フィルタ ｗ(k) 同定の例 （音響エコーキャンセラ） 出力 ＋ ｙ(k) － 誤差 ｅ(k) ネット ワーク エコーキャンセラ 適応 フィルタ 伝達特性 － ＋ エコー 同じ入力ｘ(k) に対して、同じ出力を出すフィルタ係数、 エコー消去 ｗ(k)が求まれば、ｗ(k) が未知系の特性と考えられる。 エコー発生 インパルス応答の測定法としても利用される 逆フィルタの設計 予測 ｄ(k) ｄ(k) 現在 ｄ(k) 未知系 ？ ｘ(k) ｙ(k) ＋ 適応フィルタ － ｗ(k) 誤差 ｅ(k) ・ 未知系によって変形した原信号の回復を最適化 ・ 必要に応じて、遅延や雑音を付加する（→ 逆フィルタ） ｄ(k) 現在 ｘ(k) 遅延 過去 ｙ(k) ＋ 適応フィルタ － ｗ(k) 誤差 ｅ(k) ・ 過去の信号を用いて、現在の信号を 予測できるように学習させる ・ 学習後、ｘ(k) に現在の信号を入力（遅延０）すれば、 未来が予測できる 77 予測の応用 ｄ(k) 遅延 現在 適応フィルタによるハムの除去例 1 ｄ(k) ｘ(k) 適応フィルタ ｙ(k) －＋ ｗ(k) 過去 0.5 処理前 0 -0.5 -1 誤差 ｅ(k) 適応フィルタの内部構成 1 1.5 2 2.5 3 3.5 2 2.5 3 3.5 1 適応フィルタ ＯＮ 0.5 0 -0.5 -1 0 0.5 1 1.5 sec. フィルタ部（可変ＦＩＲフィルタ）の構成 出力 ｙ(k) フィルタ部 0.5 sec. 実際には、予測できるのは、ｄ(k) に含まれる周期成分 → ｙ(k) を出力 → ｄ(k) に含まれる周期成分の抽出 （雑音に埋もれた周期信号の抽出） → ｅ(k) を出力 → ｄ(k) に含まれる周期成分の除去 （ハムキャンセラ、など） 入力 ｘ(k) 0 フィルタ係数 入力 ｘ(k) ｘ(k-1) ｘ(k-2) ｘ(k-L+1) Ｚ-1 Ｚ-1 Ｚ-1 ｗ1(k) ｗ2(k) ｗ3(k) ｗL(k) 適応アルゴリズム 誤差 出力 ｙ(k) は、入力 ｘ(k-j+1) と、係数 ｗｊ(k) ｅ(k) 出力 ｙ(k) との積和 （たたみ込み） として得られる。 時間によって 変化する 誤差パワー ２次関数の特長 L ｊ=1 ｅ(k) ＝ ｄ(k)－ｙ(k) 適応フィルタの目標： 誤差パワー J＝Ｅ［ ｅ(k)2 ］ を最小化 ２次関数 係数 ｗ Ｊ ＝ Ｅ［｛ｄ(k)－ｙ(k) ｝2］ 誤差パワー 誤差パワー ｙ(k) ＝ Σｗｊ(k) ｘ(k-j+1) 多次関数 係数 ｗ 誤差パワーが係数の２次関数であると、．．． ・ 唯一の最適解が存在し、解析的に求めることが容易 ＝ Ｅ［｛ｄ(k) ーΣｗｊ(k) ｘ(k-j+1) ｝2］ 誤差パワー Ｊ は係数 ｗｊ(k) の２次関数である （ p.138） ・ 誤差パワーが小さくなるように係数を修正すれば、 探索的に求めることも容易 78 誤差曲面の等高線表示 （Ｌ＝2） 誤差曲面 （Ｌ＝2） 誤差パワー 係数 ｗ2 ｗ(k) ｗ(k+1) 係数 ｗ1 係数 ｗ2 最適フィルタ係数 係数 ｗ1 誤差パワーの最小値を与える係数が、最適フィルタ係数 この ｗ(k) を少しずつ修正していって、最適フィルタ係数 に近づけていく。 修正の方針 → 誤差パワーが最大減少する方向 （急斜面） 係数の修正方法（適応アルゴリズム） ある点での最大減少方向 - ∂ｅ2(k) ∂ｗ 係数 ｗ2 ｗ(k) ｗ(k+1) ＝ 2･ｅ(k)･ｘ(k) 学習同定法 （最も広く利用されている） ｗ(k+1)＝ｗ(k)＋ 係数 ｗ1 ｘ(k)＝［ｘ(k)，ｘ(k-1)］T （最も簡単な適応アルゴリズム） ｗ(k+1)Tｘ(k) ＝ｗ(k)Tｘ(k) ＋ 修正は、入力 ｘ と誤差 ｅ に基づいて行われる ｄ(k) ＋ 誤差ｅ(k) ＦＩＲフィルタ 出力 ｙ(k) 同定 ◇ ＦＩＲと系が同じ特性を持つということは、 系のインパルス応答 ≒ ＦＩＲ のインパルス応答 ◇ ＦＩＲ のインパルス応答 ＝ フィルタ係数 ∴ 最適なフィルタ係数は、系のインパルス応答の推定値 （ＴＳＰなど測定用入力が利用できない場合の推定法） （ p.141） 収束速度と定常誤差の関係 収束速度 誤差パワー － ｅ(k)･ｘ(k)Tｘ(k) ｘ(k)Tｘ(k) ＝ ｙ(k)＋ｅ(k) ＝ ｄ(k) μ：ステップサイズ 適応フィルタはインパルス応答の推定手段 ｅ(k)･ｘ(k) α：ステップサイズ (0<α≦1) α＝1， β＝0 のとき ｗ(k+1)＝ｗ(k)＋2μｅ(k)･ｘ(k) 系Ｇ α ｘ(k)Tｘ(k)＋β 最適フィルタ係数 フィルタ係数の修正 （LMS 法） 入力 ｘ(k) 時刻 ｋ におけるフィルタ係数 を ｗ(k) と表す。 ｗ(k)＝［ｗ1(k)，ｗ2(k)］T 定常誤差 （トレードオフ） ステップ 大 サイズ α 早いが誤差が大きい 遅いが誤差が小さい 小 修正回数 ｋ （ p.146） 79 入力信号による収束速度の差 入力信号と収束速度 有色スペクトルの信号は収束が遅い 誤差パワー ◇ 適応フィルタの収束の早さは、 入力信号の性質に依存 白色雑音 白色信号 → 速い 有色信号 → 遅い 音声 各種高速アルゴリズム → 主として信号の白色化 修正回数 ｋ （ p.146） 適応フィルタがうまく動作する例 （適応雑音抑圧処理） Ｓ Ｇ 騒音Ｎ Ｓ＋Ｇ･Ｎ ＋ － Ｗ･Ｎ 適応フィ ルタ Ｗ Ｎ Ｇ の特性を持って騒音Ｎを消去 Ｇ1 騒音Ｎ1 騒音Ｎ2 Ｇ2 Ｎ1＋Ｎ2 Ｓ Ｓ 騒音Ｎ1 騒音Ｎ2 誤差 （ p.148～） 複数の雑音源が存在する場合 Ｓ 誤差パワーが小さくならない （適応フィルタがうまく動作しない）例 Ｓ＋Ｇ･Ｎ1 ＋ － Ｇ Ｎ2 適応フィ ルタ Ｗ 誤差 モニタしているのは、消すべき騒音Ｎ1とは異なった 騒音Ｎ2 （騒音Ｎ2から騒音Ｎ1をつくることは不可能） モニタマイク Ｍ0 で受音する雑音に 伝達特性ＧM が付加される場合 Ｓ＋Ｇ1Ｎ1 ＋Ｇ2Ｎ2 ＋ － Ｓ Ｓ＋Ｇ1Ｎ1 ＋ － Ｇ1 適応フィ ルタ Ｗ 誤差 複数の騒音Ｎ1，Ｎ2 が複数の経路Ｇ1，Ｇ2で受音されて いる場合。適応フィルタは、Ｇ1となってＮ1を消去すること はできるが、その場合、Ｎ2は消去できない。 騒音Ｎ1 ＧM 適応フィ ＧMＮ1 ルタ Ｗ Ｍ0 誤差 雑音を消去するためには、適応フィルタの特性は、 Ｇ1/ＧM となる必要があり、これは逆フィルタとなるので、 不安定な場合は、雑音を消去できない 80 適応フィルタが良好に動作するためには･･･ 誤差パワーを小さくするフィルタが 理論的に存在すること が必要条件である （その特性は不明や時変でもよい） 動作不良の場合は理論的考察が必要 今後の展望 ・ 適応信号処理は高度な処理ゆえに 計算コストがネックであった。 ・ 近年のDSP技術の進展に伴って 処理コスト解決の見通しが得られた。 ・ 今後は各種の民生機器への導入が 急速に進展するものと考えられる。 存在が保証されていれば、 ｘ と ｄ を与るだけで、自動的に、 誤差パワーを小さくする特性を見出す。 81 講習概要 （発展編） 受音系における信号処理の応用 ５．伝達関数による対象系のモデル化 ６．インパルス応答の測定法 ◇ マイクロホンアレー＊）による指向性制御 （信号の遅延時間制御） ７．逆フィルタ ◇ 音の方向検出 （相互相関関数による信号の時間差検出） ８．適応フィルタとその応用 *) 複数のマイクロホンよりなる受音装置 ９．受音系における信号処理の応用 マイクロホンアレーによる指向性制御 加算形アレー（遅延和アレー）の原理 τ ２種類の制御方針 Ｍ1 １）加算形： 目的音を強調 遅延和アレー （超指向性アレー） 遅延 2τ ｄ θ ２）減算形： 不要音（雑音）を除去 適応型アレーなど 遅延 τ Ｍ2 ＋ Ｍ3 τ＝ d sinθ / c θ：目的方向 ｃ：音速 （ p.181） （ p.174） 遅延和アレー （超指向性アレー）の指向特性 遅延を変えれば 指向性を変えることができる 通話者 Ｍ3 Ｍ2 θ ｄ 鋭い指向性 を作る Ｍ1 τ 遅延 τ ＋ 雑音 遅延 2τ 82 ディジタル技術を利用すれば 簡単に任意の遅延を実現できる アナログ遅延器のインパルス応答 １） サンプリング周期 Ｔs の整数倍の遅延の付与 （遅延 τ＝ｎＴs の場合 ｎ：整数） 信号を、時間 τ だけ遅らせるアナログ遅延器 のインパルス応答 は、δ（ｔ –τ） 信号を ｎ サンプル分 シフトすれば（遅らせれば）よい δ（ｔ） 細かい遅延制御を行うには、 サンプリング周波数を上げて、 細かくサンプリングすれば良い 遅延器 τ δ（ｔ –τ） δ（ｔ –τ） に相当するディジタルインパルス応答 （＝ δ（ｔ –τ） を ＡＤ変換した信号） を持ったディジタルフィルタを実現すればよい。 ２） ディジタルフィルタを利用した任意遅延の付与 （ p.215） 遅延量のインパルス応答の離散化 δ（ｔ） 理想 ＬＰＦ A/D ｓin {πk } πk 遅延ディジタルフィルタのインパルス応答 ｈ(k) ＝ ［sinc 関数］ ｓin {π( k -τ/Ts)} π( k -τ/Ts) k=0 以外では 0 δ（ｔ –τ） 理想 ＬＰＦ A/D ｓin {π( k -τ/Ts)} π( k -τ/Ts) k： 離散時間（整数） Ts： 標本化周波数 sinc 関数を τ/Ts サンプル シフトした信号 ＦＩＲフィルタによる時間 τ の遅延器 τ/Ts ・ Ｎタップの FIR フィルタで実現するには、0 ～ N-1 の 範囲で ｈ(k) を 切り出して、ｈ(0)，･･･，ｈ(N-1) を フィルタ係数とする。 注） τが小さい場合には打ちきり誤差が大きい → 可能なら、k＜0 の部分も係数とする。 指向特性の制御 ディジタルフィルタの係数（プログラムの変数） を変えることで、遅延が自由に制御できる ＦＩＲ フィルタ ｘ(k) ｘ(k-τ/Ts) 0 k 0 → 指向特性の自由な制御 ｔ τ/Ts 標本化周期 Ts で DA すれば、 時間τ遅れた アナログ信号 が得られる 83 遅延時間τがわかれば、 音源方向θがわかる 音源方向検出 τ １） 指向性のビームをスキャンして、 最大パワーとなる方向を検出 Ｍ1 ２） 相互相関関数を利用した方向検出 ｘ1(k) ｄ θ ３） 多チャンネル相関行列を利用した 高分解能信号処理 （p.203） Ｍ2 ｘ2(k) τ＝ d sinθ / c θ：音源方向 ｃ：音速 θ＝ sin-1（ｃτ/ｄ） （ p.197） 相互相関関数 ２つの信号 ｘ1(k) と ｘ2(k) との 相互相関関数 ｘ1(k) ｔ φ12 (τ) = ∑ x1 (k +τ ) x2 (k ) ｘ2(k) 遅延時間の検出 ｘ1(k) ２つの信号、 k τ：離散時間 （整数） ｔ ｔ ｘ1(k) をτだけ ずらしながら ２つの波形の積の 総和を計算。 重なったところで最大 （積がすべて正となる） 最大値の位置に注意 相互相関関数 φ12 (τ ) τ 0 τｓ ここではなく ｘ2(k) ｔ ｘ1(k) と ｘ2(k) の ｔ 相互相関関数 φ12(τ) のピーク値を与える τの値が τｓ τｓを与える 音源追従 マイクロホンアレーシステム １） ２つのマイクロホンで得られた信号の 相関関数φ12(τ)を計算。 ２） 相関関数が最大値をとるτの値τsを求める。 こちら ３）音源方向θsを求める。 （方向検出方法は他の方法でも可） ４）θs方向に指向性ビームを形成 方向の推定精度をあげる場合、 必要に応じて補間を行う 84 ２種類の指向性制御方式の比較 減算形アレーによる雑音の消去 １） 加算形：目的音を強調 超指向性アレー ○処理が簡単 ×大規模 雑音 遅延器 D M1 y 1 d M2 ２）減算形：不要音（雑音）を除去 適応型アレー ○小規模 ×処理がやや複雑 d sinθN 適応フィルタを用いた遅延の推定 適応形マイクロホンアレーの構成 θN + τN 可変フィルタ 雑音 ｈ1 θN y y d M2 ｈM + マイクロホン アレー τN d sinθN ＋ ｈ2 適応 フィルタ M1 フィルタ 制御部 （ p.186～） 適応形マイクロホンアレーの実験例 雑音抑圧効果のデモ 0 （目的音） （不要音） (1) 音声 ＋ 雑音 ：無指向性マイク (2) 音声 ＋ 雑音 ：適応形アレー To desired sound -10 雑音源 マイクロホン (a) Relative responce (dB) 8.5 cm -20 To undesired sound -30 -40 0 1 2 Frequency (kHz) 壁 3 4 (3) 音声 ＋ 音声 ：無指向性マイク （天気予報） （ニュース） (4) 音声 ＋ 音声 ：適応形アレー マイクロホン アレーの位置 目的音源 (b) 85