STFTとISTFTの秘密の関係
タイトルに深い意味はありません。
まあ、ここの記事で意味なんかあった試しがありませんけどね。
今回は、FFTW3ライブラリを用いたSTFTとISTFTについて、『変換後の逆変換で、元データを再現できるのか』という内容をお送りします。
ちなみに、今回実装は省きます。
需要もないですし、ソース貼り付けると重くなるので止めときます。
別に自分が書いたソースを見られるのが嫌なわけではありません。
そういうことなので、実装方面についてはあまり触れません。
面倒くさいので。
おや？ 何か見えましたか？
しかし、敢えて言いましょう。
世の中、気にした方が負けなのです。
……違うか。
STFTについて:
まあ、特にSTFTに関しては情報過多なくらいなので触れることもないでしょう。
ただ、僕の場合は元のデータの長さに合わせて、STFTする回数を動的に変える手法を使っていたので、その話でもしようと思います。
『元のデータ長に対して何回STFTが可能か』は、『いくつ紙テープを貼り付ければ期待する長さになるか』という問題に似ています。
STFTにおけるサンプリング数を紙テープの一枚あたりの長さに、サンプリングの間隔を紙テープの糊しろの長さに見立てれば、何となく分かってもらえると思います。
たとえば、サンプリング数を8、サンプリング間隔を4とすれば、それは8cmの長さの紙テープを糊しろ4cmでつなぎ合わせていくことと同じです。
さて、問題は元データの長さに合わせてSTFTの回数をどう変動させるかです。
元データ長に対してサンプリング数8、サンプリング間隔4のSTFTは何回実行できるか。
8cmの紙テープを4cmの糊しろで張り付けていき、cm以下で最も長くする場合に紙テープは何本必要か。
ここで、ならば、どうでしょう？
をそれぞれ回数と紙テープの数に見立てましょう。
紙テープは、より
です。
これは、最初の一枚を除いた結果なので、この場合、15枚でぴったり64cmになります。
STFTの回数も同じです。
すなわち、あるデータ長のデータに対して、サンプリング数
、サンプリング間隔
のSTFTを適用する場合の最大STFT回数
は、
で求められます。
は、紙テープと同様に最初の一枚分です。
これで、元データ長に対して、STFTの回数を動的に求めることができます。
ISTFTについて:
ぶっちゃけ、今回の本題はここです。
まず、前提を思い出してください。
忘れた？ 困った人ですね。
これは、FFTW3ライブラリを使っている場合を想定しているということですよ。
まったく、僕が前置きで下らないことしか言わないと思ってましたね？
概ね、その認識で合ってます。
だからこそ、ここで再確認したわけで。
いいとして。
ここからは、長さ64のデータに対して、サンプリング数8、サンプリング間隔4でSTFTしたものをISTFTすると仮定して話を進めます。
FFTW3を使っている以上、注意しなければならないことがいくつかありますね。
まず、逆変換で得られる値をSTFTしたときのサンプル数で割らなければなりません。
これは、公式マニュアルの通りです。
さて、ISTFTの逆変換のイメージは大丈夫でしょうか？
まず、STFTなので得られた周波数特性も複数あると思います。
それをそのまま戻してつなぎ合わせると、どうなるか。
ここで、前節で話した紙テープの話を思い出してください。
ISTFTした結果をそのままつなぎ合わせることは、糊しろ0で貼り付けることと同じです。
要するに、間違いです。
紙テープの糊しろ、つまりサンプリング間隔は4なので、逆変換した結果も当然4ずつ重ねて繋げていかなければなりません。
ここで、別の問題が発生します。
重ねていくということは、重複するということです。
だからと言って、逆変換で得られた結果を勝手に削ってはいけません。
この問題は、該当する部分を重複した数で割ることで解決します。
ちょっとイメージし辛いですね。
まず、紙テープを思い浮かべてください。
長さ8cmの紙テープを8分割しましょう。
これを、Cの配列に見立ててください。
このとき、最初と最後の例外を除いて、配列の添字0から3は前の結果と、4から7の4つの値は後の結果とそれぞれ重複します。
なので、基本的には2で割ればよいのです。
繰りかえしますが、逆変換の結果は全て意味を持っているので、確固たる意味がある場合を除いて、勝手に削ってはいけません。
さて、ISTFTの面倒なところは、この『重複する数』がフレキシブルな点です。
つまり、サンプリング数とサンプリング間隔の組み合わせ次第では、単純に2で割っても期待する結果を得られないということです。
ここで、『重複する数』の法則は以下となります。
※いずれも最初と最後は例外的な扱いになります。
- サンプリング数 == サンプリング間隔*2
- 重複は必ず2つ。
- 中間は必ず重複する。
- サンプリング数 < サンプリング間隔*2
- 重複は必ず2つ。
- 中間に重複しない箇所がある。
- サンプリング数 > サンプリング間隔*2
- 重複は複数。
- 中間は必ず重なる。
- 最初の方と最後の方で、重複数が段階的に変化する。
面倒くさいですよねー。
まあ、数式で求めることもできるんです。
できますが、とにかく面倒です。
この問題に対する僕のとった解決策は、カウンター配列を作ることでした。
要するに、元データ長と同じ長さの配列を作り、結果を格納する毎にカウントしていくんです。
ちょっと、想像の世界に飛び立ちましょうか。
言語はCでいきましょう。
そこ、嫌いとか言わない。
まず、ISTFTの結果格納配列resultを用意します。
次に、重複数をカウントする配列countも用意します。
これらの配列は、ともに元データの長さ分、ここでは64だけ確保します。
また、STFTの結果が15個のfftw_complex*ということはオーケーですね？
では、いきましょう。
まず、最初の逆変換で、8個の結果が返されます。これはSTFTのサンプリング数と同じです。
なので、resultの0から7番目にそれぞれ適切に結果を格納します。
合わせて、countの0から7番目もインクリメントします。
次の逆変換も同様に8個の結果を得られます。
ここで、サンプリング間隔が4であることを思い出してください。
したがって、resultの4から11番目に結果を格納することになります。
同様に、countも4から11番目をインクリメントします。
これを15回繰り返し、最後に
for(int i=0;i<64;i++) result[i] /= count[i];
みたいにすれば終わりです。
どうでしょう、イメージできました？
もちろん、最初に書いた通り、全ての結果はサンプリング数、つまりここでは8で割る必要があることも忘れてはいけません。
この方法はメモリを大量に消費しますが、一番楽な解決法だと思います。
最後に窓関数の扱いですね。
STFTの際に窓関数を用いていたのであれば、ISTFTでも同様の窓関数を使う必要があります。
STFTの際は乗算だったので、ISTFTでは結果に対して除算を適用すればよいのです。
僕の場合はWAVの音声ファイルに対して、これらの方法を使いました。
結果から言って、単純にSTFTしてISTFTしたものはオリジナルと同じでした。
diffをとろうが、cmpしようが、同じファイルという結果を得られました。
まあ、WAVに戻す処理を施す前はかなり細かいところで差異がありましたけどね。
e-30以下の領域で。
そんなわけで、『変換後の逆変換で、元データを再現できるのか』は、『ほぼできる』と結論づけていいと思います。
編集後記:
いやあ、久々に長文書いて疲れました。
まあ、今回も走り書きのメモみたいなものなので、もしかしら足りない部分があるかもしれません。
え？
『足りてる部分があるの？』……だと……？
……そうですね。
個人的には、これ以上の不足は必要ないですね。
こんな感じで。
今回は、終わりたいと思います。見つかっちゃったね！きらっ☆
また告知か……
告知ばっかりですみません。
気がつけば、まさかの2ヶ月放置でした。リアルでは色々あったのにね（微笑）
いいとして。
ご存知の方もいらっしゃるでしょうが、僕はGAEでブログを作ってます。
て言うか、作り過ぎました。
でも、流石に手が回らず、放置サイトがいくつも点在しているのはまずいです。
そこで、この度サイト群を統合するサイトを新設し、年内にも旧サイト群は閉鎖しようと思っています。
で。
統合先の新設サイトは既に作りました。ここからは、このサイトだけを拡充していきます。
代わり映えしないですが、よろしくお願いします。（愛想がつきていなければ）
"AQUA NORA"
このはてなブログも、その内まともに更新します。
それでは、またその日まで。
なんだ、ただの告知か
はい、そうです、これはただの告知です。
久しぶりなのに、もっと何か無いのかと言われそうですが、無いものは無いんです。
まあ、僕の現状に照らし合わせるなら、有るものの方が少ないですが。
さておき。
以前もご紹介しました、GAEでやってるサイトをちょっと大幅に更新しました。ここと同様、今後も更新が滞ると思いますが、暇つぶしにでもなれば幸いです。
目の毒になるかもしれませんが、自己責任でお願いします。
物置サイト
"Aqua Nora's warehouse"
それでは、本日はこれにて。
もう笑うしかないときのPython〜実行速度編〜
どーも、色々詰んでる者です。
あんまり更新してないと、僕の現在状況を知っている人にあらぬ誤解を受けて捜索願を提出されてしまう恐れがあるので、どーでもいいネタでWebを汚してやろうと思います。
有り体に言うと、ただの八つ当たりです。
まあ、バカの馬鹿話はおいといて。
今回は、実行速度に関する、ちょっとした実験結果です。
背景にあるのは、ここでも何回か登場したPythonでWAVをどーたらこーたらの話なんですが、要点をかいつまんで話すと、『Python遅くねぇ？』です。
さすがに、Cなら0.22秒でできることがPythonだと7.69秒かかるというのは、いくら仏様でも一度目でブチ切れます。
山羊、指令書食ってる場合じゃありません。
で。
当然、この無能の塊こと僕のコーディングですから、何かしら不備があるに決まっています。
そこで、何とか実行速度を上げたいという思いから、今回の実験を行ったわけです。
前置きはこの辺で。
本題に入りましょう。
2011/7/31 追記：
この記事が気になった方は、コメントも参照してください。
検証実験：
以下、testN(Nは0以上の整数)と記述されているものは、テストに用いたサンプルコードのtestN.pyを指します。
実験環境は、CPUはCore-i7、メモリ4GB、Ubuntu（64bit）です。
また、実行速度の計測はUNIX/Linuxのtimeコマンドを使いました。
appendとリスト内包表記：
test1とtest2は、リストの中身を別リストへ移す(コピーではない)プログラムです。実行速度はそれぞれ、test1が3.64秒、test2が2.29秒でした。
この結果より、可能な限りリスト内包表記を用いることが、速度向上に有益であると考えられます。
hex(ord())とbinascii.b2a_hex()：
test3とtest4は、読み込んだバイナリデータをPythonで扱える形の16進数に変換しています。Pythonでは\x00形式がとんでもなく扱い辛いので、この実験を行いました。結果ですが、test3が6.40秒でtest4が7.92秒でした。
このことから、同じ結果を得るならhex()とord()を組み合わせた方が早いことが分かります。
list型とarray型：
test5はtest2との比較用プログラムで、要素の移し替えという点では同じです。結果、test5は5.31秒かかり、test2の2.29秒に遠く及びませんでした。
この結果から、array型の利用が速度向上には繋がらないと考えられます。
arrayのappendメソッドとextendメソッド：
test5ではarrayの要素追加にappendを使いましたが、test6ではextendを使います。extendのそもそもの使い方とは違いますが、一応の比較実験です。結果は、test5の5.31秒に対して、test6は15.99秒でした。
結論は言うまでもなく、appendの方がまだマシ、ということです。
arrayのappendメソッドと＋メソッド：
test7では要素の追加に＋を使っています。気になる実行速度ですが、何と12.43秒もかかり、appendを使っていたtest5の5.31秒の足元にも及びませんでした。
この結果から、どうあがいてもappend有利のようです。
リストと文字列：
test8はtest1と同じような操作ですが文字列への移し替えという点で異なります。結果ですが、test8が3.96秒でtest1の3.51秒、さらにtest2の2.29秒と比べると、速度的には劣ることが分かりました。
コピー元と同じ長さのリストを作っておく：
test9はtest1と同様の処理を行いますが、コピー元となるリストと同じ長さのリストを作っておくという工程が入ります。移す際も、要素の入れ替えとなります。結果、test9は4.67秒でした。
ただ、こっちのほうが早くなった例があるそうで、何だかよく分からない結果になってしまいました。
テストに使ったサンプルコード：
### all
# sample.wav: 44.1kHz, mono, 16bit
# 3:56.17=10415232samples(by SoX)
def load():
fp=open("sample.wav","rb")
tmp=fp.read()
fp.close()
return tmp
raw=load()
### test1.py
txt=[]
for i in raw:
txt.append(i)
print len(txt),len(raw),txt[0],raw[0]
### test2.py
txt=[i for i in raw]
print len(txt),len(raw),txt[0],raw[0]
### test3.py
txt=[hex(ord(i)) for i in raw]
print len(txt),txt[0:30]
### test4
import binascii
txt=["0x"+binascii.b2a_hex(i) for i in raw]
print len(txt),txt[0:30]
### test5
import array
a=array.array("c")
for i in raw:
a.append(i)
print len(a),len(raw),a[0],raw[0]
### test6.py
import array
a=array.array("c")
for i in raw:
a.extend(i)
print len(a),len(raw),a[0],raw[0]
### test7.py
import array
a=array.array("c")
for i in raw:
a+=array.array("c",i)
print len(a),len(raw),a[0],raw[0]
### test8.py
txt=""
for i in raw:
txt+=i
print len(txt),len(raw),txt[0],raw[0]
### test9.py
txt=range(len(raw))
for i,j in enumerate(raw):
txt[i]=j
print len(txt),len(raw),txt[0],raw[0]
最後に：
まあ、ここまで散々やった挙句に、と言うか、散々やったために。
最近はC実装を検討しています。
だって、C早いんだもん……
まあ、世の中こんなもんです。
でも、Pythonは楽なので、今後も使いつづけてはいくと思います。
無論、趣味の領域で、ですが。
それでは、また〜。僕が生きていれば、ですが。
GaucheのFFTとFFTW3のFFTとscipy.fftpackのFFTと
それぞれの計算結果の比較テストです。
どうしてこんなことをしたかと言うと、pythonのscipyモジュールにおけるfftpackのfftが、FFTWを使っている（ラッパーである）というような記述を見かけた気がして、その確認をしてみたいと思ったからです。
確認手段ですが、FFTWであるならば別に2のべき乗個のサンプル数でなくともFFTが計算可能なはずです。まあ、その程度の安易な判断基準で調べてみようと思ったわけです。
ここ（FFTW公式ページ）に、こう書いてあります。『For example, the standard FFTW distribution works most efficiently for arrays whose size can be factored into small primes (2, 3, 5, and 7), and otherwise it uses a slower general-purpose routine.』
要するに、小さな素数に分解できるサンプル数であれば最適に計算できる。それ以外は遅いということで、どんな個数でもFFTを計算できるようになっている（と僕は思っています）ということです。
という訳で、早速実験結果を示していこうと思います。
出典1：GaucheのFFTの実装はこちらを使用しています。
出典2：FFTW3(C)の実装はこちら(リンク切れ確認しました)を参考にしています。
テストコード:
Gauche:
(define fft
(let ((pi 3.14159265358979323846)
(unshuffle
(lambda (seq)
(let loop ((rest seq)
(evens '())
(odds '()))
(if (null? rest)
(cons (reverse evens) (reverse odds))
(loop (cddr rest)
(cons (car rest) evens)
(cons (cadr rest) odds)))))))
(lambda (sequence)
(let ((len (length sequence)))
(if (= len 1)
sequence
(let ((nth-root (make-polar 1 (/ (* 2 pi) len)))
; principal len-th root of unity
(half-len (quotient len 2))
(packs (unshuffle sequence)))
(let loop ((step 0)
(root 1)
(evens (fft (car packs)))
(odds (fft (cdr packs)))
(front '())
(rear '()))
(if (= step half-len)
(append (reverse front) (reverse rear))
(loop (+ step 1)
(* root nth-root)
(cdr evens)
(cdr odds)
(cons (+ (car evens) (* root (car odds)))
front)
(cons (- (car evens) (* root (car odds)))
rear))))))))))
(define (make-test-data start len)
(let ((pi (* (atan 1) 4)))
(if (eq? len start) '()
(cons (* 0.25 (sin (/ (* start pi) 16))) (make-test-data (+ start 1) len)))))
(define (main args)
(print "test data:\n" (make-test-data 0 64) "\n")
(print "FFT result:\n" (fft (make-test-data 0 64)) "\n")
0)
Scipy(Python):
#!/usr/bin/env python
#coding: utf-8
from scipy.fftpack import fft
import cmath
def makeTestData(length):
pi = (cmath.atan(1)*4).real
return [(0.25*cmath.sin((i*pi)/16.0)).real for i in range(length)]
def printOrg(args):
for i in args:
print i
def main():
testdata = makeTestData(64)
printOrg(("test data:", testdata, "\n"))
printOrg(("FFT result:",list(fft(testdata))))
if __name__=="__main__":
main()
FFTW3(C):
#include "iostream"
#include "fstream"
#include "cmath"
#include <fftw3.h>
using namespace std;
int main(void){
/* size of array */
int N=64;
/* inut data */
double in[N];
/* output data */
fftw_complex *out;
/* plan */
fftw_plan p;
/* fike stream */
ofstream fout;
/* circular constant */
double pi;
pi=4.*atan(1.0);
/*--------------------------------------------
set initial data
---------------------------------------------*/
for(int i=0;i<N;i++)
in[i]=0.25*sin(pi*i/16);
printf("test data\n");
for(int i=0;i<N;i+=4){
for(int j=0;j<4;j++)
printf("%.17lf,",in[i+j]);
printf("\n");
}
printf("\n");
/*--------------------------------------------
make plan for out.dat
-------------------------------------------*/
out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N);
p=fftw_plan_dft_r2c_1d(N,in,out,FFTW_ESTIMATE);
/*--------------------------------------------
FFT execute and output out.dat
-------------------------------------------*/
fftw_execute(p);
fout.open("out.dat",ios::out);
for(int i=0;i<N;i++)
fout << out[i][0] << "+" << out[i][1] << "j, ";
fout << endl;
fout.close();
/*--------------------------------------------
finalize for out.dat
---------------------------------------------*/
fftw_destroy_plan(p);
fftw_free(out);
/*--------------------------------------------
myself test
---------------------------------------------*/
return 0;
}コンパイル例は、
g++ sinFFT.c -lfftw3 -lm -o sinFFT
です。
実験結果:
データは、
0.25*sin(i*pi/16) [0<=i<N] ※pi=atan(1)*4
で作りました。
64点の比較:
Gaucheの結果
test data: (0.0 0.04877258050403206 0.09567085809127245 0.13889255825490054 0.17677669529663687 0.2078674030756363 0.23096988312782168 0.2451963201008076 0.25 0.2451963201008076 0.23096988312782168 0.20786740307563636 0.1767766952966369 0.13889255825490054 0.09567085809127247 0.04877258050403215 3.061616997868383e-17 -0.04877258050403209 -0.09567085809127242 -0.1388925582549005 -0.17677669529663687 -0.2078674030756363 -0.23096988312782163 -0.24519632010080758 -0.25 -0.2451963201008076 -0.23096988312782166 -0.20786740307563636 -0.17677669529663692 -0.13889255825490054 -0.0956708580912726 -0.04877258050403218 -6.123233995736766e-17 0.04877258050403206 0.09567085809127249 0.13889255825490046 0.17677669529663684 0.2078674030756363 0.23096988312782163 0.24519632010080758 0.25 0.24519632010080766 0.23096988312782168 0.2078674030756364 0.1767766952966371 0.13889255825490057 0.09567085809127263 0.04877258050403199 9.184850993605148e-17 -0.04877258050403181 -0.09567085809127246 -0.13889255825490043 -0.17677669529663698 -0.20786740307563628 -0.2309698831278216 -0.24519632010080764 -0.25 -0.24519632010080766 -0.23096988312782168 -0.2078674030756364 -0.17677669529663712 -0.1388925582549006 -0.09567085809127265 -0.04877258050403202) FFT result: (5.721188726109846e-18 -1.8421947292037886e-16-2.9166809982044777e-16i -8.498308346471968e-16+8.0i 4.043583508252498e-16+8.076655356112082e-16i 7.608905701446265e-17+4.563343961201419e-16i 6.297497145813267e-17+1.5613079972381614e-16i -1.3970063801681657e-16-1.5616341699250225e-16i 2.330215154791825e-16-5.107134458601434e-17i 3.617815678165878e-16+4.0204732249633977e-16i 5.717554844058287e-16-1.966160023663208e-16i 4.938456040493896e-18+3.1821620729113183e-16i -1.7440997361919716e-16-2.7061971517407665e-16i 2.940044844704947e-16+3.8896004591581506e-16i 6.55154332450006e-16-3.4587607754684564e-16i 1.2289717057143745e-15-1.2375865681495943e-15i -5.969746088959934e-16-3.477681435157823e-16i -1.0530111373640578e-16+2.2204460492503146e-16i 3.674178734110533e-16-2.8814604803755036e-16i 5.416072942455187e-16+7.485047969070208e-16i -8.506982922946665e-17-1.0194425727101134e-16i 6.397406424680694e-17-3.56306875512855e-18i 3.3483894061684197e-16+9.428000023306903e-17i -1.0205118296577502e-16-6.243168134268933e-16i -2.6642515209083876e-16+4.070424721724977e-17i 2.0477232194821014e-16+6.898041510879305e-17i 7.893552629494091e-16-9.085792234317277e-17i 1.3527069786992498e-17-1.2809365211637865e-15i -8.518397658450753e-16-1.8206917554502449e-16i -1.891382459022937e-16+6.381128144919853e-17i 1.4385169300244952e-16-1.4136105444788473e-16i -7.267471740958727e-17+4.440892098500626e-15i -4.240721826793214e-16-1.9294112024003635e-17i -4.3836802112395277e-16 -4.2407218267932135e-16+1.929411202400329e-17i -1.8369701987210292e-16-8.881784197001252e-16i 1.4385169300244972e-16+1.4136105444788513e-16i -1.891382459022936e-16-6.381128144919829e-17i -8.51839765845075e-16+1.8206917554502446e-16i -4.255463544398187e-16+1.7271435842813282e-15i 7.89355262949409e-16+9.085792234317274e-17i 2.047723219482102e-16-6.89804151087929e-17i -2.664251520908385e-16-4.070424721724962e-17i 4.939773846683397e-16+5.87092241122261e-16i 3.3483894061684178e-16-9.428000023306925e-17i 6.397406424680686e-17+3.563068755128624e-18i -8.506982922946655e-17+1.0194425727101139e-16i -7.081873257584551e-16-1.871037900800844e-15i 3.6741787341105324e-16+2.8814604803755016e-16i -1.0530111373640583e-16-2.2204460492503146e-16i -5.969746088959934e-16+3.477681435157826e-16i -2.4286751921463073e-16+1.3967362279310442e-16i 6.551543324500061e-16+3.458760775468456e-16i 2.940044844704947e-16-3.8896004591581516e-16i -1.7440997361919716e-16+2.706197151740768e-16i -6.748990406600161e-17+8.296934855383177e-16i 5.7175548440582905e-16+1.9661600236632065e-16i 3.6178156781658783e-16-4.020473224963399e-16i 2.330215154791824e-16+5.107134458601419e-17i 5.291258494160125e-16-1.4007287666498568e-15i 6.297497145813317e-17-1.5613079972381606e-16i 7.608905701446285e-17-4.563343961201422e-16i 4.043583508252498e-16-8.076655356112086e-16i -2.959254581434994e-15-8.000000000000004i -1.8421947292037874e-16+2.9166809982044807e-16i)
Scipy(Python)の結果
test data: [0.0, 0.048772580504032062, 0.095670858091272445, 0.13889255825490054, 0.17677669529663687, 0.20786740307563631, 0.23096988312782168, 0.24519632010080761, 0.25, 0.24519632010080761, 0.23096988312782168, 0.20786740307563636, 0.17677669529663689, 0.13889255825490054, 0.095670858091272473, 0.048772580504032152, 3.061616997868383e-17, -0.04877258050403209, -0.095670858091272418, -0.13889255825490049, -0.17677669529663687, -0.20786740307563631, -0.23096988312782163, -0.24519632010080758, -0.25, -0.24519632010080761, -0.23096988312782166, -0.20786740307563636, -0.17677669529663692, -0.13889255825490054, -0.095670858091272598, -0.04877258050403218, -6.123233995736766e-17, 0.048772580504032062, 0.095670858091272487, 0.13889255825490046, 0.17677669529663684, 0.20786740307563631, 0.23096988312782163, 0.24519632010080758, 0.25, 0.24519632010080766, 0.23096988312782168, 0.20786740307563639, 0.17677669529663709, 0.13889255825490057, 0.095670858091272626, 0.048772580504031993, 9.1848509936051484e-17, -0.048772580504031812, -0.095670858091272459, -0.13889255825490043, -0.17677669529663698, -0.20786740307563628, -0.2309698831278216, -0.24519632010080764, -0.25, -0.24519632010080766, -0.23096988312782168, -0.20786740307563639, -0.17677669529663712, -0.1388925582549006, -0.095670858091272654, -0.04877258050403202] FFT result: [(5.7211887261098457e-18+0j), (-1.8421947292037881e-16+2.9166809982044777e-16j), (-1.7189608386965301e-15-8j), (4.0435835082524976e-16-8.0766553561120823e-16j), (7.6089057014462675e-17-4.5633439612014192e-16j), (6.2974971458132724e-17-1.5613079972381606e-16j), (2.5452196551807334e-16-6.0050017912816233e-16j), (2.3302151547918262e-16+5.1071344586014254e-17j), (3.6178156781658773e-16-4.0204732249633977e-16j), (5.7175548440582865e-16+1.9661600236632067e-16j), (1.0306143961226947e-16+2.8158535896252559e-16j), (-1.7440997361919682e-16+2.7061971517407675e-16j), (2.9400448447049455e-16-3.8896004591581501e-16j), (6.5515433245000576e-16+3.4587607754684545e-16j), (5.7396000703720268e-17+2.2204460492503131e-16j), (-5.9697460889599298e-16+3.4776814351578233e-16j), (-1.0530111373640581e-16-2.2204460492503131e-16j), (3.67417873411053e-16+2.8814604803755011e-16j), (5.7396000703720268e-17-6.6613381477509392e-16j), (-8.5069829229466307e-17+1.0194425727101136e-16j), (6.3974064246806987e-17+3.5630687551287716e-18j), (3.3483894061684173e-16-9.4280000233069054e-17j), (1.8548887372774134e-16+4.8649361844068378e-16j), (-2.6642515209083851e-16-4.0704247217249745e-17j), (2.0477232194821023e-16-6.8980415108792806e-17j), (7.8935526294940887e-16+9.0857922343172754e-17j), (-2.1786391399475917e-17+1.3807649197502463e-15j), (-8.5183976584507503e-16+1.8206917554502459e-16j), (-1.891382459022936e-16-6.3811281449198188e-17j), (1.4385169300244957e-16+1.4136105444788528e-16j), (-3.8669320914634225e-16+0j), (-4.240721826793214e-16+1.9294112024003339e-17j), (-4.3836802112395277e-16+0j), (-4.240721826793214e-16-1.9294112024003339e-17j), (-3.8669320914634225e-16-0j), (1.4385169300244957e-16-1.4136105444788528e-16j), (-1.891382459022936e-16+6.3811281449198188e-17j), (-8.5183976584507503e-16-1.8206917554502459e-16j), (-2.1786391399475917e-17-1.3807649197502463e-15j), (7.8935526294940887e-16-9.0857922343172754e-17j), (2.0477232194821023e-16+6.8980415108792806e-17j), (-2.6642515209083851e-16+4.0704247217249745e-17j), (1.8548887372774134e-16-4.8649361844068378e-16j), (3.3483894061684173e-16+9.4280000233069054e-17j), (6.3974064246806987e-17-3.5630687551287716e-18j), (-8.5069829229466307e-17-1.0194425727101136e-16j), (5.7396000703720268e-17+6.6613381477509392e-16j), (3.67417873411053e-16-2.8814604803755011e-16j), (-1.0530111373640581e-16+2.2204460492503131e-16j), (-5.9697460889599298e-16-3.4776814351578233e-16j), (5.7396000703720268e-17-2.2204460492503131e-16j), (6.5515433245000576e-16-3.4587607754684545e-16j), (2.9400448447049455e-16+3.8896004591581501e-16j), (-1.7440997361919682e-16-2.7061971517407675e-16j), (1.0306143961226947e-16-2.8158535896252559e-16j), (5.7175548440582865e-16-1.9661600236632067e-16j), (3.6178156781658773e-16+4.0204732249633977e-16j), (2.3302151547918262e-16-5.1071344586014254e-17j), (2.5452196551807334e-16+6.0050017912816233e-16j), (6.2974971458132724e-17+1.5613079972381606e-16j), (7.6089057014462675e-17+4.5633439612014192e-16j), (4.0435835082524976e-16+8.0766553561120823e-16j), (-1.7189608386965301e-15+8j), (-1.8421947292037881e-16-2.9166809982044777e-16j)]
FFTW3(C)の結果
test data 0.00000000000000000,0.04877258050403206,0.09567085809127245,0.13889255825490054, 0.17677669529663687,0.20786740307563631,0.23096988312782168,0.24519632010080761, 0.25000000000000000,0.24519632010080761,0.23096988312782168,0.20786740307563636, 0.17677669529663689,0.13889255825490054,0.09567085809127247,0.04877258050403215, 0.00000000000000003,-0.04877258050403209,-0.09567085809127242,-0.13889255825490049, -0.17677669529663687,-0.20786740307563631,-0.23096988312782163,-0.24519632010080758, -0.25000000000000000,-0.24519632010080761,-0.23096988312782166,-0.20786740307563636, -0.17677669529663692,-0.13889255825490054,-0.09567085809127260,-0.04877258050403218, -0.00000000000000006,0.04877258050403206,0.09567085809127249,0.13889255825490046, 0.17677669529663684,0.20786740307563631,0.23096988312782163,0.24519632010080758, 0.25000000000000000,0.24519632010080766,0.23096988312782168,0.20786740307563639, 0.17677669529663709,0.13889255825490057,0.09567085809127263,0.04877258050403199, 0.00000000000000009,-0.04877258050403181,-0.09567085809127246,-0.13889255825490043, -0.17677669529663698,-0.20786740307563628,-0.23096988312782160,-0.24519632010080764, -0.25000000000000000,-0.24519632010080766,-0.23096988312782168,-0.20786740307563639, -0.17677669529663712,-0.13889255825490060,-0.09567085809127265,-0.04877258050403202 FFT result: 5.72119e-18+0j, -1.84219e-16+2.91668e-16j, -1.88549e-15+-8j, 4.04358e-16+-8.07666e-16j, 7.60891e-17+-4.56334e-16j, 6.2975e-17+-1.56131e-16j, 3.45995e-16+-7.46308e-16j, 2.33022e-16+5.10713e-17j, 3.61782e-16+-4.02047e-16j, 5.71755e-16+1.96616e-16j, 7.37882e-18+2.95025e-16j, -1.7441e-16+2.7062e-16j, 2.94004e-16+-3.8896e-16j, 6.55154e-16+3.45876e-16j, 1.12907e-16+2.22045e-16j, -5.96975e-16+3.47768e-16j, -1.05301e-16+-2.22045e-16j, 3.67418e-16+2.88146e-16j, 1.12907e-16+-6.66134e-16j, -8.50698e-17+1.01944e-16j, 6.39741e-17+3.56307e-18j, 3.34839e-16+-9.428e-17j, 2.53264e-16+6.28092e-16j, -2.66425e-16+-4.07042e-17j, 2.04772e-16+-6.89804e-17j, 7.89355e-16+9.08579e-17j, -8.53516e-17+1.36312e-15j, -8.5184e-16+1.82069e-16j, -1.89138e-16+-6.38113e-17j, 1.43852e-16+1.41361e-16j, -3.31182e-16+0j, -4.24072e-16+1.92941e-17j, -4.38368e-16+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j
まあ、それぞれの結果格納様式は別として。
部分的に取り出してみました。
Gauche: 5.721188726109846e-18 -1.8421947292037886e-16-2.9166809982044777e-16i -8.498308346471968e-16+8.0i 4.043583508252498e-16+8.076655356112082e-16i Scipy(Python): (5.7211887261098457e-18+0j), (-1.8421947292037881e-16+2.9166809982044777e-16j), (-1.7189608386965301e-15-8j), (4.0435835082524976e-16-8.0766553561120823e-16j) FFTW3(C): 5.72119e-18+0j, -1.84219e-16+2.91668e-16j, -1.88549e-15+-8j, 4.04358e-16+-8.07666e-16j
以上より、ほぼ同じ結果ですね。
ちなみに、67点でのScipyとFFTW3の結果から部分的に取り出したものは、
Scipy(Python): (0.14444343859530409+0j), (0.20785712096923026-0.43878172487608225j), (2.3808744181020063-7.7202925540628877j), (-0.27391453637245922+0.95925253678366906j) FFTW3(C): 0.144443+0j, 0.207857+-0.438782j, 2.38087+-7.72029j, -0.273915+0.959253j
という結果で、これもほぼ同じです。
なので。
やっぱり、scipy.fftpackのfftはfftwを使っているのかもしれません。
おわりに:
今回は、ちょっと駆け足でした。
思いつきというか手遊び気分だったんで、何かすみません。
何か、いいことあったらいいな。ないから困ってるんだけど。
ということで。
また、次の機会に。
10000ビュー越えた？
気がつけば、2ヶ月くらい更新してませんでした。
くらいって言うか、まあちょうど2ヶ月ですが。
で、気がつけば、ページビューが10000を越えてました。
ほとんど自分だけどね☆
更新のまっっっっったくないこんな辺境でも、それくらいはなるもんなんですねぇ。
いやはや。
一万回も恥を晒しているのかと思うと、ゾックゾクす(ry
………………。
しませんよ？ 念のため。
まあ、色々無い状況で色々とありまして、現在他のサイトの手入れも放置中ですが、その内更新すると思います。
願わくば、後進ではなく更新であって欲しいものです。（他人事のように）
Pythonのなみのり！ ぱーと2
前回、と言っても相当前ですけど、Pythonでwavファイルを読み込むってのを記事にしました。
今回はその焼き直しですが、かなり仕様みたいのは買えました。
いちいちキー指定なんてしてらんねーよ。
と、自分でも思ったので。
オブジェクト指向的なプログラミングを意識したつもりなんですが、そもそもオブジェクト指向を理解していないので、全然違う実装になっていると思います。
ソースは以下。
struct.unpackのくだりが遅いのは相変わらず。
格納をリストではなく文字列にしたり、ファイルの読み込みの段階からunpackしたりして実効速度の向上を図ってみましたが、結局この方法に落ち着きました。
コメントの英語に対するツッコミはご遠慮願います。
ソースコード：
#!/usr/bin/env python #coding: utf-8 import struct class wavUtils(object): """ This class is utilities for wav file. At first, you make instance this class. Then, use "load" method for reading a wav file. """ def __init__(self): self.__raw = None self.__raw_data = None self.__raw_header = None self.__unpacked_header = None self.__unpacked_data = None self.__meta_data = None def getDataList(self): return self.__unpacked_data def getHeaderList(self): return self.__unpacked_header def getPropertyList(self): return self.__meta_data def getRawStr(self): return self.__raw def getRawDataStr(self): return self.__raw_data def getRawHeaderStr(self): return self.__raw_header def load(self, fname): """ At first, loading wav file in self.__raw. Then, separates self.__raw for header part and date part. Last, unpack binary data for int. Unpacking data needs time. """ fp = open(fname, "rb") self.__raw = fp.read() fp.close() # separate header part and data part. header_end = self.__raw.find('data')+8 self.__raw_header = self.__raw[:header_end] self.__raw_data = self.__raw[header_end:] # unpacked header and data self.__unpackRawHeader() self.__unpackRawData() def __unpackRawData(self): seek_head = 0 tmp_result = [] # counter tool's values cnt = 1 cnt_delta = len(self.__raw_data)/10 # 16 bit/sample? 8 bit/sample? if self.__unpacked_header[7] == 8: unpack_fmt = '<H' seek_delta = 1 else: unpack_fmt = '<h' seek_delta = 2 print "Now converting that binary for Integer..." while seek_head < len(self.__raw_data): tmp_result.append(struct.unpack(unpack_fmt, self.__raw_data[seek_head:seek_head+seek_delta])[0]) seek_head += seek_delta # counter tool's implements if seek_head >= (cnt_delta*cnt): for i in range(cnt): print ".", if cnt==10: print cnt*10, "% done. Converting finish!" else: print cnt*10, "% done." cnt+=1 self.__unpacked_data = tmp_result def __unpackRawHeader(self): """ This method unpack header from binary for int. And after unpacking, this makes property list. """ tmp_result = [] for current_chunk in ('RIFF','fmt','data'): tmp_value = 0 seek_head = self.__raw_header.find(current_chunk) if current_chunk == 'fmt': seek_head += 4 for seek_delta in (4,2,2,4,4,2,2): if seek_delta == 4: tmp_value = struct.unpack('<I',self.__raw_header[seek_head:seek_head+seek_delta])[0] else: tmp_value = struct.unpack('<H',self.__raw_header[seek_head:seek_head+seek_delta])[0] seek_head += seek_delta tmp_result.append(tmp_value) else: tmp_value = struct.unpack('<I',self.__raw_header[seek_head+4:seek_head+8])[0] tmp_result.append(tmp_value) self.__unpacked_header = tmp_result # Make property list. tmp_result = [] wav_head_ref = ("File size (minus 8 byte).", "fmt chunk size", "FormatID", "channel", "sampling rate", "Byte/sec", "block size(Byte/sample*channel)", "bit number(bit/sample)", "data size") for i in range(len(self.__unpacked_header)): tmp_result.append((wav_head_ref[i], self.__unpacked_header[i])) # add information of "bit rate". tmp_result.append(("bit rate (Byte/sec*8)", self.__unpacked_header[5]*8)) self.__meta_data = tmp_result def main(): wu = wavUtils() wu.load("sample.wav") print wu.getRawStr()[0:5], wu.getRawHeaderStr(), wu.getHeaderList()[7] print wu.getDataList()[100000:100100] print wu.getPropertyList() if __name__ == "__main__": main()
あとがき：
まあ、前回掲載したものよりは、使いやすいと思います。
ソースも、心なしすっきりしたように見えます。
ああ、気のせいですか、そうですか。
オブジェクト指向って、こういうのじゃないのかなぁ……
