質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
MATLAB

MATLABはMathWorksで開発された数値計算や数値の視覚化のための高水準の対話型プログラミング環境です。

MacOS(OSX)

MacOSとは、Appleの開発していたGUI(グラフィカルユーザーインターフェース)を採用したオペレーションシステム(OS)です。Macintoshと共に、市場に出てGUIの普及に大きく貢献しました。

Q&A

解決済

1回答

2177閲覧

自己相関関数のxcorr

kaeruuuun

総合スコア19

MATLAB

MATLABはMathWorksで開発された数値計算や数値の視覚化のための高水準の対話型プログラミング環境です。

MacOS(OSX)

MacOSとは、Appleの開発していたGUI(グラフィカルユーザーインターフェース)を採用したオペレーションシステム(OS)です。Macintoshと共に、市場に出てGUIの普及に大きく貢献しました。

0グッド

0クリップ

投稿2021/06/02 07:05

編集2021/06/02 07:33

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

scilab

1 2 3//--------------------------------------------------------- 4// 分析用パラメータ設定 5//--------------------------------------------------------- 6len=512; 7order=200; 8sfreq=16; 9fshift=10; 10fwidth=len/sfreq; 11 12 13[x,y]=loadwave('/Users/1.wav'); 14n_sample=y(8); 15disp(n_sample); // 音声データのサンプル数をコンソールに表示 16for i=1:n_sample 17 ts(i)=i/sfreq; 18end 19 20 21//--------------------------------------------------------- 22// Hamming窓を生成 23//--------------------------------------------------------- 24for i=1:len 25 win(i)=0.54-0.46*cos(2*%pi*i/len); 26end 27n_frame=0; 28 29//--------------------------------------------------------- 30// 短時間対数パワーを計算 31//--------------------------------------------------------- 32for st=1:sfreq*fshift:n_sample-len 33 pwr=0; 34 for i=1:len 35 x1(i)=x(i+st)*win(i); 36 pwr=pwr+x1(i)*x1(i); 37 end 38 n_frame=n_frame+1; 39 tp(n_frame)=fwidth/2+(n_frame-1)*fshift; 40 logpwr(n_frame)=10*log10(pwr); 41end 42 43//--------------------------------------------------------- 44// 短時間対数パワーを描く 45//--------------------------------------------------------- 46subplot(3,1,2); 47xlabel("Time [ms]"); 48ylabel("Speech Power [dB]"); 49plot2d(tp,logpwr); 50 51 52//--------------------------------------------------------- 53// 基本周波数の推定(自己相関法) 54//--------------------------------------------------------- 55thlogpwr=-25; 56for n=1:n_frame // n_frameは分析フレーム数 57 if logpwr(n) < thlogpwr then 58 f0(n)=0; 59 else 60 st=1+(n-1)*sfreq*fshift; // 現フレームの開始サンプル番号 61 j=1; 62 for l=0:order 63 64 [r,len]=xcorr(x1,"biased"); //自己相関関数 65 66 67 // 68 end 69 //max=0; 70 //for j=80:200 71 72 argmax=max(r(80:200)); 73 74 period = argmax/16; 75 f0(n) = 1000/period; 76 end 77end 78subplot(3,1,3); 79xlabel("Time [ms]"); 80ylabel("Fundamental Frequency [Hz]"); 81plot2d(tp,f0); 82

またこのxcorrで求めた自己相関関数でこのようなグラフが出力されたのですが私が求めたいのはこのような台形ではなく,ピークが局所的に目立つようにしたいです.
イメージ説明
ピークが目立つようにとはこの図の右の赤いグラフのような形です.
![イメージ説明説明]

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

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

jbpb0

2021/06/02 07:55 編集

> 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

2021/06/02 10:02 編集

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

2021/06/02 10:20

> 音声信号を読み込んで分析窓関数(同じものを繰り返し使うので先に計算)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

2021/06/02 10:26

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

2021/06/02 10:40

上にも書きましたが、 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

2021/06/02 10:49 編集

> 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

2021/06/03 04:20

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

2021/06/04 01:58 編集

少し勉強して考えたプログラムを書きます. 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

2021/06/04 01:51

>本当に、xの最後の512個にしか関係してないx1の自己相関関数が計算できたら、それでいいのですか? xの他の部分は、要らないのですか? 今回分析したいのは512点なので,それで大丈夫です.
kaeruuuun

2021/06/04 01:53

>いろんな「j」に対する値が全部計算されたのが自己相関関数で、xcorr()はそれを1回で計算してくれる と言うことはfor文を回す必要はないのですか?
jbpb0

2021/06/04 02:08 編集

> for文を回す必要はないのですか? 何度も言ってますけど、 https://help.scilab.org/docs/6.1.0/ja_JP/xcorr.html の「例」の計算をそのまま自分でやって、元データ「y」の形(データのサイズ)や、計算結果「c」と「ind」の形を確認してください そうしたら分かります 上記「例」では、forループ無しでxcorr()の計算1回だけして、計算結果の「c」と「ind」でplot()してますよね plot()結果のグラフもありますけど、いろんなタイムラグ(グラフ横軸)に対する自己相関関数の値が、ちゃんとグラフにプロットされてますよね
kaeruuuun

2021/06/04 02:09

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

2021/06/07 04:02 編集

> 今回分析したいのは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()とかのところも直さないといけないけど)
guest

回答1

0

ベストアンサー

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

は、

scilab

1 st=1+(n-1)*sfreq*fshift; // 現フレームの開始サンプル番号 2 j=1; 3 for l=0:order 4 5 [r,len]=xcorr(x1,"biased"); //自己相関関数 6 7 8 // 9 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自体を計算して求めるところのやり方にありますので、そこを直す必要があります

投稿2021/06/07 04:31

jbpb0

総合スコア7651

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

kaeruuuun

2021/06/07 07:56

ありがとうございます. >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

2021/06/07 08:08 編集

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

2021/06/07 08:31 編集

> タイムラグとして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

2021/06/07 09:20

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問