2010年9月26日日曜日

Pythonで少数を分数に変換する関数作った

前の記事で書いたwavファイルの速度調整用にスクリプト書いたので載せておく。

少数を分数に変換するための関数。


def Euclidean(m,n):
if m < n:
m,n = n,m
if n == 0:
return m
elif m % n == 0:
return n
else:return Euclidean(n,m%n)

def deciToFrac(deci):
i = 0
while deci % 1 != 0:
deci = deci * 10
i += 1

numer,denomi = int(deci),int(10 ** i)

gcd = Euclidean(numer,denomi)

return numer / gcd,denomi / gcd



下の関数に少数を渡すと通分された分数に変換し、分子と分母を返す。
上の関数はユークリッドの互除法。
これを使って分子、分母を通分している。


使用結果はこんな感じ。



>>> deciToFrac(2.6)
(13, 5)
>>> deciToFrac(1.2)
(6, 5)

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番目の引数で再生スピードを変更できます(2.0だったら2倍速、0.5だったら半分のスピード)。

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

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

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

2010年9月24日金曜日

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

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

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

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

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


こんな感じ。





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


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



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



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

2010年9月22日水曜日

一からデジタルフィルタ設計、Pythonで実装してみる(その2)

デジタルフィルタを一から設計して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()


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

2010年9月21日火曜日

PyAudioとmatplotlibで録音アプリ作った

今日、Twitterのタイムライン上に「JavaScriptでペイントアプリ作る。一日で」というのを行っている人がいたので、自分も何か作りたくなった。

ということでmatplotlibやらPyAudioの勉強がてら簡単な録音アプリをPythonで作ってみた。

速攻でコーディングしたせいでソースの見栄えが非常に悪いので(あとまだまだ不完全なので…)、GUIのスクリーンショットだけ載せておく。






左上のラジオボタンでサンプリングレートを、グラフ下部のスライダーをいじくると録音時間を設定できる。
Saveボタンで録音した音声をwav形式で保存、Recordボタンを押すと録音開始、という非常にシンプル(芸がないとも言う)なもの。


時間があったら、再生ボタンとかエフェクターとかフィルタとか色々実装してみたい。

こういうアプリを作ってると、GUIを自由にデザインして実装できるようになりたいなぁ、としみじみ思う。
matplotlibも悪くないけどやっぱりウィジェットの配置とかデザインとか一から全部手掛けたい。

Ajax(javascript)とかHTML5とかweb周りの技術触っときたいなー。

※参考にしたサイト
[1]wavファイルの保存処理
http://www.s12600.net/psy/python/02-1.html
[2]matplotlibのGUI周り
http://matplotlib.sourceforge.net/examples/widgets/

2010年9月19日日曜日

一からデジタルフィルタ設計、Pythonで実装してみる(その1)

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

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

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



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


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

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

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

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

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

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


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



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


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

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

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





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

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

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




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

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


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





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

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

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

※続きはコチラ