Q&A
もしかしたら、「*」と「.*」を混同してませんか？
https://jp.mathworks.com/help/matlab/ref/mtimes.html
https://jp.mathworks.com/help/matlab/ref/times.html
数学での行列の掛け算は、「*」の方です
非常に基礎的な質問です.
プログラムの作成過程で行列の計算がわからなくなったため簡単な問題に置き換えています.
matlab
1A = 2 3 1 2 3 4 4 5 6 5B = 6 7 1 2 3 8 4 5 6 9C=A.' 10 11C = 12 13 1 4 14 2 5 15 3 6 16>> B.*C 17配列のサイズがこの演算に適合しません。 18
のようにエラーが出ます.2*3の行列の片方を転置したものをかける時はどのようにしたら良いでしょうか.
もしかしたら別の質問に移るかもしれませんが,私は今このφ(y)^T*yの式で上と同じ配列のサイズが合わないというエラーでつまずいています.
φ,yともに1000×2doubleでその部分だけを取り出すとこのようなプログラムを書いていました.
それでエラーが出るので,式を簡単にして考えていました.このプログラムは間違っていますか.
matlab
1a=fai.'; 2b=y; 3z=a.*b;
matlab
1[x fs]=audioread('.wav'); 2x_length=length(x); 3y_out=zeros(x_length,2); 4w=[1 0;0 1]; 5ita=0.0001; 6Ta=1000; 7 8for n=1:x_length-Ta 9 10 x_tmp=x(n:n+Ta-1,:); 11 y=x_tmp*w.'; 12 fai=tanh(y); 13 I=eye(n); 14 y_out(n,:)=y(1,:); 15 16 w=w+ita*(I-((fai.'.*y)/Ta)).*w; 17 18end
使っているプログラムと今のデータセットを載せます.
お久しぶりです.
コメントいただきありがとうございます.
質問文を修正しました.
よろしければ確認をお願いいたします.
やりたいのは行列ではなく,φ(y)^T*yです.
であれば.*を使うのが妥当ではないでしょうか.
その場合はφ(y)^T*yのy(:,1)のように列を指定する必要があるのでしょうか.
同時に2列の計算はできないのでしょうか.
> φ,yともに1000×2doubleで
を3x2として、
fai=
a11 a12
a21 a22
a31 a32
y=
b11 b12
b21 b22
b31 b32
の場合、どのような計算をしたいのでしょうか？
ちなみに、fai .* y は
a11*b11 a12*b12
a21*b21 a22*b22
a31*b31 a32*b32
となります
また、fai’ は
a11 a21 a31
a12 a22 a32
なので、fai’ * y は
c11=a11*b11+a21*b21+a31*b31
c12=a11*b12+a21*b22+a31*b32
c21=a12*b11+a22*b21+a32*b31
c22=a12*b12+a22*b22+a32*b32
として、
c11 c12
c21 c22
となります
fai'にbをかけてa11*b11のように要素ごとに掛け算をしたいです.
念の為プログラムを載せました.最後のwの式でエラーが出ています.サイズが合わないという質問と同様のエラーが出ています.
配列のサイズがこの演算に適合しません。
エラー: (行 19)
w=w+ita*(I-((fai.'.*y)/Ta)).*w;
> fai'にbをかけてa11*b11のように要素ごとに掛け算をしたいです.
それでしたら、欲しい計算結果は、
a11*b11 a12*b12
a21*b21 a22*b22
a31*b31 a32*b32
ではないのですか？
違うなら、欲しい計算結果を、a11とかを使って書いてください
【追記】 fai.*y に変えたら、式の一番右の .*w の計算ができないな
そこを *w に変えても、式の先頭の w+ の計算ができない
やはり、fai.'*y ではないのかなぁ
(ozwkさんも回答で書いてるように)
あと、その数式について書かれてるWebページとかありませんか？
まず、私の最初のコメントに書いたURLの説明を読んでください
それを読むと分かるように、
A=[1 2 3; 4 5 6]
B=A
C=A'
の場合
A.*B
B*C
C*B
は計算できますが、
B.*C
C.*B
は計算できません
なので、
> 配列のサイズがこの演算に適合しません。
w=w+ita*(I-((fai.'.*y)/Ta)).*w;
の
fai.'.*y
は計算できません
https://www.kecl.ntt.co.jp/icl/signal/sawada/mypaper/subspace2010rev.pdf
ここに載っていることを実現したいです.
wの式はp.38の最後の式です.
38ページに「以下を収束まで繰り返す」って書いてありますから、一回で終わりではないのだと思いますよ
下記はICAのPythonコードですが、計算が収束するまで何回も繰り返してますよね
https://aidemy.net/magazine/685/#i-3
(コードの「#収束判定」のところ)
こんな感じで繰り返し計算するのではないですかね
ICAにも、いろんな計算方法があるみたいなので、上記Pythonコードは質問者さんがやりたいこととアルゴリズムは微妙に違うのかもしれませんが
for分で回してその都度更新しているプログラムのつもりですが,そうなっていないということでしょうか.
今回のfor分では信号長分だけ学習させるようにしています.
> 信号長分だけ学習させる
それだと、信号の長さで学習の回数が決まるので、十分な回数の学習がされないかもしれませんよね
https://aidemy.net/magazine/685/#i-3
のコードの「analyze」という関数内の処理を見ると、音声データ「z」を分けずに全部を一度に使って「vec_w」を計算してます
計算した「vec_w」と、一回前の「vec_w」である「vec_w_prev」との差がわずかになるまで、「z」から「vec_w」を計算するのを繰り返してます (while True: のループ)
一回前の「vec_w」である「vec_w_prev」との差がわずかになるまで繰り返し計算する、すなわち、計算を繰り返しても「vec_w」がほとんど変わらなくなる状態まで計算を繰り返すのが、
「収束まで繰り返す」
ということです
信号を分けて処理すると、信号全体が同一の素性なのか？ という問題もありますよね
信号の最初の方と最後の方で傾向に何か違いがあると、部分に分けて学習すると、最終的に学習されて決まったパラメータが、信号のあるところには最適なんだけど、他のところには適してない、となってしまうかも
信号の全体にまんべんなく適しているように学習するには、信号を分けずに全部を一度に使って学習した方がいいのではないですかね
もちろん、信号のサイズがでかすぎて、分けないと処理できない、という場合もあるとは思いますが
ozwkさんもおっしゃっているように信号を分割しないでwを更新するとよいとのことですが,分割をしないとは具体的にどういうことでしょうか.
今回問題の設定条件として最初のほうは分離できていなくてもよいということだったので分割していました.それ以前にやはり後半のほうでも少しもう片方の音が残っているようなのですが,それはどのようにしたら消えるのでしょうか.現状として音の後半のほうは出力したいほうの音が非常に大きく聞こえて一見分離できているように聞こえるのですが,音の隙間では小さく聞こえてほしくないもう片方の音が聞こえている状況です.
> 分割をしないとは具体的に
私はPDF文献をちゃんと理解しているわけではないのでよく分かりませんが、元データxと、文献の式のyの関係が、
y=x_tmp*w.'
であるというのが正しいのなら、「x_tmp」を使わず「x」を使うということになるのだと思います
なので、「for n=1:x_length-Ta」のループは不要になります
その代わりに、wが収束するまで繰り返すループを入れます
ループの最初で、たとえば
w_old = w
みたいにして更新前の「w」を保管しておいて、
w=w+...
で「w」を更新してから、
w - w_old
を計算して、それが十分に小さくなったらループを抜ける、みたいな
実際は「w」は単なる数値ではなく配列なので、
w - w_old
を、たとえば
norm(w - w_old)
のようにして、単一の数値に変えてから収束判定することになります
> もう片方の音が残っている
> どのようにしたら消える
その理由によりますよね
・PDF文献のアルゴの限界
・質問者さんが、アルゴを正しくコードに書けてない
・学習不足 (wが収束するまで学習できてない)
が考えられますけど、そのどれなのかの分析を他人に任せるのは、おかしいです
それをやりたいのは質問者さんなのだから、自分で調べてください
ただし、文献には「収束まで繰り返す」と書いてありますから、その通りに「w」が収束するまで繰り返すコードを書いて、それでも消えないかを確認しないと、アルゴの限界かどうかは分かりませんよね
「w」が収束するまで繰り返すコードを書いても消えない場合は、他の同様なアルゴのプログラムと性能を比較してみると、いいかもしれません
たとえば、
https://aidemy.net/magazine/685/
の「mix_1.wav」〜「mix_3.wav」をダウンロードして、それを質問者さんのコードで処理した結果と、上記Webページの処理結果「music1.wav」〜「music3.wav」を比べてみる
(入力データが三つになるので、その変更はちゃんとやるとして)
その結果が同等なら、おそらく実力がそのくらい (アルゴの限界) なのでしょうね
質問者さんがテストしているデータが、分離が難しいものなのだろうと思います
もし、質問者さんのコードの結果の方が劣るなら、たぶん質問者さんのコードが正しくない (文献通りに計算できてない) のだと思います
あるいは、もし、質問者さんがPythonのコードを動かせる環境をお持ちなら、
https://aidemy.net/magazine/685/
のPythonコードで、質問者さんが現状テストしているデータを「mix_1.wav」〜「mix_3.wav」の代わりに使って分離してみると、いいと思いますよ
その場合でも、
> もう片方の音が残っている
のなら、ICAによる分離の性能は、そんなものなのでしょうね
わかりました.プログラムの問題の気がしているので,もう少し考えてみます.ありがとうございます.
回答1件
0
ベストアンサー
やりたいのは行列ではなく,φ(y)^T*yです.
であれば.*を使うのが妥当ではないでしょうか.
いや、行列の積(*)だと思います
wが2x2行列で
w = w + η(I-φ(y)^t y )w
ということはη(I-φ(y)^t y )は2x2でなければなりません。
φ(y), yが共に1000x2 ということは
φ(y)^t yは2x1000行列と1000x2行列の積なので2x2行列になります
(同様の理由でIもI=eye(2)かと思います)
投稿2021/07/06 23:51
編集2021/07/06 23:59
総合スコア13553
w=w+ita*(I-((fai.'.*y)/Ta)).*w;
の一番右の
.*w
も、実は
*w
が正解だったりして
ozwkさん
質問者さんが画像で上げてる数式で、そこだけ「*」が使われてるのですが、それって何か意味があるのでしたっけ？
ご回答ありがとうございます.
w=w+ita*(I-((fai.'*y)/Ta))*w;
にするとプログラム自体はじっこうできるのですが､
sound(y_out(:,1),fs)
として片方の音源を聞くと最初のほうが音源の分離ができておらず混合した声が聞こえて,後半は音源の分離ができていて片方の音源が聞こえているという形になります.
やりたいことは二つの混合された音声を分離させるというプログラムを作成するということです.
w=w+ita*(I-((fai.'*y)/Ta)).*w;でするとただxが再生されるだけだったので,
w=w+ita*(I-((fai.'*y)/Ta))*w;
のほうが正解に近いと思います.ただ,音源の最初のがうまく分離されていないという状況です.
上はI=eye(2)にしてやってみた結果です.
このアルゴリズムでWを更新して、
最終的に求まったWに対して改めてy_out=Wxとしてy_outを求めてください
更新途中のWに対して都度y_outを求めてしまっているので最初がうまく分離できていないという結果になっています
[x fs]=audioread('.wav');
x_length=length(x);
y_out=zeros(x_length,2);
w=[1 0;0 1];
ita=0.0001;
Ta=1000;
for n=1:x_length-Ta
x_tmp=x(n:n+Ta-1,:);
y=x_tmp*w.';
fai=tanh(y);
I=eye(n);
y_out(n,:)=y(1,:);
w=w+ita*(I-((fai.'*y)/Ta))*w;
end
sound(y_out(:,1),fs)
のようにしてみたのですが,これだとうまく更新できていないということでしょうか.
endの後ろに
y_out=x*w.';
を追加してみると最初は分離できているのですが少し後ろでもう片方の音が聞こえます.
1000個ずつに分割してますが、分割せずに全体に対してw=w+…の更新式を繰り返すものではないでしょうか
繰り返し回数はwが収束するまでですが、まあとりあえず適当に大きな数でいいんじゃないでしょうか
問題の条件として最初はもう片方の音が入っていても良いという条件だったのを忘れていました.
あまり精度が良くはありませんが,後半は分割できているので,プログラムとしては問題なさそうです.ありがとうございました.
あなたの回答
tips
プレビュー
