論 文 43.60.+d
2
チャンネル近接マイクロホンアレイによる風雑音の低減
*坂 田 直 人
∗1村 上 哲 郎
∗2中 島 弘 史
∗2中 臺 一 博
∗3 ［要旨］ 風雑音は一般的に非定常な雑音であり，信号の波形レベルでの相関をもとにした処理については あまり行われていない。本論文では 2 チャンネルを近接させたマイクロホンを用いて，各チャンネルで相関 のある風雑音の収録を行い，相関の分析・風雑音の低減の二つの実験を行った。振幅・パワー・複素信号のそ れぞれについてコヒーレンス関数により相関を分析した結果，どの項目についても 125 Hz 以下で 0.3∼0.8 の相関が確認された。その相関を利用して 2 種類の線形ビームフォーマにより風雑音の低減を行い，125 Hz 以下で 3∼10 dB 程度のパワーの低減が確認された。また，従来法（パワースペクトルサブトラクション）と 提案法とでカートシス比を比較し，提案法は従来法に比べて音質の点で優位であることが確認された。 キーワード 風雑音，雑音低減，マイクロホンアレイ，ビームフォーミングWind-induced noise, Noise reduction, Microphone array, Beamforming
1. ま え が き
1.1 は じ め に 風雑音はマイクロホンに風が吹き込むことにより発 生する雑音である。屋外での音声収録において，風雑音 は音質低下の原因となっており，その改善が望まれてい る。風雑音は時間領域・空間領域で高い非定常性を持っ ていることから，パワースペクトルサブトラクション （PSS）[1]やポストウィーナフィルタリング等の非線形 処理などによる雑音除去[2–4]が行われてきた。しか し，非線形処理は非線形ひずみを引き起こし，また，音 源の位相情報を回復不能にするなどの問題がある。そ れに対し，線形処理は音質の点から非線形処理よりも 優れているが，風雑音は空間的な相関を持たないことか らあまり適用されていない。また，風雑音の低減には， 一般的にはある程度の大きさのあるウィンドスクリー ンを用いるが，その大きさにより概観を損なったり，設 置スペースの問題により使用できない場合があり，小型 な機材による風雑音を低減可能な技術が求められる。 本研究では，小型な機材を用い，従来行われてこな かった線形処理による風雑音の低減を行い，非線形処 理よりも音質劣化の少ない低減を目的とした。小型簡 易マイクロホンアレイ（2チャンネル近接マイクロホ∗Wind-induced noise reduction using a small two-channel microphone array,
by Naoto Sakata, Tetsuro Murakami, Hirofumi Nakajima and Kazuhiro Nakadai.
∗1クラリオン ∗2工学院大学
∗3ホンダ・リサーチ・インスティチュート・ジャパン （問合先：坂田直人 e-mail: Naoto Sakata@clarion.
co.jp） （2015 年 12 月 14 日受付，2016 年 6 月 6 日採録決定） ンアレイ）を用いて風雑音を収録することで，通常は 相関のないチャンネル間の風雑音信号に相関を持たせ た。そして，非線形処理よりも音質の点で有利である 死角形成型の線形ビームフォーマにより風雑音の低減 を行い，SN比改善量と発生した非線形ひずみ量を従 来法と比較した。 1.2 過去の研究との比較 風雑音については過去に様々な研究が行われている。 風雑音の基本的性質は[5] で明らかにされている。ま た，[6]は二つのマイクロホンで風雑音を収録したと きの風雑音信号のチャンネル間の相関の近似式を求め ており，2種類の風向きについて相関が近似されてい る。一つは風の流れがマイクロホンに対面している場 合（Crosswind）で，もう一つはマイクロホンの並び 方に沿っている場合（Downwind）である（図–1）。そ れぞれの相関の近似式は，exp(−7x)，exp(−3.2x)で あった。ここでxは，相関をとる周波数の波長λに対 するマイクロホンの間隔dの比である（x = d/λ）。こ こでλ = v/f であり，vは風速（m/s）を，f は計算す る相関の周波数を表す。この近似式は，野外で大きな マイクロホンアレイを用いた実験により実験値を近似 する式として得られた。風雑音の低減についても様々 な試みが行われている。[7]では三つのマイクロホンを 音響抵抗体の中に入れた装置を用いている。[2]では2 チャンネルのマイクロホンを用いて，周波数領域で最 適なゲインを適応させている。[3]では，3チャンネル のマイクロホンアレイを用い，コヒーレンス関数に基づ きノイズのパワースペクトルから目的信号のパワース ペクトルを推定している。[4]では空間的に相関のある ノイズに対する雑音低減の手法として，ポストフィルタ を用いている。[2–4]はスペクトルのゲイン調整等の非
図–1 マイクロホンと風向きの関係 図–2 2 チャンネル近接マイクロホンアレイ 図–3 TCM-370 の周波数特性 線形処理による手法であり，音質の劣化が問題となる。 本研究では，非線形処理による音質劣化の問題を解 決するため，線形処理による風雑音の低減を目指す。 通常，風雑音は空間的な相関を持たないが，通常よりも 小型のマイクロホンアレイを用いることでチャンネル 間に相関のある風雑音を収録することができる。それ により，従来あまり行われてこなかった，死角形成型 の線形ビームフォーマにより線形処理が可能となった。
2. 2 チャンネル近接マイクロホンアレイ
2チャンネル近接マイクロホンアレイは，収録される 風雑音に相関を持たせるための小型のマイクロホンア レイである（図–2(a)）。素子の直径が4 mmの小型無 指向性マイクロホン「AV-LEADER TCM-370」（図– 2(b)）を二つ直線状に並べて構成される。各マイクロ ホンとの距離は，適宜変更して用いる。以降，Lチャ ンネル及びRチャンネルを単にL及びRと記述する。 「AV-LEADER TCM-370」の周波数特性を図–3に示 す。本研究では，この機材を用いて風雑音を収録し，分 析及び風雑音の低減を行う。3. 風雑音のコヒーレンス
3.1 実 験 環 境 2チャンネル近接マイクロホンアレイにより線形フィ 図–4 コヒーレンス計測の実験環境 表–1 コヒーレンス計測の収録条件 パラメータ 設定 扇風機との距離 158 cm 風向 Crosswind，Downwind チャンネル間の距離：d 10，20，30，50，100 mm 風速：v 2.2，1.7，1.1 m/s ICレコーダのゲイン設定 30 サンプリング周波数 44.1 kHz 収録時間 30 s ルタによる風雑音の低減が可能か確認するため，風雑 音のチャンネル間の相関をコヒーレンス関数により 分析した。風雑音の収録は一般室内で行った。実験時 の環境を図–4に示す。主な使用機材は2チャンネル 近接マイクロホンアレイ（図–2），扇風機「NEOVE FTS30-T12」，ICレコーダ「TASCAM DR-05」，風 速計「CELESTRON CE48021」の四つである。収録 時の機材の設定，風向，風速等の収録条件を表–1に示 す。風速はマイクロホン地点の風速を測定した。主要 なパラメータはチャンネル間の距離，風速，風向の三 つである。 3.2 コヒーレンス関数による解析 本実験では複素・振幅・パワーの 3 種類のコヒー レ ン ス を 解 析 す る 。一 般 的 な コ ヒ ー レ ン ス 関 数 γ(Al(ω), Bl(ω))は次式で定義される。 γ(Al(ω), Bl(ω)) = A ∗ l(ω)Bl(ω) |Al(ω)|2 |Bl(ω)|2 (1) ここでAl(ω)とBl(ω)は複素スペクトルであり，ωは離 散角周波数，lはフレームのインデックス（= 1, 2, . . . , L） を表す。Al(ω)は Al(ω) のl についての平均値（= L l=1Al(ω)/L）を表す。A∗l(ω)はAl(ω)の複素共役 である。複素コヒーレンスγc(ω)は次式で計算される。 γc(ω) = γ(X1,l(ω), X2,l(ω)) (2) ここでX1,l(ω)とX2,l(ω)は，チャンネル1, 2のそれぞれに対応した入力信号の複素スペクトルである。こ こで，チャンネル1はLチャンネル，チャンネル2は Rチャンネルに対応している。X1,l(ω)は次式で計算 される。 X1,l(ω) = N−1 n=0 x1(n + Sl)w(n)e−jωn/fs (3) ここでj は虚数単位，w(n)は窓関数，x1(n)はチャン ネル1の入力信号，nは離散時間，N はフレーム長か つFFT長，Sはフレームシフト長，fs はサンプリン グ周波数である。X2,l(ω)はx2(n)に対して同様に処 理することで得られる。実験ではN = 212，S = N/2 とし，w(n) をハニング窓とした。振幅コヒーレンス γa(ω)は次式で定義される。 γa(ω) = γ(G(|X1,l(ω)|), G(|X2,l(ω)|)) (4) G(As,l(ω))はAs,l(ω)のl についての分散を維持し たまま平均値を0にする関数である（G(As,l(ω)) = As,l(ω)− As,l(ω)）。s はチャンネル番号である。ま た，パワーコヒーレンスγp(ω)は次式で定義される。 γp(ω) = γ(G(|X1,l(ω)|2), G(|X2,l(ω)|2)) (5) 複素コヒーレンスは遅延時間の推定やビームフォーミ ングの雑音低減などの線形処理による能力の評価に用 いる。振幅・パワーコヒーレンスはスペクトルサブト ラクションなどの非線形処理の能力の評価に用いる。 3.3 実 験 結 果 図–5 はv = 2.2 m/s，Downwind，d = 10 mm で の3種類のコヒーレンスの値を比較したグラフである。 複素コヒーレンスについてはγc(ω)の絶対値を表示し ている。125 Hz以下の低周波帯域で複素コヒーレンス を含む各コヒーレンスの値がおよそ 0.3∼0.8 となっ た。風雑音はこの帯域において線形相関があることが 分かる。125∼500 Hzの帯域では，パワーコヒーレン スや振幅コヒーレンスが約0.4∼0.6 となったのに対 し，複素コヒーレンスの値は0.1と小さい値となった。 これはこの帯域ではチャンネル間に位相の相関はほと んどないが，パワーや振幅だけで見れば両チャンネル で同様なパワーの発生傾向があることを示しており， 片方のチャンネルを参照したPSS 等の手法がこの帯 域での風雑音低減について有効であることを示唆して いる。1∼8 kHz では各コヒーレンスの値が非常に高 くなっているが，これについては風速ごとでの複素コ ヒーレンスの実験結果において述べる。 図–6はv = 2.2 m/s，Downwindでのマイクロホン の間隔ごとでの複素コヒーレンスの値を示したグラフ である。マイクロホンの間隔が10 mm のとき，低周 波帯域で0.7以上の高い相関が得られた。マイクロホ 図–5 3 種類のコヒーレンス関数による相関 図–6 マイクロホン間隔ごとの複素相関 図–7 風向ごとの複素相関 図–8 風速ごとの複素相関 ンの間隔が小さい方がより高い相関を得られることが 分かる。図–7はv = 2.2 m/s，d = 10 mm での実測 による風向ごとの複素コヒーレンスの値と，それぞれ の理論的な近似式を比較したグラフである。実測にお いては，風向きがDownwindの場合に低周波帯域で高 い相関が得られた。また，[6]で報告された相関の近似 式との大きなずれはなかった。 図–8はd = 10 mm，Downwind での風速別のコ
図–9 複素相関の時間遷移 ヒーレンスの値である。風速が速くなるにつれて，1∼ 4 kHz付近の相関値が，より高周波側にずれている。収 録した風雑音の試聴時，風量が上がるにつれてファン の回転による風切り音が高い音に移行していく様子が 聴き取れたことから，この帯域での相関はファンの風 切り音によるものであり，本研究で取り扱う風雑音の 相関ではないと考えられる。また，最も遅い風速であ るv = 1.1 m/sにおいて相関が他と比べ若干低くなっ ていることから，風速が速いほうが高い相関が得られ ると考えられる。 図–9にv = 2.2 m/s，d = 10 mm，Downwindでの 複素コヒーレンスの時間による変化を示す。計算にお いては，各フレームについて，そのフレームを含んだ 直後30フレームでコヒーレンスを計算した。色の濃 度がコヒーレンスの値を示しており，明るい箇所ほど コヒーレンスの値が高い。低周波に着目するため，0∼ 200 Hzを拡大した。ここで，特に50 Hz以下で0.8程 度の相関が各フレームにおいて定常的に確認でき，低 周波帯域で安定した線形処理の効果が期待できる。 以上の結果より，2チャンネル近接マイクロホンア レイを用いて風雑音を収録した場合，125 Hz以下の低 周波帯域においてチャンネル間の風雑音信号に複素相 関が得られるため，それを利用した線形処理が可能で あることが分かる。
4. 線形ビームフォーマによる風雑音の低減
4.1 風雑音混入収録信号のモデル 第3章にて，2チャンネル近接マイクロホンにて風 雑音の線形処理が可能であることが示された。本節で は，本研究での線形処理の基盤となる，風雑音の混入 した収録信号のモデルを説明する。 風雑音が音響信号のように伝搬すると仮定し，風雑 音の伝達関数を音響伝達関数と同様に定義する。風雑 音以外の雑音がないと仮定し，ただ一つの音源x(n)を 風雑音と共に収録した収録信号y(n) は次式によりモ デル化される。 y(n) = x(n)∗ h(n) + q(n) ∗ k(n) (6) ここで∗ は畳み込み演算子，nは離散時間，h(n) は 再生・収録系のインパルス応答である。q(n)は風雑音 源で，再生・収録系における音源に相当する。k(n)は 風雑音の伝達特性を表した時間領域の信号であり，再 生・収録系におけるインパルス応答に相当する。式(6) をフーリエ変換することで，周波数領域でも時間領域 と同様な関係が次式のように得られる。 Y (ω) = X(ω)H(ω) + Q(ω)K(ω) (7) 式(6)，(7)のモデルは，実際には風雑音について相関 のある低周波帯域において考えられるモデルであるが， フィルタ設計時には全帯域において成り立つものとし た理論的なフィルタの設計を行う。 4.2 処理システム 本研究における風雑音低減の処理システムについて 説明する。Lチャンネル及びRチャンネルの収録信号 l(n)，r(n)をフレーム分割し，各フレームの時間信号 lm(n)，rm(n)を得る。lm(n)，rm(n)について窓関数 w(n) をかけて高速フーリエ変換（FFT）を行い，複 素スペクトルLm(ω)，Rm(ω)を求める。Lm(ω)は次 式で表される。 Lm(ω) = FFT{lm(n)w(n), N} (8) N はフレーム長もといフーリエ変換長，FFT{x(n), N} は x(n) に つ い て N 点 で FFT を 行 う 関 数 で あ る。Rm(ω)もLm(ω)と同様に求められる。Lm(ω)， Rm(ω)に対してフィルタGL(ω)，GR(ω) を乗算し， 処理後の出力Zm(ω)を次式により得る。 Zm(ω) = GL(ω)Lm(ω) + GR(ω)Rm(ω) (9) ここでGL(ω)，GR(ω) は，風雑音が除去されるよう に求める。 フィルタGL(ω)，GR(ω)を逆フーリエ変換（IFFT） することにより，時間領域の信号に対してフレーム 分割せずにフィルタを適用することを考える。フィル タGL(ω)，GR(ω)の時間領域での表現であるgL(ω)， gR(ω)を次式により得る。 gL(n) = IFFT{GL(ω), N} (10) gR(n) = IFFT{GR(ω), N} (11) IFFT{X(ω), N}はX(ω)についてN点の逆高速フー リエ変換を行う関数である。得られたgL(n)，gR(n) をフレーム分割していない時間領域の信号l(n)，r(n) に畳み込むことで，次式のようにフレーム分割なしで フィルタを適用できる。 z(n) = gL(n)∗ l(n) + gR(n)∗ r(n) (12)ここで，z(n)は処理後の出力信号である。本研究にお いては，式（12）により，線形ビームフォーマを適用す る。本研究でのフィルタの設計については，GL(ω) = 1 としてGR(ω)のみ調整する「線形減算ビームフォー マ（LS-BF）」と風雑音の伝達関数，所望信号の伝達関 数を利用し，所望信号を維持したまま風雑音を除去す る「音源分離ビームフォーマ（SS-BF）」の2種類の方 法でフィルタを作成する。 4.3 線形減算ビームフォーマ（LS-BF） GL(ω) = 1とし，風雑音のみを入力したときに出力 が0となるようにGR(ω)を求める。Lm(ω)とRm(ω) を風雑音のみの信号とする。このとき，L(ω)とR(ω) は式(7) のモデルに当てはめると次式のように表さ れる。 L(ω) = Q(ω)KL(ω) (13) R(ω) = Q(ω)KR(ω) (14) KL(ω)，KR(ω)はL及びRの風雑音の伝達関数であ る。式(9)に式(13)，(14)を当てはめると処理後の出 力Z(ω)は次式で表される。 Z(ω) = 1・Q(ω)KL(ω) + GR(ω)Q(ω)KR(ω) (15) Z(ω) = 0としてGR(ω)について整理すると，フィル タGR(ω)は次式により得られる。 GR(ω) =−KL(ω) KR(ω) (16) 4.4 音源分離ビームフォーマ（SS-BF） 風雑音のL，Rへの伝達関数をKL(ω)，KR(ω)，所望 信号の伝達関数をHL(ω)，HR(ω)とし，L(ω)とR(ω) を所望信号と風雑音の混入した信号とする。L(ω)， R(ω)は式(7)に当てはめると次式のように表される。 L(ω) = X(ω)HL(ω) + Q(ω)KL(ω) (17) R(ω) = X(ω)HR(ω) + Q(ω)KR(ω) (18) 式(9)にこれらの条件を当てはめて整理すると，処理 後の出力Z(ω)は次式となる。 Z(ω) = X(ω) (GL(ω)HL(ω) + GR(ω)HR(ω)) +Q(ω) (GL(ω)KL(ω) + GR(ω)KR(ω)) (19) ここで，所望信号がLで収録される風雑音のない信号 X(ω)HL(ω)とし，風雑音を表す信号をRで収録され るQ(ω)KR(ω)と考えて次式のように変形する。 Z(ω) = X(ω)HL(ω) GL(ω) + GR(ω)HR(ω) HL(ω) +Q(ω)KR(ω) GL(ω)KL(ω) KR(ω)+ GR(ω) (20) 所望信号を保持し，風雑音の出力が0になるようなフィ ルタGL(ω)，GR(ω)を考えると，式(20)よりGL(ω)， GR(ω)は次の条件を満たす必要がある。 GL(ω) + GR(ω)HR(ω) HL(ω) = 1 (21) GL(ω)KL(ω) KR(ω)+ GR(ω) = 0 (22) HR(ω)/HL(ω) = HLR(ω)，KL(ω)/KR(ω) = KRL(ω) として式(21)，式(22)を連立して解くとフィルタが 次式のように得られる。 GL(ω) = −1 KRL(ω)HLR(ω)− 1 (23) GR(ω) = KRL(ω) KRL(ω)HLR(ω)− 1 (24) HLR(ω)はL を基準としたRの音響相対伝達関数， KRL(ω)はRを基準としたLの風雑音の相対伝達関 数である。 4.5 実 験 所望信号を含んだ風雑音を収録し，理想伝達関数，実 測伝達関数それぞれの場合についてフィルタを作成し てLS-BF及びSS-BFによる風雑音の低減効果を従来 法（PSS [1]）と比較して分析した。低減効果の分析に ついては，雑音抑圧量の評価として元信号及び処理信 号に対してパワースペクトル，帯域別SN比を計算し， 処理による非線形ひずみの発生量の評価としてカート シス比[8]を計算し，結果の比較を行った。 4.5.1 音声収録環境 風雑音の収録は工学院大学八王子校舎の9号館無響 室で行った。実験室内の様子を図–10に示す。主な使 用機材は3章と同様であるが，風の発生についてはサー キュレータ「D2 DAS-KJ191」を，音声の収録におい てはAD/DA変換器「M-Audio Fast Track Pro」と ノートPC「Lenovo ThinkPad L520」を用いた。収 録時の環境等の条件を表–2に示す。 機材の配置について図–11に示す。所望信号はスピー カ「BOSE 101-MM」図–12，パワーアンプ「YAMAHA MX-1」により再生した。所望信号はRWC研究用音 楽データベース[9]のポピュラ音楽No. 8を 10 秒程 度のモノラル音源に編集したものを用いた。再生時の サンプリング周波数は44.1 kHzとした。
図–10 無響室 表–2 風雑音低減実験での収録条件 パラメータ 設定 扇風機との距離 200 cm 風向 Downwind チャンネル間の距離：d 4，10，20 mm 風速：v 2.6，1.8 m/s サンプリング周波数 44.1 kHz 収録時間 10 s 図–11 収録時の機材配置 図–12 BOSE 101-MM 4.5.2 理想相対伝達関数 理想相対伝達関数は，信号がその伝達過程で遅延し かしないと仮定して設計した伝達関数であり，信号の 伝達関数の測定が困難な場合に用いることができる。 フィルタの設計にはHLR(ω)及びKRL(ω) の二つのみ 求めればよい。所望信号及び風雑音の理想相対伝達関 数HLR(ω)，KRL(ω) は次式のように計算した。 HLR(ω) = e−jωts (25) KRL(ω) = e−jωtu (26) ts はLを基準としたRの所望信号の遅延時間，tu は Rを基準としたLの風雑音の遅延時間である。ts は音 速（340 m/s），マイクロホンの間隔d，風向の関係か ら計算して求めた。tu は，複素コヒーレンス関数（式 （2））の位相情報から信号の遅延時間を推定することに より計算した。複素コヒーレンス関数の位相情報は， 計算に用いた二つの信号の周波数成分ごとの平均的な 位相差を表す。その位相差についてアンラッピング処 理を施し，周波数に伴う位相の変化に連続性を持たせ た位相差関数φ(ω) の傾きから二つの信号間の遅延時 間を推定できる。φ(ω)は次式により計算される。 φ(ω) = unwrap{arg {γc(ω)}} (27) ここでarg{} は複素数の偏角を計算する関数であり， unwrap{}は位相のアンラッピング処理を行う関数で ある。推定遅延時間τˆはコスト関数J (τ )を最小にす るτ として得られる。J (τ ) はτ を変数としたφ(ω) と直線(−ωτ) の誤差の2乗和であり，誤差最小のと きに最良の近似直線の係数としてτˆが得られる。J (τ ) は次式により定義される。 J (τ ) = P p=1 |φ(ωp)− (−ωpτ )|2 (28) ここでωp はインデックスp に対応する離散角周波数 であり，P はφ(ω) と(−ωτ) の近似に用いる上限の 離散角周波数に対応するインデックスを表す。これは φ(ω) の(−ωτ)による線形近似であるため，τˆは擬似 逆行列を用いて次式のように計算できる。 ˆ τ =−ω+φ (29) ここでω は0からωP までの離散角周波数の列ベク トルである。φは位相差を表す列ベクトルであり，そ の要素はω に対応した値(φ = φ(ω))となっている。 また，ω+ はω の擬似逆行列である。 4.5.3 音響相対伝達関数の推定 LとRの音響伝達関数を測定し，計算により音響 相対伝達関数HLR(ω)を推定した。音響伝達関数の測 定はTSP（Time-Stretched Pulse）法を用いて行っ た。TSPの応答に対して，10回の加算平均処理を行 い，HL(ω)，HR(ω)を測定した。TSPの信号長は0.2 秒とし，振幅を0.5とした。推定したHL(ω)，HR(ω) からHLR(ω)を次式により推定した。 HLR(ω) = HR(ω)HL(ω) |HL(ω)|2+ max|HL(ω)|2× rh (30) ここでrh は発散防止係数である。
4.5.4 風雑音の相対伝達関数の推定 風雑音の伝達関数は風雑音源が計測できないため， KRL(ω)を直接求めた。フレーム分割した風雑音信号 を用いて，RからLへの平均的な伝達特性を計算した。 Lm(ω)，Rm(ω)を風雑音のみの信号とする。Lm(ω)， Rm(ω)の離散角周波数ω について各フレームでの複 素振幅を縦に並べた列ベクトルをLω，Rω とする。ω についてのRを基準としたLの伝達特性をKRL(ω) で表すと，理想的にはLω は次式のように表される。 Lω = KRL(ω)Rω (31) これをRω の疑似逆行列R+ω によりKRL(ω)につい て解くと次式となる。 KRL(ω) =R+ωLω (32) 得られたKRL(ω)は，KRL(ω)Rm(ω)のLm(ω)との 誤差の2乗和を最小にする。 4.5.5 フィルタの発散防止係数 式 (23)，(24) に よ り 定 義 さ れ た SS-BF は ， KRL(ωα)HLR(ωα) = 1 となるω = ωα でフィルタ が発散するので，それを防ぐための発散防止処理を行 う。本実験においてはSS-BFでは次式のように発散 防止処理を行い，フィルタを適用した。 GL(ω) =− d(ω) ∗ |d(ω)|2+max|d(ω)|2×rg (33) GR(ω) = KRL(ω)d(ω) ∗ |d(ω)|2 +max|d(ω)|2×rg (34) ここでrg は発散防止係数である。また，d(ω)は次式 のように定義される。 d(ω) = HLR(ω)KRL(ω)− 1 (35) 4.5.6 評 価 方 法 処理結果のパワースペクトル・帯域別SN比により， 低減性能を分析する。信号x(n) のパワースペクトル Px(ω)は，x(n)のフレーム分割信号xm(n)に対して 式(8)と同様な処理で求めた複素スペクトルXm(ω) から次式により求められる。 Px(ω) = 1 M M m=1 |Xm(ω)|2 (36) ここでM はフレームの総数である。各種音源につい てパワースペクトルを求め低減性能を周波数領域で確 認した。また，下限周波数f1 から上限周波数f2 にお ける帯域のSN比の計算については，信号のパワース ペクトルとf1，f2 に対応するその離散角周波数ωpの インデックスp1，p2（ωp1 = 2πf1，ωp2 = 2πf2）に より次式で求めた。 SNR(p1, p2) = 10 log10 p2 p=p1Ps(ωp) p2 p=p1{Ps+n(ωp)− Ps(ωp)} (37) こ こ で Ps(ω) は 所 望 信 号 の パ ワ ー ス ペ ク ト ル ， Ps+n(ω) は 所 望 信 号 と 共 に 風 雑 音 を 含 む 信 号 の パ ワ ー ス ペ ク ト ル で あ る 。式（37）に お い て ， p2 p=p1{Ps+n(ωp)− Ps(ωp)} が0以下になる場合が あるが，その場合は便宜的にSNR(p1, p2) = 0として， 計測不能として扱った。 カートシス比は，非線形ひずみの発生量としての評 価に用いられる指標である[8]。カートシス比κratioは 処理前と処理後のカートシスκori，κprocにより求めら れる。それぞれのカートシスκは次式により求めた。 κ = M m=1 P p=1 Cmp4 MP M m=1 P p=1 C2 mp MP 2 (38) ここでCmpは，カートシスを計算するパワードメイン 時間周波数領域信号であり，m，pはそれぞれフレーム， 離散角周波数のインデックスであり，M，P はそれぞ れフレーム及び離散角周波数インデックスの最大値で ある。実際の計算時には，式(8)と同様な処理で求め た複素スペクトルの絶対値を2乗して求めた。式(38)
より得られたκori，κprocから，カートシス比κratioを
次式により計算した。 κratio = κproc κori (39) κratio が大きい値を示すほど，発生した非線形ひずみ 量が大きく，ミュージカルノイズによる音質劣化が疑 われる。本実験においては，κori を風雑音低減前の信 号（風雑音＋所望信号）のカートシス，κproc を各種手 法による処理後の信号のカートシスとしてカートシス 比を計算した。 4.5.7 実験時パラメータ 実験時のフレーム長などのパラメータについて説明 する。低減処理は，風速2.6 m/sの環境で収録した音 源について行った。フレーム長は 2,048，シフト長は 1,024，FFT長はフレーム長と同じとし，窓関数はハ ニング窓を用いた。これらのパラメータは，PSSの実 行条件においても同様である。また，発散防止係数は rh = 10−4，rg= 10−2 とした。理想相対伝達関数は， Downwindについてはts= 0，tu =−1.9×10−3とし， Crosswindについてはts=−1.1765 × 10−5，tu= 0 として計算した。SN比は全帯域（0∼22,050 Hz）と，
風雑音が主な帯域（10∼500 Hz）と音声信号が主な帯域 （500∼10,000 Hz）の3種類について計算した。PSS の実行においては，先頭20フレームの風雑音のみの 区間により雑音パワーを推定し，低減処理を行った。 具体的には，パワードメインでスペクトル減算を行い， 観測信号の位相を付与した上で時間領域信号に戻した ものを再度短時間フーリエ変換してパワースペクトル などの分析を行った。また，PSSの実行の際，パワー 減算により負の値となる周波数成分については0を代 入した。 4.6 実験結果と考察 4.6.1 理想伝達関数 理想伝達関数による風雑音の低減処理の結果を示す。 図–13は，風向Cr（Crosswind）の風雑音に対して 理想伝達関数を用いたときの処理結果のパワースペク トルである。図中の“Original”は，処理する前のLの 信号のパワースペクトルである。横軸は周波数を，縦軸 は各周波数成分の相対パワーレベルを表す。PSSでは 500 Hz以下で 8 dB程度風雑音のパワーが低減した。 また，1,000 Hz以上の帯域では，広い範囲でパワーが 上昇し，Originalからのスペクトルの概形の崩れが見 られた。これは，PSSにおいてフレームの接合の際に 発生する，不連続成分による影響であると考えられる。 LS-BFでは全体的に，Originalよりもわずかにパワー が上昇し，雑音レベルの改善が見られなかった。これ は，理論伝達関数の設計誤差によるものであると考え られる。SS-BFでは500 Hz以下の低周波で15 dB程 度の低減が確認できた。また，1,000∼8,000 Hzの音 声帯域では，パワーの減衰が確認された。これは，伝 達関数の測定誤差によるものであると考えられる。 図–14は，風向 Dw（Downwind）の風雑音に対し て理想伝達関数を用いたときの処理結果のパワースペ クトルである。PSSはCrの場合と同様な形状となっ た。LS-BFでは低減が確認されなかった。これは理想 伝達関数の設計誤差によるものであると考えられる。 SS-BFでは 50∼1,000 Hzの周波数帯で Original か ら5 dB程度低減した。 表–3に理想伝達関数による各風向，手法でのSN比 の改善量を示す。表中の“帯域”はSN比を計算した 帯域であり，“LS-BF”，“SS-BF”，“PSS”は，各手法 についてのSN比の改善量を示す。LS-BFではどの場 合においても，SN比は改善できなかった。SS-BFで は，特にCrの場合において全帯域で13 dB程度SN 比が改善した。PSSでは，10∼500 Hzの風雑音帯域 で6 dB程度の改善が確認された。風向Crにおいて， SS-BFのSN比改善量がPSSの場合よりも 6 dB程 度上回った。 図–13 処理後のパワースペクトル（理想，Cr） 図–14 処理後のパワースペクトル（理想，Dw） 表–3 SN 比改善量（理想伝達関数） 風向 帯域 LS-BF SS-BF PSS Cr 0∼22,050 Hz −2.7 dB 12.9 dB 6.4 dB Cr 10∼500 Hz −2.8 dB 13.3 dB 6.5 dB Cr 500∼10,000 Hz −2.5 dB 3.3 dB −1.8 dB Dw 0∼22,050 Hz −0.8 dB 1.6 dB 7.7 dB Dw 10∼500 Hz −0.8 dB 1.5 dB 7.8 dB Dw 500∼10,000 Hz −1.8 dB −1.1 dB 1.8 dB 図–15 処理後のパワースペクトル（実測，Cr） 4.6.2 実測伝達関数 実測伝達関数による風雑音の低減処理の結果を示す。 図–15は，風向Crの風雑音に対して実測伝達関数 を用いたときの処理結果のパワースペクトルである。 250∼1,000 Hzの帯域で，LS-BF，SS-BFともに処理 前のスペクトルからほとんど変化がなかった。LS-BF では125 Hz以下でパワーが5 dB程度低減された。ま た，2,000∼4,000 Hzの音声帯域でパワーが10 dB程 度減衰した。SS-BFでは250 Hz以下でパワーが3 dB 程度増加した。図–16は，風向Dwの風雑音に対して 実測伝達関数を用いたときの処理結果のパワースペク
図–16 処理後のパワースペクトル（実測，Dw） 表–4 SN 比改善量（実測伝達関数） 風向 帯域 LS-BF SS-BF PSS Cr 0∼22,050 Hz 2.9 dB −2.3 dB 6.4 dB Cr 10∼500 Hz 3.1 dB −2.2 dB 6.5 dB Cr 500∼10,000 Hz 0.3 dB 0.2 dB −1.8 dB Dw 0∼22,050 Hz 2.2 dB −1.0 dB 7.7 dB Dw 10∼500 Hz 2.3 dB −1.1 dB 7.8 dB Dw 500∼10,000 Hz 0.4 dB 1.4 dB 1.8 dB トルである。LS-BF では，250∼1,000 Hz の帯域で 処理前のスペクトルからほとんど変化がなかった。ま た，風向Cr と同様に，125 Hz以下でパワーが 5 dB 程度低減され，2,000∼4,000 Hzの音声帯域でパワー が10 dB程度減衰した。SS-BFでは 62.5 Hz以下で パワーが3 dB程度増加した。また，125∼250 Hzの 帯域で，3 dB程度のパワーの低減が確認された。 表–4に実測伝達関数による各風向，手法でのSN比の 改善量を示す。0∼22,050 Hz，10∼500 Hzの帯域では Cr，DwともにLS-BFでのみ改善が確認され，その値 は3 dB程度となった。SS-BFではどちらの帯域につ いてもSN比が 2 dB程度減少した。500∼10,000 Hz の帯域では意味のある改善量は得られなかったが，各 手法でSN比が減少することはなかった。風向Dwで は，おおよその場合においてLS-BFのほうがSN比を 改善しているが，500∼10,000 Hzの音声帯域において はSS-BFによる改善量の方が大きな値を示した。0∼ 22,050 Hz，10∼500 Hzの帯域ではSS-BFではどちら の帯域についてもSN比が2 dB程度減少した。500∼ 10,000 Hzの帯域では SS-BF においてSN比の改善 が確認され，その値は1.5 dBとなった。 全体として，PSSによる処理ではスペクトルの形状 を崩す傾向が見られたが，LS-BF，SS-BFでは形状を 維持して処理された。表–5は各手法の処理信号につい て計算されたカートシス比である。PSSのカートシス 比がCrで2.09，Dwで4.38となっているのに対し， LS-BF，SS-BFにおいてはCr，Dwでの最大のカー トシス比の値はそれぞれ1.51，2.07であり，PSSに よる処理よりも非線形ひずみの発生量を低く抑えられ 表–5 各手法のカートシス比 風向 伝達関数 LS-BF SS-BF PSS Cr 理想 1.51 0.26 2.09 Cr 実測 0.35 0.76 2.09 Dw 理想 0.87 1.69 4.38 Dw 実測 0.96 2.07 4.38 図–17 処理前の音源のスペクトログラム ていることが確認された。これより，LS-BF，SS-BF はPSSよりも音質の点で優位であることが分かった。 4.6.3 人の音声信号に対する影響 提案手法が音声信号に与える影響を確認するため，人 の音声信号を所望信号とした風雑音の低減実験を行っ た。基本的な条件は4.5.1項に示したものと同様であ る。ただし，低減対象となる風雑音を含む音声信号は， 風雑音（風向Dw，v = 2.6 m/s）のみの信号に対して， 人の音声信号を伝達関数を畳み込んで加算することで 擬似的に生成した。このとき，音声信号と風雑音信号 のSN比が−30 dBとなるように混合させた。音声信 号については，25歳男性が日本語の母音である「あい うえお」を発話して収録したものを用いた。各母音は ピッチについて変化がないように注意し，また，各母 音を区切って発声した。 図–17は処理前の音源のパワースペクトログラムで ある。濃度は相対パワーレベルを表す。直流成分周辺 の低周波に強いパワーが存在し，また，各母音の周波数 成分を確認できる。図–18はPSS処理後のパワースペ クトログラムである。濃度のばらつきから，各フレー ムで雑音の推定誤差が大きいことが分かり，各フレー ムでスペクトルの形状が崩れていることが確認された。 図–19は実測伝達関数を用いたSS-BFによる処理後の パワースペクトログラムである。図–17と比較すると， はっきりとした改善は確認されなかったが，PSSに見 られた，各フレームでのスペクトルの形状の崩れは確 認されなかった。処理信号のカートシス比は，SS-BF は1.06，PSSは3.47であり，非線形ひずみの発生量 はSS-BFのほうがより少ないことが確認された。
図–18 PSS 処理後のスペクトログラム 図–19 SS-BF（実測）処理後のスペクトログラム
5. 結
論
2チャンネル近接マイクロホンアレイを用いて風雑 音を収録し，各チャンネルの信号の相関を分析した。 その結果，パワー，振幅，複素信号のどの項目について も125 Hzで0.3∼0.8の相関が得られた。複素相関を 利用した風雑音の低減では，2種類の線形ビームフォー マにより 125 Hz以下で 3∼10 dB程度のパワーの低 減が確認された。また，従来法（PSS）と提案法のカー トシス比を計算して比較した結果，提案法は従来法よ りも音質の点で優位であることが確認された。 文 献[ 1 ] S. F. Boll, “Suppression of acoustic noise in speech using spectral subtraction,” IEEE Trans. Acoust. Speech Signal Process.,27, 113–120 (1979).
[ 2 ] M. Yoshida and T. Oku, “Wind noise reduction device,” U. S. Patent, US 8,428,275 B2 (2013). [ 3 ] Y. Chung, “Rejection of flow noise using a
co-herence function method,” J. Acoust. Soc. Am., 62, 388–395 (1977).
[ 4 ] K. Kumatani, B. Raj, R. Singh and J. McDonough, “Microphone array post-filter based on spatially-correlated noise,” Proc. Interspeech 2012 (2012).
[ 5 ] R. Raspet, J. Webster and K. Dillion, “Frame-work for wind noise studies,” J. Acoust. Soc. Am., 119, 834–843 (2006).
[ 6 ] F. Shields, “Low-frequency wind noise correlation in microphone arrays,” J. Acoust. Soc. Am., 117, 3489–3496 (2005).
[ 7 ] D. Herman, “Wind noise rejection apparatus,” U. S. Patent, US 8,391,529 B2 (2013).
[ 8 ] Y. Uemura, Y. Takahashi, H. Saruwatari, K. Shikano and K. Kondo, “Automatic optimization scheme of spectral subtraction based on musical noise assessment via higher-order statistics,” Proc. Int. Workshop Acoustic Echo and Noise Control (2008). [ 9 ] M. Goto, H. Hashiguchi, T. Nishimura and R.
Oka, “RWC music database: Popular, classical, and jazz music databases,” Proc. 3rd Int. Conf. Music Information Retrieval (ISMIR 2002 ), pp. 287–288 (2002). 坂田 直人 2014年工学院大卒，2016 年同大学院修 士了，同年クラリオン（株）入社，現在に 至る。レトロゲームやその音楽に興味を持 つ。趣味はゲーム音楽コンサートの鑑賞。 日本音響学会会員。 村上 哲郎 2016年工学院大卒，同年同大学院入学， 現在在学中。ディジタル信号処理に関する 研究に従事。趣味は音楽制作。日本音響学 会会員。 中島 弘史 1994年工学院大卒，1996 年同大学院修 士了，2006 年博士（工学）。1996 年日本音 響エンジニアリング ，2006 年（株）ホン ダ・リサーチ・インスティチュート・ジャ パンを経て，2011 年より工学院大情報学 部准教授，現在に至る。恩師，東山三樹夫 先生の後継として 2 代目数理音響学研究 室を開設。日本音響学会会員。 中臺 一博 1993年東大卒，1995 年同大学院修士了。 NTT，NTT コムウェア，JST ERATO 北 野プロジェクトを経て，2003 年より（株） ホンダ・リサーチ・インスティチュート・ ジャパン，現在に至る。2003 年博士（工 学）。2006 年東工大客員准教授，2011 年早 大客員教授，2016 年東工大特定教授。ロ ボット聴覚，実時間情報統合，音環境理解 の研究に従事。人工知能学会理事，日本ロボット学会，情報 処理学会，ヒューマンインタフェース学会，IEEE 各会員。
