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

ディジタル信号処理のための Python 入門
第 5 回. フィルター
赤部晃一
AHC Lab
YYYY-MM-DD
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
1 / 16
フィルター
信号から特定の範囲の周波数を除去
HPF
1.5
1.5
1.0
1.0
0.5
0.0
0
0.5
1000
4000
0.0
0
5000
1.5
1.5
1.0
1.0
0.5
YYYY-MM-DD
1000
2000
3000
Frequency [Hz]
4000
5000
4000
5000
BEF
2.0
H
H
2000
3000
Frequency [Hz]
BPF
2.0
0.0
0
LPF
2.0
H
H
2.0
0.5
1000
2000
3000
Frequency [Hz]
4000
5000
Koichi Akabe
0.0
0
(AHC-Lab)
1000
2000
3000
Frequency [Hz]
2 / 16
理想的なフィルターのインパルス応答
遮断周波数 ωc の低域通過フィルターの場合
n = 0 のとき:
ωc
h[n] =
π
n ̸= 0 のとき:
h[n] =
=
1
2π
∫
ωc
ejΩn dΩ
−ωc
[
]ω
1 ejΩn c
2π jn −ωc
1 ejωc n − e−jωc n
2π
jn
sin(ωc n)
=
πn
インパルス応答が無限に続き、かつ非因果的（n < 0 で値を持つ）
なため、実現できない。
=
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
3 / 16
インパルス応答の打ち切り（FIR フィルター）
インパルス応答を 2N + 1 の長さで打ち切った場合

sin[ωc (n−N )]

(0 ≤ n ≤ 2N, n ̸= N )
 π(n−N
)
ω
c
h[n] =
(n = N )
π


0
(n > 2N )
因果性を確保するため、中心を N とした。
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
4 / 16
インパルス応答の打ち切り結果
遮断周波数 3000Hz, 標本化周波数 10000Hz (ωc = 0.6π)
N=8
0
20
H [dB]
H [dB]
20
40
60
40
60
80
100
80
0
N=50
20
0
2000
4000
6000
Frequency [Hz]
8000
120
0
10000
2000
N=200
10
4000
6000
Frequency [Hz]
8000
10000
8000
10000
N=4000
0
0
10
20
30
H [dB]
H [dB]
20
40
50
40
60
60
70
80
0
▶
▶
▶
YYYY-MM-DD
80
2000
4000
6000
Frequency [Hz]
8000
10000
0
2000
4000
6000
Frequency [Hz]
N が小さいほどリアルタイム性が高いが、周波数応答が理想
的なフィルターとは異なる。
N が大きくなると理想的なフィルターに近づくが、リアルタ
イム性が低くなる。
遮断周波数付近にリップル（波）がある。
Koichi Akabe
(AHC-Lab)
5 / 16
バタワースフィルター
バタワース・ローパスフィルター
1
{
(
H(s) = ∏
N −1
s − ωc exp 12 +
k=0
(2k+1)π
2N
)}
振幅特性
1
(
|H(jω)| = √
1+
ω
ωc
)2N
ここで、N はフィルターの次数
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
6 / 16
バタワース・ローパスフィルターの特性
wc = 1, N = 1, 2, 3, 4
0
H [dB]
−20
−40
−60
−80
−100 -1
10
100
Frequency [rad/s]
101
通過域が平坦（リップルが無い）。傾斜は約 −6N [dB/oct]
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
7 / 16
バタワースフィルターのディジタル化
1+sT /2
双一次変換による近似 z = esT ≈ 1−sT /2 を用いる。近似による誤
差を考慮して、遮断周波数の補正が行われることもある。
SciPy によるバタワースフィルターの実装
http://docs.scipy.org/doc/scipy/reference/generated/
scipy.signal.butter.html
from scipy . signal import butter , lfilter
# :::
# x : signal
# n : order
# fc : critical frequency
# fs : sample frequency
# low pass / high pass
b , a = butter (n , fc / fs , " lowpass " )
# band pass / band stop
b , a = butter (n , [ fc1 / fs , fc2 / fs ] , " bandpass " )
y = lfilter (b , a , x )
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
8 / 16
第一種チェビシェフ多項式
−1 ≤ x ≤ 1 の場合、
Tn (x) = cos(n cos−1 x)
漸化式に展開すると
T0 (x) = 1
T1 (x) = 0
Tn (x) = 2xTn−1 (x) − Tn−2 (x)
展開された Tn の定義域を拡張し、x > 1 に対しても考える。（解
析接続）
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
9 / 16
第一種チェビシェフ多項式の出力
10
8
C
6
4
2
0
−1.0
YYYY-MM-DD
−0.5
0.0
Koichi Akabe
0.5
x
(AHC-Lab)
1.0
1.5
2.0
10 / 16
第一種チェビシェフフィルター
第一種チェビシェフ多項式を利用したフィルター
第一種チェビシェフ・ローパスフィルター
振幅特性:
1
|H(jω)| = √
1 + ϵ2 Tn2
(
ω
ωc
)
（伝達関数は割愛）
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
11 / 16
第一種チェビシェフ・ローパスフィルターの特性
ϵ = 1, wc = 1, N = 2, 3, 4, 5
0
G [dB]
−20
−40
−60
−80
−100
10
-1
10
0
10
1
Frequency [rad/s]
バタワースフィルターに比べ勾配が急。通過域にリップルがある。
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
12 / 16
SciPy における IIR フィルターの実装
バタワースフィルター、第一種・第二種チェビシェフフィルター、
楕円フィルターなどがある。
詳細は SciPy のドキュメントを参照
http://docs.scipy.org/doc/scipy/reference/signal.html#matlabstyle-iir-filter-design
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
13 / 16
音声ファイルへの書き出し
音声信号 y をチャンネル数: 1, 量子化ビット数: 16（2 バイト）,
標本化周波数: 44100Hz で出力する例
wav = wave . open ( " ./ upsample . wav " , " wb " )
wav . setnchannels (1)
wav . setsampwidth (2)
wav . setframerate (44100)
wavdata = [ min ( max ( int ( yy * 2**15) , -2**15) , 2**15 -1)
for yy in y ]
wav . writeframes (
struct . pack ( " h " * len ( wavdata ) , * wavdata ) )
wav . close ()
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
14 / 16
問題: 標本化周波数変換
11025Hz で標本化された音声信号 sample.wav の標本化周波数を
以下の手順によって 44100Hz に変換せよ。
1. sample.wav を読み込み、配列化せよ。
2. 配列の各要素を、インデックスが 44100/11025 倍となる位置
に移動し、空となった要素は 0 で埋めよ。
3. 遮断周波数 11025/2Hz（ナイキスト周波数）のローパスフィ
ルターを設計し、上で作成した信号に対して適用せよ。
4. ローパスフィルター通過後の信号を標本化周波数 44100Hz の
音声ファイルとして upsample.wav に保存せよ。
5. 両者の音声を比較せよ。
6. この方法によって標本化周波数を変換できる理由を考察せよ。
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
15 / 16
付録: 配列インデックスの引き伸ばしと 0 埋め
配列 x の要素を、インデックスが scale 倍となる位置に移動し、
空となった要素は 0 で埋める。
y = [0.0] * int ( len ( x ) * scale )
ref_pref = -1
for i in range ( len ( y ) ) :
ref = int ( i * scale )
if ref != ref_pref :
ref_pref = ref
y [ i ] = x [ ref ]
YYYY-MM-DD
Koichi Akabe
(AHC-Lab)
16 / 16

Download Report