So-net無料ブログ作成
検索選択

「A SMARTERWAY TO FIND PITCH」でピッチ検出 [Math]

音声の波形からピッチを検出するアルゴリズム にピッチ検出の論文が
解説されていて、 scilab でやってみました。

Normalized Square Difference なるものを計算して、
次のアルゴリズムで音程を決定するとのこと。
1.極大値をリストアップ
2.その中から傾きが正で0を跨ぐ部分から、負で0を跨ぐ部分の間で最大のものに絞る。
ただし、傾き生で0を跨ぐが、途中で切れているものも候補に入れておく
3.最大値を nmax として、k * nmax のもので、最初にでてくるもの(周波数が小さいもの)
を選ぶ。k は 0.8~1.0 で要調節。

440Hzの正弦波のNormalized Square Difference を描いたものがこちら。
計算あってるのか?-1~1 になっているからたぶんあってる。

sinwave_nsdf.png

上記のアルゴリズムを適用すると、たしかに1/440=0.00227 が取れそうな感じだ。
もっといろんな音で確かめたい。

使った scilab コードはこちら。


nice!(0)  コメント(0)  トラックバック(0) 
共通テーマ:音楽

scilab で自己相関関数 [Math]

信号 x = {...} があったらその自己相関関数を求めるには、
ざっくりと以下のようになる。

パワースペクトルは元信号をフーリエ変換して求める。
P=|F(x)|^2

自己相関関数のフーリエ変換でもパワースペクトルが求まる
(ウィーナー=ヒンチンの定理と言うそうな)。
P=F(r)/N

なのでパワースペクトルをフーリエ逆変換すると自己相関関数が得られる。
自己相関をは直接求めると計算コストがかかるので、x を fft して P を
求め、それをフーリエ逆変換して求めるのが一般的。

みたいなことが、教科書とかに書いてある。あ、そうなんですか…と。
ピンとこないので scilab で見てみました。

sinwave_corr.png

fft+ifft で求めた自己相関関数と、ごりごり計算した自己相関関数が
ぴったりと重なりました。fft+ifft の方はパワースペクトルと同じく
真ん中で対象になっているので、データを半分にして表示しました。

使った scilab スクリプト
funcprot(0);
clear all;

//////////////////////////////////////////////////////////////////////
// parameters
//////////////////////////////////////////////////////////////////////
sampleNum = 1024;   // sample number
sampleFreq = 44100; // sampling frequency
waveFreq = 440;     // wave frequency

//////////////////////////////////////////////////////////////////////
// main
//////////////////////////////////////////////////////////////////////
// create time-axis
ts = [0:sampleNum-1] * 1/sampleFreq;
// create sin wave
x=sin(2 * %pi * waveFreq * ts);
x(sampleNum/2+1:sampleNum)=zeros(1,sampleNum/2);
// create freq-axis
fs = [0:sampleNum-1] * sampleFreq/sampleNum;
// data number for fft
N = sampleNum;
// fft
f_fft = fft(x);
power = f_fft.*conj(f_fft);
// correlation by fft
r=ifft(power);
taus=[0:N-1] * 1/sampleFreq;
// correlation by hand
M=N/2;
r2=zeros(1,M);
taus2=[0:M-1] * 1/sampleFreq;
for m=1:M+1
    for n=1:M-m+1
        r2(m)=r2(m)+x(n)*x(n+m-1);
    end;
end;

subplot(211);
plot2d(ts,x);
title = sprintf("A %dHz wave is sampled with sampling rate=%dHz",waveFreq,sampleFreq);
xtitle(title,"time (sec)","arbitral");

subplot(212);
plot(taus(1:M),r(1:M),taus2,r2);
xtitle("autocorrelation computed by two different way", "tau(s)", "auto correlation");
legend(["fft" "straight"], 1, "ur");



ピッチ検出に必要なバッファサイズ [Math]

ピアノの最低音の周波数は27.5Hzとのこと。
( http://d.hatena.ne.jp/session_oyaji/20070520/1179666560 )
これの周期を図れるくらいバッファに貯めると
どんな感じになるのかを scilab で見てみました。

44.1KHz で 4096 個サンプリングすると、
275Hz, 55Hz, 110Hz,... はこんな感じになります。

低い周波数を対象にすると結構長い事サンプリングが必要だなーと。
サンプリング周波数を下げれば良いですが、同時に高い周波数も
対象にしたいですし。悩ましい。
かといってまだ何も実装していないので、4096 このサンプリングが
どれほどの負荷なのか知りませんけど。時間にして 0.1 秒くらいは
サンプリングしないとダメなんですね。

harmonics.png

使った scilab のコードは以下。
funcprot(0);
clear all;

//////////////////////////////////////////////////////////////////////
// parameters
//////////////////////////////////////////////////////////////////////
sampleNum = 4096;   // sample number
sampleFreq = 44100; // sampling frequency
baseWaveFreq = 27.5;// base wave frequency

//////////////////////////////////////////////////////////////////////
// main
//////////////////////////////////////////////////////////////////////
// create time-axis
ts = [0:sampleNum-1] * 1/sampleFreq;
// create sin wave of frequency 
// 27.5, 55, 110, 220, 440, 880, 1760, 3520 Hz
powNum = 4;
pows=[0:powNum-1]';
fs=baseWaveFreq*(2^pows);
waves = sin(2 * %pi * fs * ts);

// plot
for i=1:powNum,
    subplot(int(powNum/2)*100+20+i)
    plot2d(ts,waves(i,:));
    legendtext = sprintf("%f Hz",fs(i));
    legends(legendtext,1,"ur");
end


scilab で fft [Math]

前回の正弦波を fft してみた。
fft の分解能は サンプリング周波数/サンプルの数 だそうで、
今回は 44100 Hz/1024=約43.1Hz
なので、 430 Hz のデータが最大なので、そこから + 43.1 した
430 〜 473 Hz の範囲が最大なんじゃないかと。であれば 440Hz
はこの範囲なのであってるなと。

sinwave_fft.png

scilab スクリプトはこんな感じ。
funcprot(0);
clear all;

//////////////////////////////////////////////////////////////////////
// parameters
//////////////////////////////////////////////////////////////////////
sampleNum = 1024;   // sample number
sampleFreq = 44100; // sampling frequency
waveFreq = 440;     // wave frequency

//////////////////////////////////////////////////////////////////////
// main
//////////////////////////////////////////////////////////////////////
// create time-axis
ts = [0:sampleNum] * 1/sampleFreq;
// create sin wave
f=sin(2 * %pi * waveFreq * ts);
// data number of data for fft
N = sampleNum/2;
// create freq-axis
fs = [0:sampleNum] * sampleFreq/sampleNum;
// fft
f_fft = fft(f);
power = f_fft.*conj(f_fft);

subplot(211);
plot2d(ts,f);
title = sprintf("A %dHz wave is sampled with sampling rate=%dHz",waveFreq,sampleFreq);
xtitle(title,"time (sec)","arbitral");

subplot(212);
plot2d("nl",fs(1:N),power(1:N));
xtitle("powerspectrum", "frequency(Hz)", "arbitral");

[maxPower,maxIndex] = max(abs(power(1:N)));
maxFreq=fs(maxIndex);
printf("max frequency range=%f~%f\n",maxFreq, maxFreq + sampleFreq/sampleNum);



scilab で正弦波を描く [Math]

scilab で音声の解析とかしたい。
久しぶりなのでリハビリを兼ねて正弦波を描く練習。

440Hz(ラの音)を44.1KHzでサンプリングしたものを1024個取り出して表示。
周期は 1/440=0.0023 だからたぶん合ってる。
440HzSignWave.png

scilab のスクリプトはこんな感じ。
funcprot(0);
clear all;

//////////////////////////////////////////////////////////////////////
// parameters
//////////////////////////////////////////////////////////////////////
sampleNum = 1024;   // sample number
sampleFreq = 44100; // sampling frequency
waveFreq = 440;     // wave frequency

//////////////////////////////////////////////////////////////////////
// main
//////////////////////////////////////////////////////////////////////
ts = [0:sampleNum-1] * 1/sampleFreq;
f=sin(2 * %pi * waveFreq * ts);
plot2d(ts,f);
title = sprintf("A %dHz wave is sampled with sampling rate=%dHz",waveFreq,sampleFreq);
xtitle(title,"time (sec)","arbitral");


二項分布から拡散方程式 [Math]

二項分布では n ステップ後に位置 x にいる確率は p(x,n) で表しました。
この確率から差分方程式を作って拡散方程式が導かれます。せっかく p(x,n)
前回計算したので、拡散方程式も見てみないとモッタイナイのでここに
書いておきます。なお、内容は検索で出てきた立命館大学の講義資料
まんまです。

今回 x は変数変換後の変数として使うので、変換前は m と書いておきます。
なので、n ステップ後に位置 m にいる確率を p(m,n) と今回は表記します。

ステップ n+1 で m にいる確率 p(m,n+1)は以下のようになる。


(ステップ n で m-1 にいる確率)*(右に移動する確率)と
(ステップ n で m+1 にいる確率)*(左に移動する確率の和。

これを変形して(こんな変形思いつかないけど、賢い人は変形できるんでしょうねぇ)


変数変換で極限をとって連続にする。




連続版の確率密度分布を f(x,t) とすると、区間 Δx の確率は次のように表せる。


差分方程式を f(x,t) で表すと

微分の形に変形して

こんな変形絶対できないんですけど。まぁ、天才アインシュタイン
考えてノーベル賞を取った理論だそうなので、できなくて当然。

を保ったまま
そんなこと言われても、ああそうですかとしか。

最終的に得たのが次の式で、Dを拡散係数と呼ぶ。


元の講義ノートから変数の文字を変えたりして結局見にくくなった感が。。。
まぁこれを書いてわかった気になるのがひとつの目的なので。

二項分布でランダムウォーク [Math]

何かと出てくるランダウウォークですが、Scilab でシミュレーションして
理解した気になりたいところですが、こればかりは式の導出を避けては通れ
なさそうなので、導出をトレースしてわかった気になろうと思います。

原点からステップ毎に確率1/2で -1 or 1 を選んで移動する現象があるとします。
ステップ回数が多いとだんだん原点から離れていきます。この現象の何が大事
なのかというとこのサイトにとてもわかりやすくまとまっていて次の3つが
挙げられています。

1. n ステップ後に x にいる確率は正規分布になる
2. 分布の幅が√n で広がっていく
3. 原点にいる確率は 1/√n で減っていく

この3つを見ておきたいと思います。

式の導出は一応自分でも手計算で確認しつつ、こちらにまとめました。

nステップ後に位置mにいる確率



式変形により正規分布の形に近似した式

a はステップ毎の移動距離。xは位置(x=a*m)。

平均0、分散μの正規分布


正規分布と比べると

と対比できるので、「分布の幅が√n で広がっていく」のがわかります。

また、原点にいる確率は

これから「原点にいる確率は 1/√n で減っていく」というのも見て取れます。
なるほど。


スターリングの公式 [Math]

最近個人的な理由で、ランダムウォークを理解するため、
二項分布から正規分布を求める式変形を理解しようとしています。
その中で出てくるスターリングの公式です。

ここで書くまでもなく、ググれば良質な情報がいくらでも出てきますが、
理解した気になるために、ここに書いておきます。



両辺の対数をとる



log(1+x)をマクローリン展開 [Math]

二項分布から正規分布を求めようとしています。
その中で log(1+x) のマクローリン展開が出てくるのでメモっておきます。

マクローリン展開はテイラー展開の x=0 の周りで展開するバージョン

テイラー展開


マクローリン展開は a=0 なので


n=3 以降は O(3) とかやって捨てるので n=2 までが必要です。
そのためには x=0 のときの log(1+x) の微分と二回微分が出てくるので計算しておきます。

log(1+x)の微分

x=0 のときは 1

log(1+x)の二回微分

x=0 のときは -1


あと log(1+x) の x=0 のときのゼロ階微分は 0。
これらを使って、log(1+x) をマクローリン展開



3乗の項以降は O(3) と書いて、最終的には捨てました。

Scilab でコイン投げの二項分布 [Math]

Scialb でシミュレートしてわかった気になるシリーズ第三弾。
コイン投げを、n 回投げて、 x 回表が出る確率の分布が二項分布。

1.n回コインを投げる
2.表が出た回数を数える
....というのを1000回繰り返して分布を描く。

n=5
binomialx_5.png

n=100
binomialx_50.png

n が大きくなるにつれて、平均=n*p、分散=n*p(1-p) ただしp=1/2 の正規分布に近づいている。
なんか嘘っぽいけど。

Scilab のコード。Github にもおいてあります。

clear;

// sampling num
samplenum=100
// number of repetition
repeat=1000;

// coin toss
R=grand(repeat,samplenum,'uin',0,1);

scf(1);
sums=sum(R,'c');
histplot([0:1:samplenum],sums)

n=samplenum;
p=1/2;

avg=n*p;
s=sqrt(n*p*(1-p))
x=[0:0.1:samplenum];
plot(x,(1/(s*sqrt(2*%pi)))*exp(-((avg-x).^2)./(2*s^2)));



この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。