自己相関関数のxcorr

Question

scilabで自己相関関数を求めたいです.
私が実現したいのはx1(t)とそこからlだけ遅れたx1(t+l)の自己相関関数を求めてそれをr(j)に代入するということです.
scilabのxcorr関数で実現できないかということでやってみているのですが,自己相関関数の定義がうまくできません.

//---------------------------------------------------------
// 分析用パラメータ設定
//---------------------------------------------------------
len=512;
order=200;
sfreq=16;
fshift=10;
fwidth=len/sfreq;


[x,y]=loadwave('/Users/1.wav'); 
n_sample=y(8);
disp(n_sample); // 音声データのサンプル数をコンソールに表示
for i=1:n_sample
    ts(i)=i/sfreq;
end


//---------------------------------------------------------
// Hamming窓を生成
//---------------------------------------------------------
for i=1:len
    win(i)=0.54-0.46*cos(2*%pi*i/len);
end
n_frame=0;

//---------------------------------------------------------
// 短時間対数パワーを計算
//---------------------------------------------------------
for st=1:sfreq*fshift:n_sample-len
    pwr=0;
    for i=1:len
        x1(i)=x(i+st)*win(i);
        pwr=pwr+x1(i)*x1(i);
    end
    n_frame=n_frame+1;
    tp(n_frame)=fwidth/2+(n_frame-1)*fshift;
    logpwr(n_frame)=10*log10(pwr);
end

//---------------------------------------------------------
// 短時間対数パワーを描く
//---------------------------------------------------------
subplot(3,1,2);
xlabel("Time [ms]");
ylabel("Speech Power [dB]");
plot2d(tp,logpwr);


//---------------------------------------------------------
// 基本周波数の推定（自己相関法）
//---------------------------------------------------------
thlogpwr=-25;
for n=1:n_frame // n_frameは分析フレーム数
    if logpwr(n) < thlogpwr then
        f0(n)=0;
    else
        st=1+(n-1)*sfreq*fshift; // 現フレームの開始サンプル番号
        j=1;
        for l=0:order
           
               [r,len]=xcorr(x1,"biased"); //自己相関関数
           
 
            //
        end
        //max=0;
        //for j=80:200
           
           argmax=max(r(80:200));
            
        period = argmax/16;
        f0(n) = 1000/period;
    end
end
subplot(3,1,3);
xlabel("Time [ms]");
ylabel("Fundamental Frequency [Hz]");
plot2d(tp,f0);

またこのxcorrで求めた自己相関関数でこのようなグラフが出力されたのですが私が求めたいのはこのような台形ではなく,ピークが局所的に目立つようにしたいです.

ピークが目立つようにとはこの図の右の赤いグラフのような形です.
![説明]

質問しましたが,もう少し勉強してから質問修正します.

Comments

[jbpb0]

> x1(t)とそこからlだけ遅れたx1(t+l)の自己相関関数を求めてそれをr(j)に代入する

の「l」と「j」は、どのような関係ですか？

for n=1:n_frame
で「n」を変えながらループを回していますが、その中で使われてる「x1」は毎回同じものですけど、それは意図通りですか？

そこよりもずっと上の
for st=1:sfreq*fshift:n_sample-len
のループ内で、ループの毎回
for i=1:len
x1(i)=...
でx1に値を上書き代入しているので、x1にはstのループの最後の回の計算値だけが入ってますが、それは意図通りですか？

[kaeruuuun]

基本的にxcorrで自己相関関数を定義したところ以外は意図通りです.
音声信号を読み込んで分析窓関数(同じものを繰り返し使うので先に計算)x1を生成しています.
そのx1(t)を用いてx1(t)とx1(t+j)の自己相関関数を求めたいわけです.
>「l」と「j」は、どのような関係ですか？
申し訳ございません.
私の勘違いで上に書いたように自己相関関数r(j)をx(t)とx(t+j)から求めたいということです.

[jbpb0]

> 音声信号を読み込んで分析窓関数(同じものを繰り返し使うので先に計算)x1を生成しています.
そのx1(t)を用いてx1(t)とx1(t+j)の自己相関関数を求めたいわけです.

for n=1:n_frame
のループで、「n」を変えながらたくさんの回数xcorr()を計算してる理由がわかりません
だって、x1は1セットしかなくて、その1セットをループで毎回使っているのだから、「n」をいろいろ変えながら計算しても、x1から得られる自己相関関数はやはり1セットしかないです
ループを回す必要ないですよね？
1回だけxcorr()で自己相関関数を計算したらいいのではないですか？

「n」によってx1が変わるのならわかりますが、コードはそうなってないので

[kaeruuuun]

でもjを1からorderまで変化させてそれぞれで自己相関関数を求めたいで,1回ではダメだと思います.

[jbpb0]

上にも書きましたが、
for st=1:sfreq*fshift:n_sample-len
のループ内で、x1には毎回上書き代入されるので、そのループが終わった時点では、「st=n_sample-len」という条件で計算された結果だけがx1に入ってます

具体的には、
for i=1:len
のループで「i」が1〜512(=len)変わるので、
i=1の時：x1(1)=x(1+n_sample-512)*win(1)=x(n_sample-511)*win(1)
〜(略)〜
i=512の時：x1(512)=x(512+n_sample-512)*win(512)=x(n_sample)*win(512)
だけがx1に入ってます

つまり、x1には、x(元データ)の中の、最後の512個のデータだけが反映されてます
それ以外のデータは、x1には関係ありません

> そのx1(t)を用いてx1(t)とx1(t+j)の自己相関関数を求めたい

本当に、xの最後の512個にしか関係してないx1の自己相関関数が計算できたら、それでいいのですか？
xの他の部分は、要らないのですか？

[jbpb0]

> jを1からorderまで変化させてそれぞれで自己相関関数を求めたいで,1回ではダメ

いろんな「j」に対する値が全部計算されたのが自己相関関数で、xcorr()はそれを1回で計算してくれる

一つの「j」だけで計算されたのは自己相関関数ではない
(元データを正規化してたら、それは自己相関係数)

https://teratail.com/questions/341675
にも書きましたけど、
https://help.scilab.org/docs/6.1.0/ja_JP/xcorr.html
の「例」の計算をそのまま自分でやって、元データ「y」の形(データのサイズ)や、計算結果「c」と「ind」の形を確認してください
そうしたら、私が言ってる意味がわかると思います

[kaeruuuun]

コメントありがとうございます.今すぐに理解が難しそうなので時間がある時にゆっくり確認してまた返事をさせていただきます.申し訳ございません.

[kaeruuuun]

少し勉強して考えたプログラムを書きます.
st=1+(n-1)*sfreq*fshift; // 現フレームの開始サンプル番号
j=1;
for l=0:order
//
r(j)=0;
// この部分を穴埋め
for t = 1:st-l
y(t+l)= circshift(x1(t),l);
[r(j),lags(j)]=xcorr(x1(t),y(t+l));
end
j=j+1;
end
x1とlだけ遅らせたyとの相関を求めれば良いのではないかと考えました.

実行されたファイルの 90 行目

xcorr: 引数 #2 の値が間違っています: 1 より大きい値を指定してください.

全て配列の引数は1から開始にしているのにこのようなエラーが表示されました.

追記
色々改良しているうちにわからなくなってきました.

[kaeruuuun]

>本当に、xの最後の512個にしか関係してないx1の自己相関関数が計算できたら、それでいいのですか？
xの他の部分は、要らないのですか？

今回分析したいのは512点なので,それで大丈夫です.

[kaeruuuun]

>いろんな「j」に対する値が全部計算されたのが自己相関関数で、xcorr()はそれを1回で計算してくれる

と言うことはfor文を回す必要はないのですか?

[jbpb0]

> for文を回す必要はないのですか?

何度も言ってますけど、
https://help.scilab.org/docs/6.1.0/ja_JP/xcorr.html
の「例」の計算をそのまま自分でやって、元データ「y」の形(データのサイズ)や、計算結果「c」と「ind」の形を確認してください
そうしたら分かります

上記「例」では、forループ無しでxcorr()の計算1回だけして、計算結果の「c」と「ind」でplot()してますよね
plot()結果のグラフもありますけど、いろんなタイムラグ(グラフ横軸)に対する自己相関関数の値が、ちゃんとグラフにプロットされてますよね

[kaeruuuun]

わかりました.matlabにある例でもやって見たりしてるのですが,理解が追いついていないです.勉強します.

[jbpb0]

> 今回分析したいのは512点なので,それで大丈夫です

それなら、何で
for n=1:n_frame
のループを回しているのですか？

「n」が何であろうと、x1は(512個のデータが入ってる)1セットしかないのだから、x1から自己相関関数を計算した結果はループの毎回で同じです
同じ結果にしかならないのだから、ループを回して何回も自己相関関数を計算する意味無いですよね

【追記】 for n=1:n_frameのループが終わった時点で、f0に入ってる値は2種類しかありません
グラフが台形になるのは、そのためです

1つは、「if logpwr(n) < thlogpwr」が成立した場合の「f0(n)=0」
もう一つは、上記if文が成立しない場合の「f0(n) = 1000/period」

後者は、これまでに何度も書いてるように、ループの毎回で同じx1から計算したrからperiodを(max()とか使って)計算してるので、periodはループのどの回も同じにしかならないので、f0は全部同じ値です

結果、f0は、0と、x1から計算した値の2種類だけです

質問者さんがやりたいことを実現するには、
https://teratail.com/questions/342569
にも書きましたけど、x1を計算してるところを直さないとダメです
そこ直さないで、xcorr()とかmax()とかのところだけ直しても、ダメです
(もちろん、xcorr()とかmax()とかのところも直さないといけないけど)

Answer 1

x1(t)とそこからlだけ遅れたx1(t+l)の自己相関関数を求めてそれをr(j)に代入する

は、

        st=1+(n-1)*sfreq*fshift; // 現フレームの開始サンプル番号
        j=1;
        for l=0:order

               [r,len]=xcorr(x1,"biased"); //自己相関関数


            //
        end

で、ちゃんと計算できてるはず (st, j, lはx1には無関係なので要らないけど、有っても計算結果には害はない)

lenには自己相関関数のタイムラグが入ってるので、そこから必要なところを探して、rからそこだけ取り出してください

たとえば、x1の要素数が512個の場合、lenには
-511, -510, -509...509, 510, 511
が入ってると思うので、欲しいのがタイムラグ0〜のものならば、rの先頭からではなく、512番目以降のを取り出すことになると思います
len(512)に0が入ってれば、r(512)がタイムラグ0の自己相関関数の値です

　.

このxcorrで求めた自己相関関数でこのようなグラフが出力されたのですが私が求めたいのはこのような台形ではなく

の原因は、x1からの自己相関関数の計算のやり方ではなく、x1自体を計算して求めるところのやり方にありますので、そこを直す必要があります

Comments

[kaeruuuun]

ありがとうございます.
>lenには自己相関関数のタイムラグが入ってるので、そこから必要なところを探して、rからそこだけ取り出してください

たとえば、x1の要素数が512個の場合、lenには
-511, -510, -509...509, 510, 511
が入ってると思うので、欲しいのがタイムラグ0〜のものならば、rの先頭からではなく、512番目以降のを取り出すことになると思います
len(512)に0が入ってれば、r(512)がタイムラグ0の自己相関関数の値です

これがよくわかりません.lenがx(t)とx(t+len)の自己相関関数を求めた時のlenですよね？

[kaeruuuun]

上のコメントは無視してください.
なんとなくわかりました.タイムラグとしてlenには負の値が入っているので,タイムラグが0からのが必要であれば,r(512)から取り出すということですね.ちなみにこれはいつの段階で行えばいいでしょうか.
rの最大値の位置を求める時なのか,x1の自己相関関数を求めたあとなのかがよくわかりません.
r(80)からr(200)の最大値を求めたければ,r(592)とr(712)にしたら良いということでしょうか.

[jbpb0]

> タイムラグとしてlenには負の値が入っているので,タイムラグが0からのが必要であれば,r(512)から取り出す

そうです

> いつの段階で行えばいいでしょうか.
rの最大値の位置を求める時なのか,x1の自己相関関数を求めたあとなのか

r(自己相関関数)の用途が、最大の位置が分かればいいだけ(それ以外には使わない)なら、どちらでもいいです
x1からrを計算した後、rの最大の位置を求めるところまでのどこかでやれば、どこでやっても結果同じですよね
・x1からrを計算した直後に、r(1:511)は捨ててしまう (タイムラグが負の範囲を捨てて、正の範囲だけ残す)
・max()に入れる際に、rの範囲指定に512を足す
のどちらでも、自分が後日コード見たときに、意図が読み取り易いと思う方で (明日の自分は他人ですから、なんでそうしたのかなんて覚えてません)

> r(80)からr(200)の最大値を求めたければ,r(592)とr(712)にしたら良いということでしょうか.

そうです
ただし、512, 80, 200といった数値は、今後変わるかもしれないので、コード中に決め打ちで数字を書くよりも、変数(lenとか)から計算するようにした方がいいですよ

[kaeruuuun]

わかりました.ありがとうございます.