「AR（拡張現実）で何が変わるのか?」とか「食える数学」とか

最近読んだ2冊の本がとても面白かったので感想などをポツリと。

「AR（拡張現実）で何が変わるのか?」



「東のエデン」とのコラボ作品などで有名なAR3兄弟によるARの解説本。
AR3兄弟自らが作った作品紹介をおりまぜつつ、「ARとは何か」、「ARで何が変わるのか」、「ARの未来は」といったことに関して言及していく。

すごく興味深かったのが「都市ガスはARなんじゃないか？」というお話。
ぱっと聞いただけでは何を言ってるんだかよくわからない。
いわく、本来無味無臭のガスに匂いを付加した都市ガスというのは、ガスの「性質」みたいなものを拡張している、という観点から考えると、AR（拡張現実）なんじゃないか、と。
他にも「風鈴」とか、変わったところだと「ディズニーランド」なんてものまでそういったARに含まれるんじゃないかと。

それまで持っていた「カメラ画像に合った画像をリアルタイムで表示させるための技術」といった狭い考えを覆し、ARの本質に触れられた（ような気がした）一冊でした。


「食える数学」



数学者とエンジニア、二つの職業を経験した著者が数学の実用性を身近な例を使って紹介していく、というもの。
紹介されているのはSuica、暗号、スパムフィルタリングなど、普段何気なく使っている技術。
これらの技術にどのように数学が応用されているかについて分かりやすく解説されている。
数式もほとんど出てこないし、文章も柔らかく、しかし表現が完結で短めなので、サクサク読めます。
「数学に興味はあるけどニガテ…」といった人でも難なく読み進められると思います。

だからといって理数系が得意な工学畑の人が読んでツマラナイかといったらそんなことは全然無い。
というかむしろ、工学畑の人が見るべきなのでは？とか思ったり。

著者さんは「武器としての数学」といったことを述べられています。
欧米に追い越せ、追い付けだった昔は目標、目的が明確に定まっていたが今はそうでは無い。
自分で目標を設定し、そこにたどり着くすべを自ら模索しなければならない、と。
そういった時に必要になってくるのが数学。

この「武器としての数学」を工業大学を出た自分がうまく使いこなせているのか、と自問してみると、かなりアヤシイ…。
あくまで「手段」と割りきってる節があり、証明をよく読まないこともしばしば。
こうなると本質を理解していないために、文字通り「付け焼刃」になってしまう…。

本の中では「数学は体育」とも述べられている。
食らいつかなければ身につかないと。
真に「自らの武器」とすべく、表面を追うだけではなく、しっかりと身につけなければなーと反省させられる部分も。

数学に詳しい人も、そうでない人にもオススメな万人向けの一冊。

Bezier曲線

OpenGLで曲面を描画する方法について調べている。
その中でベジエ曲線とかB-splite曲線とかいろいろ出てきた。

文章や図を眺めてるだけだけじゃイマイチ理解が深まらないので、Procceing上でベジエ曲線描画するプログラムを自力実装してみた（注：Processingにはベジエ曲線を描画するための関数、bezier()があるので純粋に描画したいだけならそれ使った方が早いです）。

int w = 500;
int h = 500;

void setup(){
  smooth();
  noStroke();
  fill(200,0,0,200);
  size(w,h);
  float[][] pArray ={{50.,200.},
                {200.,100.},
                {300.,400.},
                {450.,200.}};
  fill(0,200,0,200);
  Bezier(pArray);
  

  fill(200,0,0,200);
  for(int i = 0;i < pArray.length;i++){
    ellipse(pArray[i][0],pArray[i][1],10,10);
  }
}

void draw(){
}

int combination(int n,int k){
  int c = 0;
  
  c = factorial(n) / (factorial(k) * factorial(n - k));
  
 return c; 
}

int factorial(int n){

  if(n == 0){
    return 1;
  }
  
  int f = 0;
  
  if(n == 1){
    return n;
  }
  else{
    f = n * factorial(n - 1);
  }
  return f;
}

float Berstein(int i,int N,float t){
  float b = 0;
  b = combination(N,i) * pow((1 - t),(N - i)) * pow(t,i);
  return b;
}

void Bezier(float[][] pArray){
  final int numT = 1000;
  final int x = 0;
  final int y = 1;
  
  for(int ti = 0;ti <= numT;ti++){
    float t = ti / float(numT);
    float px = 0;
    float py = 0;
    for(int i = 0;i < pArray.length;i++){
      float b = Berstein(i,pArray.length - 1,t);
      px += b * pArray[i][x];
      py += b * pArray[i][y];
    }
    ellipse(px,py,10,10);
  }
}

出力結果はこんな感じ。

OpenGLのキーボード入力

OpenGL触ってて、色々と引っかかったのでメモ。

・OpenGLで矢印とかファンクションキーから入力を受け付けたい場合は、glutSpecialFunc()に特殊キーのコールバック関数を設定する必要がある(C++の場合。Pythonの場合は特殊キーがどうとか考えなくてOK←追記：PyQtでOpenGLウィジェットを用いる場合の話です。OpenGL単体で使う場合は必要)。

※こんな感じ

void specialKeyFunc(unsigned char key,int x,int y){
 if(key == GLUT_KEY_LEFT){ 
  //process
 }
 else if(key == GLUT_KEY_RIGHT) {
  //process
 }
}

glutSpecialFunc(specialKeyFunc);

・文字を描画する際にはglutBitmapCharacter()を利用する（ココとかココを参考に）。

続OpenCV（光学迷彩的な何か

引き続きOpenCV。
光学迷彩的な何かを作ってみた。

やってることは至極単純で、事前に撮影した室内の画像とwebカムからキャプチャした画像の差分を計算しているだけ。


↑こんなん。左上が処理後、右上がキャプチャ画像、左下がもとの背景画像。
デコ部分の処理が甘く浮き上がってしまっている…。


それにしても、1000円かそこらのwebカムと無料のソフトウェアライブラリでここまで遊べるとは…。
引き続きいろいろ試してみたい所存でございます。

OpenCVが楽しすぎる件について

C++の勉強がてらOpenCVを触っているんだけど、これがかなり楽しい。

とりあえず、Webカムからリアルタイムに画像取得→画像から顔認識して顔部分の色を変化させる、くらいまでは簡単に出来た。




※左がカメラキャプチャした元画像、右が顔認識した部分の色を変化させたもの


一応下に使用したC++のソースを載せておく。

#include <cv.h>
#include <highgui.h>
#include <ctype.h>

int main (int argc, char **argv)
{
  CvCapture *capture = 0;
  IplImage *frame = 0;
  IplImage *out = 0;
  double w = 320, h = 240;
  int c;

  if (argc == 1 || (argc == 2 && strlen (argv[1]) == 1 && isdigit (argv[1][0])))
    capture = cvCreateCameraCapture (argc == 2 ? argv[1][0] - '0' : 0);

  cvSetCaptureProperty (capture, CV_CAP_PROP_FRAME_WIDTH, w);
  cvSetCaptureProperty (capture, CV_CAP_PROP_FRAME_HEIGHT, h);

  CvHaarClassifierCascade* cascade = cvLoadHaarClassifierCascade( "C:/OpenCV2.0/data/haarcascades/haarcascade_frontalface_default.xml",cvSize(1,1));
  CvMemStorage* storage = cvCreateMemStorage(0);

  cvNamedWindow ("Capture", CV_WINDOW_AUTOSIZE);
  cvNamedWindow("Face detection", CV_WINDOW_AUTOSIZE);
  while (1) {
    frame = cvQueryFrame (capture);
    cvShowImage ("Capture", frame);

 out = frame;

    CvSeq* objects = cvHaarDetectObjects(
      frame,
      cascade,
      storage,
      1.1,
      2,
      0,
      cvSize(30,30));

    for(int i = 0;i < (objects ? objects->total : 0);i++){
     CvRect* r = (CvRect*) cvGetSeqElem(objects,i);
        cvSetImageROI(out,cvRect(r->x,r->y,r->width,r->height));
        cvAddS(out,cvScalar(200),out);
        cvResetImageROI(out);
    }

 cvShowImage("Face detection",out);

    c = cvWaitKey (2);
    if (c == '\x1b')
      break;
  }

  cvReleaseCapture (&capture);
  cvDestroyWindow ("Capture");
  cvDestroyWindow ("Face detection");

  return 0;
}

※おまけ


お約束。

EclipseとCDTとPythonとOpenCVとわた(ry

今日も引き続きEclipseの環境構築。

OpenCVなぞをインストールしていた。

手順に関してはココを参考にしておけば間違いない。

とにかく問題になったのはOpenCVの最新バージョン2.2をインストールしてしまったこと。
上記サイトでは2.0をインストールしていて、そのとおりにやってみると確かに動く。

これが2.2をインストールしてしまうと、CDTでもPythonでも動かなかった（少なくとも自分の使っているマシンでは）…。

ソフトウェアモジュールはできるだけ安定して動くバージョンをインストールしよう、というお話。

EclipseとCDTとOpenGLと私

Eclipse+CDT+OpenGLメモ

・glutの設定はココに書いてある通り行うべし
・OpenGLプログラムを作成する際、描画関数内で頂点情報を指定したあと、swapBufferes()を忘れずに記述すること
※2010/12/12/追記:swapBuffers()は無くてもOK

ProcessingでAndroidアプリ開発

最近Processingにハマっている。

ProcessingはMITメディアラボで開発されたプログラミング言語。
メディアアートやインフォデザインなんかによく使われているらしい。

画像やインタラクションに強いため「動くモノ」、「触れるモノ」を簡単に作れるので触っててすごく楽しい。


そんなProcessingでAndroidアプリを作る方法を発見した。

…が、上記サイトの通りやってもうまく行かなかったので解決方法についてメモ。

とりあえずサイトの通りの手順だと

1.AndroidSDKをインスコ
2.Android互換バージョンのProcessingをココからゲッツ。解凍する。
3.Processingを立ち上げ、コーディング
4.コーディングし終わったら動作テスト
5.メニューバーの中にあるAndroid>Android Mode(もしくはCtrl+D)を選択。ファイルダイアログが出てくるのでAndroidSDKが保存してあるフォルダを選択する。

…ここで問題が発生した。
「Android core.zipがダウンロード出来ない」というエラーメッセージが出てきてうまくいかない。

ちょっと調べたところ、ココからファイルをダウンロードする部分がうまく動いていない様子。
上記ファイルを直接ダウンロード→Sketchフォルダに保存するとうまくいった。


…が、いざ実機にインストール、と思ったら自分が使っている端末のOSのバージョンがAndroid1.6で対応していなかった orz（エミュレータではキチンと動作しました）


まぁそれはいずれNexus Sを買ったときにでも(いつになることやら…)。


冒頭にも書きましたが、Processingはビジュアルやインタラクションに強く図形の描画、マウス、キー入力の取得が容易にできるプログラミング言語。
ということで、ゲームやデザイン系アプリのちょっとしたプロトタイプ開発に向いているんじゃないかな、と思います。

HT 03-Aの加速度センサ

HT 03-Aでプログラミングしていたときのこと。
EclipseのAndroid LogCat中に下記のような記述を発見。

10-23 14:55:42.172: DEBUG/SensorManager(4725): found sensor: AK8976A 3-axis Accelerometer, handle=0

"AK8976A"
「これはもしや、加速度センサの型番では？」
と思い、ググってみたら正解だったようで、旭化成の3軸加速度、電子コンパスICがヒットした。
データシートによるとサンプリング間隔は10~15msらしい。
つまり、SENSOR_DELAY_FASTESTで計測しても上記のサンプリング間隔より早く計測することは出来ない、ということ。

加速度センサ使ってアプリ作る人はこの点、ちょっと気を付けておいた方がいいかもです。

OEDEC2007の資料がかなり良さげ

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

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

PythonWebめぐり

ウェーブレット変換の勉強中。
Web巡回して資料を探している最中にPythonでオシロ作ってるプロジェクトを発見した（一番上のヤツ）。

これは凄い…。
PyQtとかいうGUIモジュール使って作ってるらしい。

リッチなGUIをPythonで作りたいと思ってたからちょっと勉強してみたい。
これを参考にして今まで作ってきたアプリをゴージャス（爆）に改良することを目標にやってみるか。
ちょっと調べてみた感じここら辺が勉強用に丁度よさそう。

#肝心のウェーブレット変換も同じページに載ってたのでいじくってみよう。

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.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()

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

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/

一からデジタルフィルタ設計、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結果を表示するまでの流れでプログラムが終了してしまうので、これをリアルタイムで出来るようにしたい。