目次 信号処理工学 Ⅰ 第 回 : ディジタルフィルタ 電気通信大学電子工学専攻電子知能システム学講座 ディジタルフィルタ ( 復習 ) FIR フィルタの補足 IIR フィルタの設計 IIR フィルタの実現 FIR フィルタと IIR フィルタの比較 最後の課題 長井隆行 ディジタルフィルタ ( 復習 ) 線形位相 FIR フィルタの補足 FIR フィルタ フィードバックがない インパルス応答が有限 伝達関数が分子だけ N n h( n) n IIR フィルタ フィードバックがある インパルス応答が無限 伝達関数に分母がある ak : 分子の次数 k N : 分子のフィルタ長 : 分母の次数 bk : 分母のフィルタ長 k h (n) : フィルタ係数 : 次数 ( の数 ) N : フィルタ長 ( の数 ) FIR フィルタ IIR フィルタ 線形位相 FIR フィルタの つのタイプ (No. のスライドも参照 ) タイプ タイプ タイプ タイプ タイプ 対称性 偶対称 奇対称 N 奇数 偶数 奇数 偶数 jω jω N / 周波数特性 n n n)co )/ π j ω )/ n π j ω N / β ( n)co n ω n β ( n ω 設計できないフィルタ なし ( 全て可能 ) HPF, BSF (LPF と BPF が可能 ) LPF, HPF, BSF (BPF のみ可能 ) LPF, BSF (HPF と BPF が可能 )
設計できないタイプのフィルタ 問題 () 前項の表で設計できないタイプのフィルタとはどういう意味か? ( 例 ) タイプ の線形位相 FIR フィルタでは HPF( ハイパスフィルタ ) が作れない タイプ 線形位相 FIR フィルタで作れるのが BPF だけであることを確かめよう π にしたいところだが H ( ω) N / β ( n)co n ω n N / H ( π ) β ( n)co n ω n 絶対に になってしまう H ( ω) )/ n タイプ 線形位相 FIR フィルタ HPF がだめな場合は BSF もだめ π π 6 バンドパスフィルタ (BPF) 解答例 () IIR フィルタの性質 BPF 以外のフィルタを作るためには H () ) / n H ( π ) n で ( に近い値 ) になる必要がある しかし実際には H (), H ( π ) )/ となってしまい 常に である 従って BPF 以外は作ることができない BPF であれば 問題ない フィードバックがある ( 伝達関数に分母がある ) インパルス応答が無限 全ての極が単位円内にある必要がある ( 設計時に安定性を考慮する必要がある ) 完全な線形位相特性が実現できない 近似的に満たす ( 満たすように設計する ) 設計や実装がFIRフィルタに比べ複雑 FIRフィルタに比べ少ない遅延器で急峻 ( 良好 ) な周波数振幅特性を得ることができる 7 8
9 IIR フィルタの設計法 IIR フィルタの設計は FIR フィルタに比べ難しい 振幅特性 位相特性 安定性を考慮する必要がある アナログフィルタ ( 安定性を満たすもの ) をディジタルフィルタに変換する よく使われる 基準ローパスフィルタ ( アナログフィルタ ) 双一次 - 変換 an H ( ) b N a b L a フィルタ係数を直接決定する L b c L K L d N cn d N N 所望の周波数特性と位相特性との誤差を最小とするように係数を決定する ( 最適化問題 ) 難しい c clar all アナログフィルタ バターワースフィルタ 次数が N で カットオフ周波数 (-db) が [rad/] の基準フィルタ 細かいことは ATLAB にまかせてしまう %[Z,P,K] buttap(n) は N 次の正規化されたプロトタイプ Buttrworth アナログ % ローパスフィルタの零点 極 およびゲインを出力します このフィルタは左半平面 % で 単位円周上に N 点の極をもち 零点をもちません N; % 次数 のバターワースフィルタを求める [Z,P,K] buttap(n); % [NU,EN] ptf(z,p,k ) は 伝達関数 % NU() % H() -------- % EN() % を作成します ベクトル Z は零点の位置 ベクトル P は極の位置 スカラ % K は利得です ベクトル NU と EN は それぞれ の降ベキ順に並べら % れた分子と分母の係数が出力されます [NU, EN] ptf(z,p,k); % ラプラス変換の周波数応答 W:/:; % 角周波数 [rad/] から [rad/] まで / きざみ Hfrq(NU,EN,W); % 振幅特性のプロット plot(w,ab(h)) p_.m バターワース特性 ( 基準ローパスフィルタ ) [ rad / ] アナログフィルタ 周波数 ( 振幅 ) 特性の変換 チェビシェフフィルタ clar all 次数が N で カットオフ周波数 (-db) が [rad/] の基準フィルタ 細かいことは ATLAB にまかせてしまう %[Z,P,K] chbap(n,rp) は 通過帯域に Rp db のリップルをもつ N 次の正規化され % たプロトタイプ Chbyhv I 型アナログローパスフィルタの零点 極 ゲインを出 % 力します Chbyhv I 型フィルタの遮断帯域応答は 可能な限り平坦になります N; % 次数 のチェビシェフフィルタを求める [Z,P,K] chbap(n,); % 次数は N でカットオフ周波数 [rad/](-db) で指定 % [NU,EN] ptf(z,p,k ) は 伝達関数 % NU() % H() -------- % EN() % を作成します ベクトル Z は零点の位置 ベクトル P は極の位置 スカラ % K は利得です ベクトル NU と EN は それぞれ の降ベキ順に並べら % れた分子と分母の係数が出力されます [NU, EN] ptf(z,p,k); % ラプラス変換の周波数応答 W:/:; % 角周波数 [rad/] から [rad/] まで / きざみ Hfrq(NU,EN,W); % 振幅特性のプロット plot(w,ab(h)) p_.m [ rad / ] チェビシェフ特性 ( 基準ローパスフィルタ ) 基準ローパスフィルタ H () の振幅特性を所望特性 H ˆ ( ) に変換する LPF LPF H () H ˆ ( ) H () ω c LPF BPF H () H ˆ ( ) ω ( ) H ωc [rad/] ( ) H ωb ω LPF HPF ω c LPF BSF H () H ˆ ( ) ω [rad/] ω ω ω ω ω ω ( ) H H ˆ ( ) ω ω ω ωω [rad/] b b [rad/] ω c ˆ ω b H ( ) H ω ω ω ω
周波数 ( 振幅 ) 特性の変換 双 次変換 LPF HPF H () H ˆ ( ) [rad/] clar all %HPF に変換した伝達関数 ( 手計算 ) NU[ ]; EN[ 6.8.96.9.8 ]; H ( ) ˆ H ( ) H.... 基準バターワースフィルタ (p_.mの結果) 6.8.96.9.8 アナログフィルタ ( 領域 ) からディジタルフィルタ ( 領域 ) へ変換 安定性を考慮 対 の対応 T 変換 ( 写像 ) 双 次 - 変換 % ラプラス変換の周波数応答 W:/:; % 角周波数 [rad/] から [rad/] まで / きざみ Hfrq(NU,EN,W); % 振幅特性のプロット plot(w,ab(h)) %ATLAB の関数 lphp が使える %NU[ ]; %EN[.... ]; %[NUT,ENT] lphp(nu,en,); p_.m H() の根 ( 極 ) が 平面の左半平面になければならない ( 必要十分条件 ) 平面 平面 H() の極が 平面の単位円内になければならない ( 必要十分条件 ) clar all 双 次変換の例 p_.m の例 (HPF) アナログカットオフ周波数 [rad/] /π[h].8[h] サンプリング周波数を [H] として変換.8π/.[rad/] 近辺がカットオフ周波数になるはず %ATLAB の関数 lphp が使える ( 周波数変換 ) NU[ ]; EN[.... ]; [NUT,ENT] lphp(nu,en,); ˆ H ( ) H 6.8 ˆ.8.9.8.78.8.96.78..9.9..8.8. IIR フィルタ設計まとめ 今までの手順を順番に行う 所望の特性 ( 仕様 ) 次数 カットオフなど アナログ基準ローパスフィルタ 周波数変換 LPF or BPF or HPF or BSF カットオフ周波数分子の次数 N 分母の次数 ( サンプリング周波数 ) バターワース or チェビシェフ所望の次数のものを ATLAB で作る LPF LPF or BPF or HPF or BSF ω ( ) H c tc. % 双 次変換 [NUd,ENd] bilinar(nut,ent,); % サンプリングレートは H % 周波数特性の計算 ( 域 ) [H,W]frq(NUd,ENd,); figur plot(w,ab(h)) % 振幅特性のプロット p_.m 6 双 次変換 ディジタルフィルタ T
IIR フィルタ設計の例 ATLAB を使う 所望特性 ( 仕様 ) サンプリング周波数 :8k[H] フィルタタイプ:BPF カットオフ周波数: f k[ H] f k[ H] 次数 :( 基本アナログフィルタ ) 基本アナログフィルタ: バターワースフィルタ次数 :( ディジタルフィルタ ) ω π [ rad / ] ω π [ rad / ] clar all IIR フィルタ設計の例続き N; % 次数 のバターワースフィルタを求める [Z,P,K] buttap(n); [NU, EN] ptf(z,p,k); % 伝達関数に変換 % 周波数特性の変換 %ATLAB の関数 lpbp が使える omga_**pi; % カットオフ ( 左側 ) k*pi[rad/] omga_**pi; % カットオフ ( 右側 ) k*pi[rad/] omga_b(omga_-omga_); omga_qrt(omga_*omga_); [NUT,ENT] lpbp(nu,en,omga_,omga_b); % 双 次変換 [NUd,ENd] bilinar(nut,ent,8); % サンプリングレートは 8kH 7 ω H ˆ ( ) ω ωω ω ω ω ω b [rad/] % 周波数特性の計算 ( 域 ) [H,W]frq(NUd,ENd,,8); figur plot(w,ab(h)) % 振幅特性のプロット figur plot(w,unwrap(angl(h))) % 位相特性のプロット 8 p_.m 9 問題 () 次のような仕様の IIR フィルタを作ってみよう ローパスフィルタ 次数 N サンプリング周波数 [H] カットオフ周波数 [H] 基本アナログバターワースフィルタ H ( ) IIR フィルタの実現 フィルタ ( システム ) の一般形 IIR フィルタ a k k b k k H H num dno ( ) ( ) H num () ( ) X ( ) X ( ) H 直接型 I H dno ( ) num dno No.9 p. システムの縦属
IIR フィルタの実現続き H num ( ) は入れ替えが可能 H dno ( ) dno H num ( ) H H dno ( ) num () 問題 () 直接形 II の伝達関数が直接型 I と同じになることを確かめてみよう a a a b b v(nt ) a b b a a 直接型 I 直接型 II X ( ) X ( ) H num ( ) dno 遅延器を共通化 直接型 II 差分方程式をたてて 変換する 伝達関数 ( 直接型 IIではv(nT ) を使った差分方程式をつたてる ) FIR フィルタと IIR フィルタの比較 FIR フィルタと IIR フィルタの比較 ( ローパスフィルタ ) FIRフィルタ (N: 遅延器 )Rmアルゴリズムを使用 IIRフィルタ (N: 遅延器 ) 双 次変換を使用 線形位相特性ではない完全な線形位相特性 第 回課題 音声のノイズ除去 * 定常なノイズの混合した音声信号 (ampl.wav~ ampl.wav) のノイズを除去する * ノイズのみの信号 (noi.wav) からノイズの性質をある程度知ることができる * 基本的にどのような方法を用いてもよい * 信号は全て 6kH サンプリングである * 分からない人は次の手順に従って処理してみる FIR IIR () ノイズのみのサンプル (noi.wav) を読み込み FFT することでノイズの周波数帯域を調べる () ノイズの周波数帯域を阻止域とするディジタルフィルタを設計する () ノイズの混合した音声を設計したディジタルフィルタでフィルタリングする
具体的な手順 具体的な手順続き より具体的な手順 ( 必ずしもこれに従う必要はない ) () ノイズのみのサンプル (noi.wav) を読み込み FFT し 振幅特性 ( 絶対値 ) をプロットする (noi.wav が サンプルなので FFT も ポイントで行う ) () プロットした図を見て 振幅がほぼゼロとなる部分とゼロではない部分に分割する ( ナイキスト周波数までを考えればよく 境目はそれほど厳密でなくてもよい ) () 分割した境目の周波数 [H] を計算する ( 周波数と ポイント FFT のインデックスとの関係は?) () 計算した周波数を角周波数に直し ( π) さらにサンプリング周波数 (6) で割る (T とした場合の正規化角周波数 [rad/c] が得られる ) 次のページに続く 6 () 計算した正規化周波数の値を使って FIR フィルタを設計する ( 窓関数法 p_.m でも Rm 法 p_6.m でもよい フィルタを設計する関数では ナイキスト周波数を として周波数を指定することに注意が必要 ) ノイズがゼロに近いところを通過域 ゼロでないところを阻止域とすればよい フィルタ長は 程度でよい (6) 設計したフィルタの周波数特性をプロットする 意図した周波数特性になっているかどうかチェックする ( ノイズの周波数帯域を除去するような特性になっているか?) (7) ノイズの混合した音声信号 ( ampl.wav~ampl.wav のうち一つを適当に選ぶ ) を読み込み プロットしてみる (8) 読み込んだ音声信号を設計したディジタルフィルタに入力し ( たたみ込みの計算 ) 出力された信号をプロットする (9) プロットした結果を見てノイズが除去されていることを確かめる できれば音を聞くとよい もし 結果があまりよくなければ () のフィルタ設計のパラメータを変えて何度か試してみる 提出方法 提出は前回同様メールで送ってください 締め切り :7 月 7 日 ( 金 ) : プログラムと結果をレポートとしてまとめてメールで提出 pdf 形式にする ( 体裁は自由 ) pdf にできなければ Sword でもよい ヒントプログラム :kadai_hint.m p@appl..uc.ac.jp 7 メールの件名に学籍番号 - 課題番号を半角で書いてください ( 例 ) 今回の課題を 6 の人が送る場合 6-