OEDEC2007の資料がかなり良さげ

音声処理について調べてる途中で、OEDEC2007のセッション資料を発見したのでメモ。
これが中々良さげ（以前Blogに掲載したAMDFの分かりやすい説明なんかも載っていた）。

音声処理以外にも映像やらデータ圧縮技術なんかについても掲載されているので非常に興味深い。
時間が出来たらじっくり読んどこう。

PyAudioで早送り再生とかディレイとか

PyAudioでwavファイルを早送り再生したりディレイで再生したりしてみた。

早送りやディレイは音楽の波形データを時間軸上で伸縮させることにより実現できる。
縮ませると早く、引き伸ばすとゆっくり再生できる。

すでに音楽データを数値の配列として扱うことが出来ているわけだから、データ間を補間したり逆にデータを間引くことによりこれらを実現することができる。

試しにコーディングしてみた。

def delay(inp,rate):
    outp = []

    for i in range(len(inp)):
        for j in range(rate):
            outp.append(inp[i])

    return array(outp)

def fast(inp,rate):
    outp = []

    for i in range(len(inp) / rate):
        outp.append(inp[i * rate])

    return array(outp)

引数rateは何倍に遅くしたい（or早くしたい）かを決定するためのもの。

そのまま使うと整数倍でしか再生スピードを調整できないけど、二つを組み合わせると自在に変更可能（1.1倍再生ならfast(data,11)とdelay(data,10)の組み合わせ、とか）。

簡単なコードだけど試してみたら結構面白かった。
是非お試しあれ。

※追記：上の二つの関数を再生スピード変更関数として一元化、PyAudioでwav再生するプログラムを作ってみた。

from scipy import *
import pyaudio
import wave
import sys

def changePlaySpeed(inp,rate):
    outp = []
    for i in range(int(len(inp) / rate)):
        outp.append(inp[int(i * float(rate))])
    return array(outp)

chunk = 1024

wf = wave.open("hoge.wav", 'rb')

p = pyaudio.PyAudio()

# open stream
stream = p.open(format =
                p.get_format_from_width(wf.getsampwidth()),
                channels = wf.getnchannels(),
                rate = wf.getframerate(),
                output = True)

# read data
data = wf.readframes(chunk)

# play stream
while data != '':
    stream.write(data)
    data = wf.readframes(chunk)
    data = frombuffer(data,dtype = "int16")
    if data != '':
        data = changePlaySpeed(data,1.8)
    data = int16(data).tostring()


stream.close()
p.terminate()

changePlaySpeed()の２番目の引数で再生スピードを変更できます(2.0だったら2倍速、0.5だったら半分のスピード)。

ゆっくり再生すると変な音が混ざった感じに聞こえる。
あまりイクナイ。

再生している音源が悪いのかしら？
データの補完方法に問題があるような気もするなぁ。

とりあえず試してみました、くらいの感覚で。

グラフの対数表示とか研究室のイベントとか

今日は所属研究室の初回ゼミ&歓迎会でした。

ということで酔ってます。

今日の更新はいつにも増して内容が薄いです。

周波数特性グラフをデシベル表示にしてみました、というだけ。


こんな感じ。





前の記事に掲載した内容とほぼ変わらす。


変更点は各周波数に対するゲインを求める式を以下のように変えたのみ。

>>> 20 * log10(f[i] / fmax)

最近リアルタイム音声処理の本とか読み始めたのでそちらについても書きたいなー。

一からデジタルフィルタ設計、Pythonで実装してみる（その２）

デジタルフィルタを一から設計してPythonで実装しよう！（タイトルまんま）
…という企画の続き。

残る作業は
・数式内の定数を決める
・差分方程式を元にコーディング

たったこれだけの簡単なお仕事。


まず、数式内の定数を定める。

これは一番初めに決めた「信号から1000Hz以下の周波数成分をブっこ抜く」、という部分とデジタルフィルタの伝達関数から求めていく。

まずデジタルフィルタの伝達関数Hd(z)から周波数特性を得るため、z=e^(j2πfT)とする。
そしてその関数を用いて周波数特性を求めグラフにしたものがこちら。





遮断周波数(振幅が最大値に比べて√2/2になる周波数)が1000Hz付近になるようなRとCになっている。
R=1000、C=0.2uで決定。


これを差分方程式に代入してIIRフィルタとしてコーディングしてみる。

#coding:utf-8

from pylab import *
from scipy import *

def LPF(array):
    output = []
    
    R = 1000.
    C = 0.0000002
    Ts = 1 / 50000.

    for i in range(len(array)):
        if i == 0:
            output.append(1/(R*C)*Ts*array[0])
        else:
            output.append(1/(R*C)*Ts*array[i] + exp(-1/(R*C)*Ts)*output[i-1])

    return output

Fs = 50000.
Ts = 1 / Fs
t = arange(0,1,Ts)
f = 500

s1 = sin(2 * pi * f * t)
s2 = sin(2 * pi * 2 * f * t)
s3 = sin(2 * pi * 5 * f * t)

plot(t[0:200],LPF(s1)[0:200])
plot(t[0:200],LPF(s2)[0:200])
plot(t[0:200],LPF(s3)[0:200])
xlabel("time(s)")
legend(("500Hz","1000Hz","2500Hz"))

show()

結果がコチラ。





周波数特性と比較してみるとその通りにフィルタリングされている。

実際に周波数特性を求めたり、フィルタを作って分かったこと。
それは、アナログのローパス→インパルス不変法というデジタルフィルタの作り方だと、良い特性を持つLPFが作れないということ。
遮断周波数を過ぎてもすぐに減衰しない。

ということで実用的なデジタルフィルタを構築するには他の方式も勉強しておく必要があるようだ。

Pythonでフィルタリングツール

昨日の続き。
Python&matplotlibで移動平均フィルタリングアプリを作ってみた。
…とは言っても、かなりやっつけ仕事。

一応形にはなったので、画像とソースを張っておく。



左側がフィルタリング前の波形、右側がフィルタリング後の波形。
左下のスライダーがノイズ成分の強度調整、右下のスライダーが移動平均フィルタの強度調整。


以下が、ソースコード。

#coding:utf-8

from pylab import *
from matplotlib.widgets import Slider, Button

#FIR filter
def FIR(array,FIRCoefficient):
    output = []

    for i in range(len(array)):
        temp = 0
        for j in range(len(FIRCoefficient)):
            temp += array[i - j] * FIRCoefficient[j]
        output.append(temp)
        
    return output

#Time axis
t = arange(0.0,1.0,0.0001)
f = 10
s = []

#Noise cofficient
nc = 1

#Noisy sine wave
for i in range(len(t)):
    s.append(sin(2 * pi * f * t[i]) + nc * randn())
    
ax = subplot(121)
l, = plot(t,s)
axis([0,1,-10,10])
subplots_adjust(left=0.25, bottom=0.25)
xlabel("time(s)")
ax.set_position([0.1,0.5,0.35,0.3],"original")


#Filtered wave
FIRCoefficient = []
FIRlength = 10
for i in range(FIRlength):
    FIRCoefficient.append(1/ float(FIRlength))
s2 = FIR(s,FIRCoefficient)
ax2 = subplot(122)
m, = plot(t[:len(s2)],s2)
axis([0,1,-10,10])
ax2.set_position([0.55,0.5,0.35,0.3],"original")


#Noise controll
axcolor = 'lightgoldenrodyellow'
axnoise = axes([0.1,0.3,0.32,0.03], axisbg = axcolor)
snoise = Slider(axnoise,"Noise",0.0,10.0,valinit = 1)

def update(val):
    global s
    nc = snoise.val
    s = []
    for i in range(len(t)):
        s.append(sin(2 * pi * f * t[i]) + nc * randn())
    s = array(s)
    l.set_ydata(s)
    FIRCoefficient = []
    for i in range(FIRlength):
        FIRCoefficient.append(1/ float(FIRlength))
    s2 = FIR(s,FIRCoefficient)
    m.set_ydata(s2)
    draw()
snoise.on_changed(update)


#Filter length controll
axfilter = axes([0.55,0.3,0.32,0.03], axisbg = axcolor)
sfilter = Slider(axfilter,"Filter",5,30,valinit = 10)

def filtering(val):
    global FIRlength
    FIRlength = int(sfilter.val)
    FIRCoefficient = []
    #Moving average filter
    for i in range(FIRlength):
        FIRCoefficient.append(1/ float(FIRlength))
    s2 = FIR(s,FIRCoefficient)
    m.set_ydata(s2)
    draw()
sfilter.on_changed(filtering)

show()

昨日作ったアプリと早いところ統合したいなぁ。
でも、その前にどっちのアプリも個別にブラッシュアップしとかないと。

一からデジタルフィルタ設計、Pythonで実装してみる（その１）

これまで、デジタルフィルタの構成イメージを掴むため実際にPythonで実装してみた。

けど、本格的にフィルタを構築する場合、しっかりとした設計が必要になる。

ということで、今日はデジタルフィルタの設計について。



まず、大まかな設計の流れについて。
信号からある特定の周波数成分だけ抽出（あるいは除去）したいというニーズがあるとして、デジタルフィルタを設計するときは以下の６ステップを踏む必要がある。


(1)抽出したい成分の周波数を特定する
　目的とする信号とノイズを区別するため。

(2)目的に合わせてフィルタの種類（ローパス、ハイパス、バンドパスetc）を選択する
　高周波のノイズ除去ならローパス、FMラジオの電波受信みたいに特定の狭い周波数帯域だけ欲しい場合はバンドパス、など。

(3)選択したフィルタをアナログフィルタで実現した場合の伝達関数を求める
　この段階で得られた数式はs領域のもの。

(4)(3)で求めた伝達関数をデジタルフィルタの伝達関数に変換
　s-Z変換を用いて離散時間領域のシステムに変換する。

(5)(4)で得た伝達関数から差分方程式を導出
　伝達関数を式変形した後、z逆変換して導出。

(6)差分方程式を実際にコーディング


こんな感じ。
今までにやってきたフィルタの実装は(6)にあたる（設計がいかに大変な作業なのかがよく分かる…）。



練習がてら、上の６ステップに基づいて「信号から1000Hz以下の周波数だけ取り出すフィルタ」というものを１から設計してみる。


(1)周波数の特定
　上の題目から、取り出したい信号成分は1000Hz以下のもの。それ以上はノイズであるという前提条件を立てる。

(2)フィルタの選択
　「1000Hz以下の周波数だけ」とのことなので、ローパスフィルタを選択する。

(3)伝達関数導出
　回路理論の教科書を引っ張り出してローパスフィルタのアナログ回路構成と回路方程式、そして伝達関数を調べてみる。





図下部の数式は、回路方程式を立てる→ラプラス変換→H(s)導出、という流れになっている。
これで伝達関数H(s)が手に入った。

(4)s-Z変換を行う。これにはいくつか方法があるんだけど今回は「インパルス不変法」を使ってみる。
インパルス不変法はその名の通り、アナログフィルタのインパルス応答の標本値を利用する変換方式。

詳細の説明はコチラにゆずるとして、変換の過程はこんな数式。




めでたく、デジタルフィルタの伝達関数Hd(z)が手に入った。

※ここではオーソドックスな方法である、という理由からインパルス不変法を使った。けど本当は、前に決定した部分を振り返りながら、慎重にs-Z変換手法を選択しなければいけない。各変換手法にはそれぞれ得手、不得手があるからだ。
今回はローパスフィルタを試しに実装しているけど、実はインパルス不変法ではハイパスフィルタやバンドパスフィルタは実現できない。
ということで、今回はサクっと設計手法を決めているけど、実際の設計では１ステップずつしっかりと進めていく必要がある。当たり前の話ではあるけれど、非常に重要かつ忘れがち…。


(5)差分方程式導出
差分方程式を求めるために、デジタルフィルタの伝達関数を式変形の後、z逆変換する。
…というと難しそうに聞こえるけど簡単で、ちょこちょこっと式変形を行うだけ。





最後に出てきた数式の中のvi(k)がFIRフィルタで実装する部分、vo(k-1)がIIRフィルタで実装する部分となる。

ここまで来たらあと少し。
いくつかの定数の値を決め、コーディングするのみ。

長くなってしまったので続きはまた今度。

※続きはコチラ

PythonでIIRフィルタ

Pythonで信号処理シリーズ。
IIRフィルタ編。

ホントはscipyのsignalモジュールにFIRフィルタとかIIRフィルタを実現するためのクラスだかメソッドだかがあるので、それを使えばいいと言う話なのだけど。
キチンと仕組みを理解してフィルタを使う、プログラミングの勉強の一環として、個人的な趣味、などの理由からフィルタをPythonで実装してみようという趣旨。


前置きはこのあたりにして、本題。


前の記事にも書いたけどIIRフィルタは出力のフィードバックをそのシステム内部に含むデジタルフィルタ。

出力のフィードバックを含むというのは、プログラムの観点から言うと「フィルタ自身が過去の出力を保存しており、適宜それを使用する仕組みを持っている」ということを意味している。

これに対して、FIRフィルタは「過去の入力を保存しておき、それを適宜使用する」といったものだった。
このFIRフィルタにフィードバック機構を備え付けたものがIIRフィルタだ。
イメージとしては

FIRフィルタ+フィードバック機構=IIRフィルタ

といったような感じ（と、自分は捉えてる）。


それを踏まえてココを参考にIIRフィルタを以下のような形で設計してみた。





中央の加算器を挟んで左半分がFIR部分、右側がIIR部分となっている。
今回はIIRフィルタの影響のみを調べたいので、FIRフィルタの係数はインプットをそのまま出力する形となっている。


んでソースは下のようなものに。

#coding:utf-8

from scipy import *
from matplotlib.pylab import *

def IIR(array):
    output = []
    
    #係数格納行列
    FIRCoefficient = [1,0,0]
    IIRCoefficient = [0.7,0.2,0.1]
    
    for i in range(0,len(array)):
        temp = 0

    for i in range(len(array)):
        temp = 0
        for j in range(len(FIRCoefficient)):
            temp += array[i - j] * FIRCoefficient[j]
            if len(output) > j:
                temp += output[i - j - 2] * IIRCoefficient[j]
            
        output.append(temp)

    return output

Fs = 100.
t = arange(0,1,1/Fs)
f = 5

inp = []

#入力信号生成(正弦波+ノイズ)
for i in range(len(t)):
    inp.append(sin(2 * pi * f * t[i]) + rand())

outp = IIR(inp)

plot(t,inp,"b")
plot(t,outp,"r")
xlabel("time(s)")
show()

出力結果がコチラ。





…発振してしまっている。

どうやらIIRフィルタの係数設定はかなり難しいみたい。
FIRフィルタのときは直感的に係数を決めて、目で見て分かるくらいの効果があるフィルタを作ることができた。
でもIIRの場合、そううまくはいかないようで。

係数の設定や構成など、設計段階で考慮しなければいけないことが多いみたいだ。
ということで次回はデジタルフィルタ設計の勉強をすることにするか。

PyAudioでピッチ検出

これまで、信号処理のアルゴリズムと音声入力を取り扱う方法について色々勉強してきた。
その中からPyAudioとAMDFアルゴリズムを組み合わせて、マイクから音声を拾って音程を探り当てる、ということをやってみた。

まず、マイクに向かって「ド」の音階を出すつもりで発声してみる。





拡大してみると





それらしい波形になっている。


これにAMDFをかけたものがコチラ。





おおよそ0.0078(s)のあたりにピークが来ている。
これが声の基本周期だと考えると周波数は、1 / 0.0078 ≒ 128 (Hz)となる。

このサイトによると低い「ド」の音は131(Hz)とのこと。
そんなに大きくハズレてはいない模様。

少しホッとした。


現時点では、録音してAMDF結果を表示するまでの流れでプログラムが終了してしまうので、これをリアルタイムで出来るようにしたい。

信号処理の勉強その２

引き続き信号処理の勉強。

デジタルフィルタについて調べていた。
以下ざっくりとしたまとめ。

(1)デジタルフィルタとは
・フィルタとは信号から特定の周波数成分を抜き出す、あるいは除去するためのモノ
・離散値を処理するフィルタをデジタルフィルタという
・デジタルフィルタには線形フィルタと非線形フィルタがある（現時点ではこういうくくりがある、というぐらいで留めておく）
・デジタルフィルタは3つの基本的な構成要素から成る
　ー遅延素子：1クロック分、信号成分を遅延させて出力する素子
　ー加算器：入力された信号成分を足し合わせて出力する素子
　ー乗算器：入力された信号成分を定数倍して出力する素子

(2)FIR(Finite Inpulse Responce)フィルタ
フィードバックループが無いデジタルフィルタ。
インパルス応答が有限個のパルスで表されるからFIRフィルタ（らしい）。

(3)IIR(Infinite Inpulse Responce)フィルタ
フィードバックループを持つデジタルフィルタ。
インパルス応答が無限に続くインパルス列になるのでIIRフィルタ（らしい）。

全部ココを参考にして2、3行にまとめてみた（ザックリしすぎな感は否めない）。


ということで毎度お馴染みPythonでFIRフィルタを実装してみる。

ざっくり下のような感じで設計。
四角が遅延素子、三角が乗算器（内側の数字が係数）、丸が加算器。




んで、コーディングした結果がこちら。

#coding:utf-8

from scipy import *
from matplotlib.pylab import *

#FIRフィルタ
def FIR(array):
    output = []
    
    for i in range(2,len(array)):
        output.append(array[i] * 0.7 + array[i-1] * 0.2 + array[i-2] * 0.1)

    return output

Fs = 100.
t = arange(0,1,1/Fs)
f = 5

inp = []

#入力信号生成(正弦波+ノイズ)
for i in range(len(t)):
    inp.append(sin(2 * pi * f * t[i]) + rand())

outp = FIR(inp)

plot(t,inp,"b")
plot(t[2:len(t)],outp,"r")
xlabel("time(s)")
show()

出力結果がコレ。




青が入力信号で赤が入力をフィルタリングした信号。

少し分かりづらいけれど、よく見るとフィルタリングを行う前と後で波形のギザギザが少なくなっていることが分かる。これは高周波のノイズがフィルタによって除去されたことを意味する。
ごくごく簡単なローパスフィルタ（低周波数の信号のみ通すフィルタ）が完成した。

ということで遅延素子を用いるとローパスフィルタを作ることが出来る、ということがこれで分かった。

AMDF

AMDF(Average Magnitude Deference function)を用いたピッチ検出を試してみましたよ、というお話。

AMDFは自己相関、相互相関関数のようにピッチ検出に用いられるアルゴリズムのひとつ。

（※2010/10/12追記:OEDEC2007のプレゼン資料にAMDFの分かりやすい説明が載ってます）

ココとかココ（こちらはIEEEの認証が必要）を参考にAMDFをPythonで実装してみた。

def AMDF(array,frameLength):
    amdf = []
    length = len(array)
    for i in range(frameLength):
        temp = 0
        for j in range(length):
            temp += abs(array[j] - array[j - i])
        amdf.append(temp / length)
    return amdf

この関数が正しく動作するか確認するため、下のような波形を用意する。



周波数100Hzの方形波にノイズを加えたもの。


この波形にAMDFを適用して計算した結果がこちら。



0.01sに最小値のピークが出ている。
相互相関関数と違って、AMDFでは最小値となる点が周期性を表していることが分かる。

信号処理の勉強

SciPyも無事導入出来たということで、信号処理を勉強し始めた。
ということで調べたことをメモ。

・自己相関関数
時間領域において波形の周期性を調べる際に用いられる数学的手続き。
その波形自信の畳み込み積分を計算することによって求められる。
音声のピッチ検出などに応用されている。

・相互相関関数
2つの波形の類似性を調べるための数学的手続き。
自己相関関数と同様ピッチ検出などに応用されている。
2つの波形の畳み込み積分を計算することによって求められる。

どちらもピッチ検出に用いられているけれど、相互相関を使った方が自己相関を使うよりも計算負荷が低くいのでより良いとのこと。

毎度お馴染みPythonだったら、2つともnumpy.correlate()で計算可能。

ただし、自己相関関数を求める場合は下記のようにmodeに"same"を指定する必要あり。

>>> correlate(signal,signal,"same")

ただ使うだけじゃ面白くないので相互相関関数を自分で実装してみた。
CCFは(Cross Correlation Function：相互相関関数)の略。

def CCF(array1,array2):

    if len(array1) < len(array2):
        array1,array2 = array2,array1
        
    ccf = []
    
    for i in range(len(array1) - len(array2)):
        temp = 0
        for j in range(len(array2)):
            temp += array2[j] * array1[j - i]
        ccf.append(temp / len(array2))
        
    return ccf

この関数を用いて10Hzの正弦波（サンプリング周波数1000Hz）と、その正弦波の先頭から1周期分の波形の相互相関を求めてみる。



青が元の正弦波、赤が正弦波の先頭1周期、緑がそれら二つの波形の相互相関。

相互相関の値が0,100,200…と100の倍数で最大値をとっている。

これはつまり、この正弦波が100ms(10Hz)の周期の波形であること示している。
…ということでキチンと計算できているようだ。