ディジタル信号処理のためのPython入門

ディジタル信号処理のための Python 入門
第 2 回. 畳み込み積分とインパルス応答
赤部晃一
AHC Lab
YYYY-MM-DD
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
1 / 11
畳み込み積分
定義域が (−∞, ∞) の 2 つの関数 f と g の畳み込み積分は
∫ ∞
(f ∗ g)(t) =
f (τ )g(t − τ )dτ
−∞
離散の場合
(f ∗ g)[n] =
∑
f [m]g[n − m]
m
畳み込み積分のフーリエ変換は、2 つの関数のフーリエ変換の積
[f ↔ F , g ↔ G] −→ [(f ∗ g) ↔ F G]
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
2 / 11
インパルス関数（デルタ関数）
次を満たす関数 δ
(f ∗ δ)(t) = f (t)
インパルス関数の性質（連続: ディラックのデルタ関数）
∫ ∞
δ(x)dx = 1
−∞
δ(x) = 0 (x ̸= 0)
離散信号では、次の「クロネッカーのデルタ関数」が使われる
{
1 (x = 0)
δ[x] =
0 (x ̸= 0)
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
3 / 11
インパルス応答 h
あるシステムにインパルス信号 δ[t] を入力した際の出力信号
あるシステムのインパルス応答 h[t] と信号 x[t] の畳み込み積分
(h ∗ x)[t] は、あるシステムに信号 x[t] を入力した際の出力信号
y[t] である。
y[t] = (h ∗ x)[t]
YYYY-MM-DD
±
?
h
x
?
y
Koichi Akabe
(AHC-Lab)
4 / 11
無限インパルス応答（IIR）
ある時点における出力信号が、過去の出力信号に依存するような
システムのインパルス応答。
y[t] =
N
∑
bm x[t − m] +
m=0
N
∑
an y[t − n]
n=1
x
b0
z-1
a1
y
b1
z-1
a2
b2
z-1
a3
ここで z −1 は 1 サンプル分の遅延を表す。
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
5 / 11
有限インパルス応答（FIR）
ある時点における出力信号が、有限時間の入力信号のみに依存す
るようなシステムのインパルス応答。
y[t] =
N
∑
bm x[t − m]
m=0
x
b0
z
y
-1
b1
z-1
b2
出力が有限個の入力で「必ず」安定するため、広く利用される。
（IIR システムは発散する場合がある。）
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
6 / 11
音声ファイルの読み込み
wave モジュールを使って読み込み、numpy を使って数値データ化
する
import wave
import numpy
def main () :
wav = wave . open ( " ./ sample . wav " , " rb " )
rate = wav . getframerate ()
nframes = wav . getnframes ()
raw = wav . readframes ( nframes )
x = [ xx / 2**15 for xx in numpy . frombuffer ( raw , "
int16 " ) ]
# :::
wav . close ()
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
7 / 11
SciPy における FFT の実装
第 1 回で実装した FFT は、サンプル数が 2 のべき乗でなければ実
行できない。
SciPy にはこの問題に対処した FFT 関数が存在する。（ただし 2
のべき乗の際が最も速い）
今回の問題ではこの関数を利用しても良い。
http://docs.scipy.org/doc/scipy/reference/generated/scipy.ﬀtpack.ﬀt.html
from scipy . fftpack import fft
def main () :
# :::
H = fft ( h )
# :::
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
8 / 11
問題: 有限インパルス応答に基づくフィルター
次のインパルス応答を示すシステムを考える。

sin[2π(n−8)c/f ]

(n ̸= 8, 0 ≤ n ≤ 16)

π(n−8)
h[n] =
2c/f
(n = 8)


0
(17 ≤ n)
ただし、c = 4000[Hz]、f = 11025[Hz]（音声信号の標本化周波数）
▶
▶
▶
畳み込み積分を実装せよ。
h[n] の区間 [0,16] をプロットせよ。
sample.wav の先頭の 11025 サンプル x[n] を切り出して FFT
を実行し、その結果を dB 表示でプロットせよ。
単位変換: [dB] = 20 log10 X
▶
▶
YYYY-MM-DD
x[n] と h[n] の畳み込み積分に対して FFT を行い、その結果
を dB 表示でプロットせよ。
h[n] を [0,11025) の範囲で FFT した結果を dB 表示でプロッ
トし、このシステムの概要を説明せよ。
Koichi Akabe
(AHC-Lab)
9 / 11
実行結果
0.8
60
0.6
40
20
h
X [dB]
0.4
0.2
0.0
−0.2
0
−40
2
60
4
6
8
n
10
12
14
−60
0
16
4000
6000
8000
Frequency [Hz]
10000
2000
4000
6000
8000
Frequency [Hz]
10000
−20
0
H [dB]
Y [dB]
20
−20
−40
−60
−40
−60
YYYY-MM-DD
2000
0
40
−80
0
0
−20
−80
2000
4000
6000
8000
Frequency [Hz]
10000
Koichi Akabe
0
(AHC-Lab)
10 / 11
付録: 畳み込み積分のコード
f と g の合計の長さを持つ、結果を格納するリスト conv を作成
し、各要素について積を計算。
conv = [0.0] * ( len ( f ) + len ( g ) )
for n in range ( len ( conv ) ) :
for m in range ( n + 1) :
if m >= len ( f ) :
continue
if n - m >= len ( g ) :
continue
conv [ n ] += f [ m ] * g [ n - m ]
if 文を毎回行うのは大変。これを range に内包すると、
conv = [0.0] * ( len ( f ) + len ( g ) )
for n in range ( len ( conv ) ) :
for m in range ( max (n - len ( g ) +1 ,0) , min (( n +1) , len ( f ) ) ) :
conv [ n ] += f [ m ] * g [ n - m ]
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
11 / 11

Download Report