卒
業 研 究 報 告
題 目
DSP を用いた音声処理
指 導 教 員矢野 政顕 教授
報 告 者 学籍番号 : 1060251 氏名 : 野本 善介平成
18 年 2 月 21 日
高知工科大学 電子・光システム工学科
目次
第
1 章 はじめに
... 1第
2 章 ディジタル信号処理
... 2 2.1 ディジタル信号処理の基礎 ... 2 2.1.1 ディジタル信号処理とは... 2 2.1.2 ディジタル信号処理システムの基本構成... 4 2.2 DSPとは... 6 2.3 スペクトル解析 ... 6 2.4 DFT... 6 2.4.1 DFTの定義... 6 2.4.2 DFTの性質... 7 2.5 FFT ... 8 2.6 FFTアルゴリズム... 8 2.6.1 2 を基底とする周波数間引きFFTアルゴリズム... 8 2.7 窓掛け... 10 2.7.1 窓関数 ... 14第
3 章 音声処理システム
... 16 3.1 音声処理システムの構成 ... 16 3.2 FFTによるスペクトル解析 ... 16第
4 章 スペクトル解析のシミュレーション
... 19第
5 章 まとめ
... 26謝辞
... 27参考文献
... 28付録
... 29第
1章 はじめに
現在、我々の身の回りには様々な認証システムがある。例えば、手の指紋を読み 取って認証する指紋認証、眼球の黒目に現れる皺のパターンを識別して認証を行う 虹彩認証、手のひらや指の静脈のパターンを識別して認証を行う静脈認証が挙げら れる。他にも、撮影した顔の画像から、傾きや位置を検出して補正し、特徴点(眼 の中心、唇の端など)の位置や点同士の距離などを計測し認証を行う顔認証などが ある。これらには、画像を取り込み処理を行う画像処理の技術が使われている。し かし、高い認証率を得るには高速な処理装置や大量の記憶装置が必要となり、シス テム構築が高価になる。 比較的安価に実現可能な音声による認証システムとしてDSP(Digital Signal Processor)による音声認証システムが注目されている。それにはまず音声というア ナログ信号をディジタル信号に変換する必要がある。ディジタル信号に変換する理 由として、近年、コンピュータ技術や集積回路の飛躍的な発展により様々な分野で ディジタル化が進んでいる。また、ディジタル信号処理よってアナログ回路では設 計困難で複雑なシステムを、高精度でかつ安価に実現できる。それに伴い、ディジ タル信号の演算をいかに高速にかつ効率よく処理するかが重要となってくる。そこ で、乗算などの演算を高速かつ効率よく処理できるDSPを用いることで、よりよい 音声認証システムを構築できると考えられる[1][2]。 本研究の目的は、DSPを用いたディジタル信号処理技術の基本を学び、音声認証 として実現できる音声処理システムを検討することである。 本論文の第2章ではディジタル信号処理について説明する。第3章では、音声処理 システムの構成について説明する。第4章では、スペクトル解析のシミュレーショ ン結果について説明する。第5章で全体をまとめる。第
2章 ディジタル信号処理
2.1 ディジタル信号処理の基礎
2.1.1 ディジタル信号処理とは
我々の身の回りの物理的な現象の多くは、時間と共に変化し連続的な値をとるア ナログ信号の形でとらえられる。ディジタル信号処理とは、このアナログ信号をデ ィジタル的に処理することである。 近年、ディジタル信号処理は幅広い分野で利用されている。利用されている例と して、ディジタル・オーディオ機器や携帯電話、通信、音声、画像、制御、計測な どの分野がある。ディジタル信号処理が様々な分野で応用される理由は、色々と利 点があるからである。しかし、欠点もある。その利点と欠点を以下に示す[1]。 <利点> ・ プログラムとして表現できる処理は、原理的にどんなものでも実現でき、ア ナログ信号処理の得意でない非線形の処理、適応処理が可能。 ・ 情報圧縮・誤り訂正・暗号化などが容易であり、データの蓄積・伝送に都合 がよく、付加価値が期待できる。 ・ ビット長を長くすれば高精度化が実現できる。 ・ 温度・湿度や経年による変化がないので、安定した品質が実現できる。 ・ 小型化・低電力化・低コスト化できる。 ・ 仕様の変更などに柔軟に対処できる。 <欠点> ・ リアルタイム処理を行う場合、一つの入力サンプル値に対して次のサンプル 値が入力される前に処理を完了しなければならない。そのため、扱える信号 の帯域が制限される。 ・ 十分なビット長が確保できない場合、ディジタル化や演算で生じる誤差を避 けられない。 ・ 簡単な処理の場合でも、A/D・D/Aコンバータなどの周辺機器が必要になるディジタル信号処理を行うためには、まずアナログ信号をディジタル信号に変換 する必要がある。ディジタル信号とは、アナログ信号をある一定の間隔(サンプリ ング間隔)Tごとにサンプリング(標本化)した離散的信号のことである。図2-1の(a) にアナログ信号、(b)にサンプリングされた信号を示す[1][3]。
(a)
(b)
図 2-1 アナログ信号とサンプリングされた信号2.1.2
ディジタル信号処理システムの基本構成
図2-2にディジタル信号処理システムの基本的な構成を示す。
入力されたアナログ信号をサンプリングする前に低域通過フィルタ(LPF:Low Pass Filter)によって帯域を制限をする。次に、A/D C (A/Dコンバータ)によってサ ンプリングと量子化を行いアナログ信号をディジタルシステムが扱える形に変換 する。
アナログ
信号
ディジタル
システム
A/D C
LPF
D/A C
LPF
アナログ
信号
図 2-2 ディジタル信号処理システムの基本構成 サンプリングを行う際に、サンプリング間隔が小さければ元のアナログ信号を正 確に復元できるが、単位時間当たりのデータ量が増えてしまう欠点がある。サンプ リング間隔を決める基準に式(2.1)に示すサンプリング定理を用いる。このとき、元 の信号の周波数をf0、サンプリング周波数をfsとする。 fs ≧ 2f0 (2.1) このサンプリング定理を満たさずにサンプリングを行うと、元の信号を正確に復 元できない。サンプリング定理を満たしたサンプリングと満たしていないサンプリ ングの例を図2-2に示す。サンプリング周波数を10kHzとし、図2-2(a)に3kHzの信 号、(b)に7kHzの信号、(c)、(d)に(a)、(b)それぞれをサンプリングした信号を示し、 (e)、(f)には(c)、(d)それぞれを復元した信号を示す。(a) (b) (c) (d) (e) (f) 図 2-2 サンプリング例 図2-2 より、サンプリング定理を満たしている 3kHz の信号は正確に復元されて いるが、サンプリング定理を満たしていない7kHz の信号は正確に復元されていな いことが分かる。3kHz、7kHz それぞれに復元された信号は同じ波形になっている。 このように、サンプリング定理を満たしていない場合に復元された信号が元のアナ ログ信号とは違う別のアナログ信号として復元される現象をエイリアシング (aliasing)と呼ぶ[1]。
2.2
DSP とは
DSP(Digital Signal Processor)とはディジタル信号の演算に特化した、ディジタ ル信号処理専用のマイクロプロセッサのことをいう。ディジタル信号処理では主に 式(2.2)に示す積和の計算を行う。
α
1χ
1+
α
2χ
2+
L
+
α
nχ
n (2.2) DSPは式(2.2)のような積和の計算を高速に実行できるように設計されている。そ して、DSPはマイクロプロセッサの一種なのでどのような処理を行うかはプログラ ムで決まる。プログラムを書き換えれば、異なった処理を行うことが可能となるた め、DSPを用いることの最大のメリットとなる。 例えば、ある方式で作られたシステムを新しい方式に変更しなければならなくな った時、そのシステムが DSP を用いていればプログラムを変更するだけで簡単に 新しい方式に対応できる[1]。2.3 スペクトル解析
スペクトル解析とは、ある信号に含まれる周波数成分の分布を求めることである。 本研究では、FFT を用いたスペクトル解析を行う。2.4 DFT
2.4.1 DFT の定義
DFT (Discrete Fourier Transform : 離散的フーリエ変換)のことをいう。 ある連続時間系の時系列信号を g(t)、g(t)のフーリエ変換を G(ω)とすると、g(t) とG(ω)はそれぞれ式(2.3)、(2.4)のようにあらわされる。
∫
∞ ∞ − = ω ω ω π G j t d t g ( )exp( ) 2 1 ) ( (2.3)G(
ω
) =∫
−∞∞g(t)exp(− jω
t)dt (2.4) ディジタル信号処理では、信号g[t]をある一定のサンプリング間隔Tでサンプリ ングした時系列信号を扱う。サンプリングした信号をg[n]とし、g[t]を周期 NT の 周期信号と仮定すると式(2.3)、(2.4)は式(2.5)、(2.6)のようにあらわされる。∑
− = − = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = 1 0 1 , , 1 , 0 , 2 exp ] [ 1 ] [ N k N n N nk j k G N n g π L (2.5) 1 , , 1 , 0 , 2 sin 2 cos ] [ 2 exp ] [ ] [ 1 0 1 0 − = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − =∑
∑
− = − = N k N nk j N nk n g N nk j n g k G N n N n L π π π (2.6) このような変換を離散的フーリエ変換と呼ぶ[1]。 式(2.6)の複素指数関数を簡潔にあらわすために、式(2.7)を使う。)
/
2
sin(
)
/
2
cos(
)
/
2
exp(
N
j
N
N
j
W
Nπ
π
π
−
=
−
=
(2.7) WNを回転因子と呼ぶ。WNを使うと式(2.6)は式(2.8)のように変換できる[5]。
∑
− = − = = 1 0 1 , , 1 , 0 , ] [ ] [ N n nk N k N W n g k G L (2.8)2.4.2 DFT の性質
DFT の性質として、信号g[n],(0≦n≦N -1)が実数の場合、G[k]とG[N -k]は式 (2.9)のような複素共役になる。G
[
k
]
=
G
[
N
−
k
]
*
(2.9) *は複素共役をあらわす。G[k]を実部と虚部にわけて式(2.9)の関係をあらわした ものを式(2.10)に示す。{
}
{
}
{
[
]
}
Im
{
[
]
Im
]
[
Re
]
[
Re
k
N
G
k
G
k
N
G
k
G
−
−
=
}
−
=
(2.10) このとき、G[k]の虚部について Im{G[0]}と Im{G[N /2]}は必ず 0 になる。 よって、実信号のDFT を求める場合、式(2.8)の計算をk = 0, 1,…, N /2 について 計算すればよい。また、実信号g[n]が偶関数の場合、もとめた DFT の虚部がすべ て0、奇関数の場合は実部がすべて 0 になる[1][5]。2.5 FFT
FFT (Fast Fourier Transform : 高速フーリエ変換)とはDFTを効率よく実行す るアルゴリズムのことである。DFTを式(2.6)の定義通りに計算すると、計算量は N2に比例し、Nが少し大きくなると計算量が膨大になる。計算量が増えると計算時 間も増え、重大な問題となる。この問題を解決する為に提案されたアルゴリズムが FFTである。 FFT を用いてスペクトル解析を行う場合、二つ気を付けなければならない点があ る。一つは、DFT の対象となる信号は時間について連続なアナログ信号をサンプ リングしていることである。これについては、2.1.2 節で述べているサンプリング 定理を満たしていれば解消される。もう一つは、DFT のデータが有限長であるこ とである。この問題点については2.7 節にて述べる[1]。
2.6 FFT アルゴリズム
FFTアルゴリズムには色々ある。よく使われるのは 2 を基底とする(radix-2)アル ゴリズムと呼ばれる。データ数が 2M(M :整数)になるときに適用できるFFTアル ゴリズムである。2.6.1 節にて 2 を基底とするアルゴリズムを使った周波数間引き (decimation-in-frequency)FFTアルゴリズムについて述べる。2.6.1 2 を基底とする周波数間引き FFT アルゴリズム
2 を基底とする周波数間引きFFTアルゴリズムは、図 2-3 に示す基本演算の組み 合わせによって構成される。この基本演算はバタフライ演算(butterfly operation) と呼ぶ。N =22=4 として周波数間引きFFTアルゴリズムの例を図 2-4 に示す。図 2-3 バタフライ演算
a
b
A
B
k NW
-1
k NW
b
a
B
b
a
A
)
(
−
=
+
=
図2-4 周波数間引き FFT アルゴリズム 第1 段目]
0
[
G
]
2
[
G
]
1
[
G
]
3
[
G
]
0
[
g
]
1
[
g
]
2
[
g
]
3
[
g
-1 0 4W
-1 0 4W
-1 1 4W
-1 0 4W
第2 段目 図2-4 を例に、処理手順を説明する。まず、1 段目ではデータを 2 つのグループ に分ける。次に、N /2 だけ離れたところのデータ同士をバタフライ演算する。 2 段目では、1 段目で分けられたグループをさらに 2 つに分ける。そして、N /4 だけ離れたデータ同士をバタフライ演算して、得られたG[k]が DFT の結果となる。 この例では、2 段の演算によって結果が得られたがデータ数が 2M(M :整数)の ときは、M段の演算によってDFTの結果が得られる。 図2-4 より注意する点として、入力データg[n]は上から昇順に並んでいるが、求めら れた G[k]は昇順に並んでいない。図 2-4 のような G[k]の並び方をビット逆順(bit reversal)と呼ぶ。よって、この周波数間引き FFT アルゴリズムを使用する場合、 ビット逆順に並んでいるデータ G[k]を入力データ g[n]と同様に昇順に並び換える 必要がある。N =8 の場合の通常の並び方とビット逆順の並び方の対応を表 2-1 に 示す[5]。表2-1 ビット逆順対応表(N =8) 通常 ビット逆順 10 進数 2 進数 2 進数 10 進数 0 000 000 0 1 001 100 4 2 010 010 2 3 011 110 6 4 100 001 1 5 101 101 5 6 110 011 3 7 111 111 7
2.7 窓掛け
窓掛けとは、無限長のデータの一部分を取り出し有限長のデータにする操作のこ とである。通常、フーリエ変換は無限長の信号に対して定義される。しかし、DFT は有限の長さのデータから計算される。よって、DFT の計算を行う場合、窓掛け を行い無限に続くデータを有限のデータにする必要がある。窓掛けの過程を図 2-5 に示す[1]。 図 2-5 において、本来の信号に窓掛けを行い(b)の窓の内部の信号のみを取り出 す。そして、(c)のように窓の外側は、窓の内部の信号が周期的に続いているものと 仮定する。(c)の信号波形でわかるように窓の内側と外側の境目で信号が不連続とな っている。これは、信号の周期が窓の幅の整数倍になっていないためである。例え ば、3Hz の信号を幅が 1 秒の窓で窓掛けしたとする。信号の周期Tは式(2.11)のよ うになる。 3 1 1 = = f T (2.11)(a) 元の信号 時間 時間 (b) 窓掛けした信号 時間 (c) DFT で求める信号 図 2-5 窓掛けと DFT で求める信号の様子 したがって、Tを整数倍していけば窓の幅と一致し、窓の内部と外部の境目が不 連続にならない。信号の周期が窓の幅の整数倍になっている場合となっていない場 合、DFT で振幅スペクトルを求めた時にどういった影響がでるのか図 2-6、図 2-7 に示す。窓の幅を1 秒とし、図 2-6(a)に 3Hz の信号に窓掛けをした場合、(b)に 3.2Hz の信号に窓掛けをした場合を示す。図2-7(a)と(b)に DFT で求めた図 2-6(a)、(b)そ れぞれの振幅スペクトルを示す[1]。
データ長 (1 秒) (a) 3.0Hz の信号 データ長 (1 秒) (b) 3.2Hz の信号 図 2-6 窓掛けをした信号 図2-6 より、3.0Hzのスペクトルは、周波数成分が 3Hz のところのみ存在して いるので元の信号を正しく表しているといえる。しかし、3.2Hz のスペクトルでは、 周波数成分が複数存在している。よって、元の信号を正しくあらわしているとはい えないが、3.2Hz に一番近い 3Hz の周波数成分が最も大きくあらわれている。 一般的にスペクトル解析を行う場合、信号の周期の整数倍と窓の幅が一致すると は限らない。したがって、元の信号のスペクトルと DFT で得られるスペクトルは 多少異なったものになる。スペクトルの違いを小さくするために、窓関数(window function)を用いる[1]。
0 10 20 (dB) -40 -20 0 (Hz) (a) 3.0Hz のスペクトル 0 10 20 (dB) -40 -20 0 (Hz) (b) 3.2Hz のスペクトル 図 2-7 DFT で求めた振幅スペクトル結果
2.7.1 窓関数
窓関数とは、取り出したデータに重み付けを行う関数のことである。窓掛けされ た信号の DFT は、信号のスペクトルと窓関数のスペクトルを周波数領域において 畳み込みしたものである。スペクトル解析では主に4 つの関数が使われる。それぞ れの窓関数の式を以下に示す。また、それぞれの窓関数のスペクトルを図2-8 に示 す。この時、Lは窓の幅とし、取り出したデータの数に対応する[1]。 1. 方形(rectangular)窓 (2.12) ⎩ ⎨ ⎧ ≤ ≤ − = それ以外 , 0 1 0 , 1 ] [n n L wR 2. ハニング(Hanning)窓(
)
⎩
⎨
⎧
−
≤
≤
−
=
れ以外
そ
,
0
1
0
,
/
2
cos
5
.
0
5
.
0
]
[
n
n
L
n
L
w
Nπ
(2.13) 3. ハミング(Hamming)窓(
)
⎩
⎨
⎧
−
≤
≤
−
=
それ以外
,
0
1
0
,
/
2
cos
46
.
0
54
.
0
]
[
n
n
L
n
L
w
Mπ
(2.14) 4. ブラックマン(Blackman)窓(
)
(
)
⎩
⎨
⎧
−
+
≤
≤
−
=
それ以外
,
0
1
0
,
/
4
cos
08
.
0
/
2
cos
5
.
0
42
.
0
]
[
n
n
L
n
L
n
L
w
Bπ
π
(2.15) 図2-8 より、4 つのスペクトル全て、中央部に大きい山が存在する。そして、大 きい山の左右に小さい山が存在している。大きい山はメインローブ(main lobe)、小 さい山はサイドローブ(side lobe)と呼ばれている。メインローブの幅ができるだけ 狭く、サイドローブの大きさができるだけ小さいものが、理想的な窓関数とされて いる。しかし、図2-8 からわかるように、メインローブの幅が狭いとサイドローブ が大きい、サイドローブが小さいとメインローブの幅が広いという具合に、一方が よければ他方が悪いという関係にある。 FFT と用いてスペクトル解析を行う場合、解析の目的や信号の性質(a) 方形窓 (b) ハニング窓
(c) ハミング窓 (d) ブラックマン窓 図 2-8 窓関数のスペクトル
第
3章 音声処理システム
3.1 音声処理システムの構成
本研究では、入力された音声の振幅スペクトル解析を DSP で行う音声処理シス テムを検討した。音声処理システムのブロック図を図3-1 に示す。サンプリング
窓掛け
FFT
絶対値の計算
終了
入力信号
図 3-1 音声処理システムの構成 主な流れは、入力された信号をサンプリング定理に基づいて定めたサンプリング 周波数によってサンプリングする。次に、窓関数を用いて窓掛けを行う。窓を掛け て取り出した信号を FFT によって、周波数成分に分ける。最後に絶対値を計算し て終了する。3.2 FFT によるスペクトル解析
元となる音声信号は、周波数の異なる複数の信号を重ね合わせた擬似信号として、 音声処理のシミュレーションを行う。程だからである。また信号の周期の整数倍と窓の幅が一致させられない場合に、一 般的にハミング窓がよく使われるからである。 したがって、標本化した信号g[n]と式(2.14)に示したハミング窓WMの乗算によっ て求められた式(3.1)のデータに対してDFTを行う。
g
[
n
]
×
W
M[
n
],
n
=
0
,
1
,
L
,
N
−
1
(3.1) FFT のプログラムをできるだけ高速なものにするには、回転因子とビット逆順の 並べ替えを考慮する必要がある。回転因子とビット逆順の計算についてはあらかじ め値を計算しておき、FFT の計算を実行する時に呼び出すという形で使うようにす ると高速なプログラムが実現できる。 まず、回転因子 について述べる。 の値はk =0,1,・・・,(N / 2)-1 の範囲の値 になる。そして、 の虚部-sin(2πk /N)と実部 cos(2πk /N)は式(3.2)のような 関係になる。 k N W WNk k N W(
)
⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − N N k N k N k 2 /4 cos 2 2 cos 2 sin π π π π (3.2) よって、cos(2πk /N)をk =0,1,・・・,(3N / 4)-1 の範囲で計算しておけば十分であ る。 次に、ビット逆順の並べ替えについて述べる。ビット逆順の並び替えの方法にも 色々な方法があるが、本研究では図3-2 に示す方法を使う[5]。このとき、N =4 と する。10
00
00 にN /21=10 を加える⎭
⎬
⎫
10
00
00、10 それぞ れ にN /22=01 を加える⎭
⎬
⎫
11
01
図 3-2 ビット逆順の求め方(N =4 の場合)入力信号 g[n]が実信号の場合、FFT アルゴリズムの 1 段目の計算を減らすこと ができる。まず、図2-3 のバタフライ演算のBについて考えてみる。このとき、a, bを複素数として、Bを実部と虚部に分けると式(3.3)のようにあらわされる。
{ }
{
}
(
)
{
}
(
)
{ }
B
{
a
b
}
(
k
N
)
{
a
b
}
(
k
N
N
k
b
a
N
k
b
a
B
/
2
cos
Im
/
2
sin
Re
Im
/
2
sin
Im
/
2
cos
Re
Re
π
π
)
π
π
×
−
+
×
−
=
×
−
−
×
−
=
(3.3) しかし、1 段目の入力はすべて実数なので、Im{a-b}=0 となる。よって、式(3.3) は式(3.4)とあらわされる。{ }
{
}
(
)
{ }
B
{
a
b
}
(
k
N
)
N
k
b
a
B
/
2
sin
Re
Im
/
2
cos
Re
Re
π
π
×
−
=
×
−
=
(3.4) したがって、計算量を半分にすることができる。また、Aの出力の虚部も0 にな るため計算量を減らすことができる。 他の段においても、計算量を減らすことのできる箇所は存在する。しかし、計算 量を減らすためのif 文などの条件文は必要となり、かえって実行時間が増加する可 能性がある。したがって、ここでは他の段についての説明は省略する[5]。 最後に、FFTを行い得られたDFTの結果G[k] の絶対値|G[k]|2を求める。信号 g[k]が実信号ならば、G[k]は式(3.5)になるので絶対値はk =0,1,・・・,N / 2 について求 めればよい。G
[
N
−
k
]
=
G
[
k
]
*
(3.5) 絶対値の計算を式(3.6)に示す。さらに、結果を dB で表示する場合の計算を式(3.7) に示す。{ }
(
Re
[
]
)
(
Im
{ }
[
]
)
,
0
,
1
,
,
/
2
]
[
k
2G
k
2G
k
2k
N
G
=
+
=
L
(3.6)10
log
10G
[
k
]
2,
k
=
0
,
1
,
L
,
N
/
2
(3.7) 以上のことを考慮し、FFT によるスペクトル解析を行うプログラムを作成する。 作成したプログラムを付録の<リスト1>に示す。第
4章 スペクトル解析のシミュレーション
3 章で作成したスペクトル解析のプログラムが正しく動作するかを、シミュレー ションを行い確認する。前述したように、入力信号は周波数の異なる複数の信号を 重ね合わせた擬似信号とする。 本研究では、サンプリング定理を満たしている擬似信号PS1、PS2 とサンプリン グ定理を満たしていない擬似信号PS3 の計 3 種類の擬似信号を作成した。サンプリ ング周波数は256Hzとし、それぞれの擬似信号の周波数を以下に示す。 擬似信号PS1 5Hz、13Hz、20Hz、32Hz、41Hz、56Hz、67Hz、72Hz、89Hz 擬似信号PS2 7Hz、15Hz、21Hz、30Hz、39Hz、46Hz、67Hz、83Hz、98Hz 擬似信号PS3 180Hz、193Hz、206Hz、214Hz、224Hz、235Hz、250Hz、423Hz、 486Hz まず、擬似信号PS1 についてのシミュレーション結果を示す。図 4-1 に擬似信号 PS1 の波形を示し、図 4-2 に窓掛けされた状態のPS1 の波形を示す。PS1 のスペク トル解析結果を図4-3 に示す。 図 4-1 擬似信号PS1 の波形次に、擬似信号PS2 についてのシミュレーション結果を示す。図 4-4 に擬似信号
PS2 の波形を示し、図 4-5 に窓掛けされた状態のPS2 の波形を示す。PS2 のスペク
トル解析結果を図4-6 に示す。
最後に、擬似信号PS3 についてのシミュレーション結果を示す。図 4-7 に擬似信
号PS3 の波形を示し、図 4-8 に窓掛けされた状態のPS3 の波形を示す。PS3 のスペ
クトル解析結果を図4-9 に示す。
図 4-3 と図 4-6 のサンプリング定理を満たしているPS1、PS2 それぞれのスペク トル解析結果を見ると、擬似信号として重ね合わせた周波数成分のところにスペク トルが存在している。したがって、正しくスペクトル解析が行えたことがわかる。 しかし、サンプリング定理を満たしていないPS3 のスペクトル解析結果の図 4-9 を見ると、重ね合わせた周波数成分とは別の周波数成分のところにスペクトルが存 在している。したがって、正しくスペクトル解析されず、エイリアシングが発生し ていることがわかる。 これらの結果から、サンプリング定理を満たしている信号については正しくスペ クトル解析を行うことが可能といえる。
第
5章
まとめ
本研究では、音声認証として実現できる DSP を用いた音声処理システムの検討 を行った。DSP を用いることにより、ディジタル信号処理の基礎について学習する ことができた。また、スペクトル解析のプログラムを作成することにより、音声処 理の基礎やC 言語など様々な知識を身につけることができた。 今後の課題として、擬似信号ではなく本来の音声信号のスペクトル解析を行う必 要がある。また、スペクトル解析を正しく行うことができたので、認証システムの 実現の可能性の検討を行う必要がある。さらに、今回作成したプログラムでは演算 にかかる時間や演算量が大きくなっているので、演算の処理方法に改良を加え演算 時間・演算量を小さくすることが必要である。謝辞
本研究を進めるにあたり、日ごろより懇切丁寧な御指導と御鞭撻くださいました 高知工科大学工学部 電子・光システム工学科 矢野 政顯 教授に深く感謝いたしま す。 また、御助言を頂くとともに日頃からお世話になりました原 央 教授、橘 昌良 助教授、山本 真行 講師、植田 和憲 講師、他各先生方に厚く御礼申し上げます。 さらに、本研究を進めるにあたり御助言と御指導頂ました山岡 大祐氏をはじめ、 石元 啓一氏、国沢 千寿氏、日頃からお世話になった矢野研究室、原研究室、橘研 究室の皆様方に心から感謝いたします。参考文献
[1] : 三上 直樹 : 「ディジタル信号処理の基礎」 CQ出版社 (2003) [2] : 高橋 進一、池原 雅章 : 「ディジタルフィルタ」 培風館 (1999) [3] : 辻井 重男、鎌田 一雄 : 「ディジタル信号処理」 昭晃堂 (1996) [4] : 谷萩 隆嗣 : 「ディジタル信号処理と基本理論」 コロナ社 (1996) [5] : 三上 直樹 : 「C言語によるディジタル信号処理入門」 CQ出版社 (2002)付録
ここでのプログラムは、擬似信号にPS1のスペクトル解析の場合のものである。 <リスト 1> スペクトル解析のプログラム #include <math.h> #include <stdio.h> #define dNUM 256#define dNUM2 (dNUM/2) #define dNUM3 ((dNUM*3)/4) const float PI = 3.1415926536; const float s0 = 1.0/32768.0; float wn[dNUM3];
short br[dNUM];
float u0[dNUM], uW[dNUM], xr[dNUM], xi[dNUM], xp[dNUM2],xp2[dNUM2]; float Ps[dNUM];
float SD[dNUM2]; int c;
inline void swap(float *a, float *b);
void fftable(float wn [], short br [], int D_NUM) {
int i, d_half, ne, jp; float arg;
arg = (float)6.283185307/D_NUM; for (i=0; i<((D_NUM*3)>>2); i++) {
wn[i] = (float)cos(arg*i); }
d_half = D_NUM>>1; br[0] = 0;
for (ne=1; ne<D_NUM; ne=ne<<1) {
for (jp=0; jp<ne; jp++) br[jp+ne] = br[jp] + d_half; d_half = d_half >>1;
} }
void ffReal (float xR[], float xI[], float wn [], short br [], int D_NUM) {
float xtmpR, xtmpI, xRD2;
int j, jnh, k, jxC, jxS, ne, d_half, d_half2; d_half = D_NUM>>1; jxC = 0; jxS = D_NUM>>2; for (j = 0; j < d_half; j++) { jnh = j + d_half; xtmpR = xR[j]; xR[j] = xtmpR + xR[jnh]; xI[j] = 0.0; xtmpR = xtmpR - xR[jnh]; xR[jnh] = xtmpR*wn[jxC++]; xI[jnh] = xtmpR*wn[jxS++]; } d_half = d_half>>1;
for (ne=2; ne<(D_NUM>>2); ne=ne<<1) {
d_half2 = d_half<<1;
for (k=0; k<D_NUM; k=k+d_half2) { jxC = 0; jxS = D_NUM>>2; for (j=k; j<(k+d_half); j++) { jnh = j + d_half; xtmpR = xR[j]; xtmpI = xI[j]; xR[j] = xtmpR + xR[jnh]; xI[j] = xtmpI + xI[jnh]; xtmpR = xtmpR - xR[jnh]; xtmpI = xtmpI - xI[jnh];
jxS = jxS + ne; } } d_half = d_half>>1; } xRD2 = xR[0] - xR[1] + xR[2] - xR[3]; for (k=0; k<D_NUM ; k=k+4) { xtmpR = xR[k]; xtmpI = xI[k]; xR[k] = xtmpR + xR[k+1] + xR[k+2] + xR[k+3]; xI[k] = xtmpI + xI[k+1] + xI[k+2] + xI[k+3]; xR[k+2] = xtmpR + xI[k+1] - xR[k+2] - xI[k+3]; xI[k+2] = xtmpI - xR[k+1] - xI[k+2] + xR[k+3]; } xR[1] = xRD2; xI[1] = 0.0; for (j=0; j<dNUM; j++) if (j<br[j]) { swap(&xR[j], &xR[br[j]]); swap(&xI[j], &xI[br[j]]); } }
inline void swap(float *a, float *b) { float tmp; tmp = *a; *a = *b; *b = tmp; } void SD_F(float xx[]) { int i,h,m=0; for( h=0; h<dNUM2; h++) {
}
for(i=0; i<dNUM2; i++) { if (xp2[i-1]<xp2[i]) { if (xp2[i]>xp2[i+1]) { SD[m]=i; m++; } else { SD[m]=0; } } else { SD[m]=0; } } } void main(void) { int i,h,n; float radius1 = 0.0; float wM, pi2N; float f[9] = {5.0, 13.0, 20.0, 32.0, 41.0, 56.0, 67.0, 72.0, 89.0}; for(i=0; i<dNUM; i++)
{
Ps[i] = 0.0; }
for(n=0; n<9; n++) {
for(i=0; i<dNUM; i++) {
} } pi2N = (float)(2.0*PI)/(float)dNUM; fftable(wn, br, dNUM); for (h=0; h<dNUM; h++) { u0[h] = s0*(float)Ps[h]; wM = (float)(0.54 - 0.46*cos(pi2N*h)); uW[h] = u0[h]*wM; xr[h] = uW[h]; }
ffReal(xr, xi, wn, br, dNUM); for (h=0; h<dNUM2; h++) { xp[h] = xr[h]*xr[h] + xi[h]*xi[h]; } for (h=0; h<dNUM2; h++) { xp[h] = (xp[h] == 0) ? (float)-100.0 : (float)(10.0*log10(xp[h])); } SD_F(xp); }