ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

混合信号の抽出

この例では、rica を使用して混合オーディオ信号を分離する方法を示します。前置白色化が前処理ステップとして含まれている場合、rica を使用して独立成分分析 (ICA) を行えます。ICA モデルは次のとおりです。

ここで、 行 1 列の混合信号のベクトル、 行 1 列のオフセット値のベクトル、 列の混同行列、 行 1 列の元の信号のベクトルです。まず、 が正方行列であるとします。 がわかっている場合は、元の信号 をデータ から復元できます。

関数 rica を使用すると、混同行列 または平均 が不明な場合でも、この復元を実行できます。rica は、与えられた観測値 、... から元の信号 、... を抽出します。

データの読み込み

MATLAB® に含まれている 6 つのオーディオ ファイルを読み込みます。各ファイルをトリミングして 10,000 個の標本にします。

files = {'chirp.mat'
        'gong.mat'
        'handel.mat'
        'laughter.mat'
        'splat.mat'
        'train.mat'};

S = zeros(10000,6);
for i = 1:6
    test     = load(files{i});
    y        = test.y(1:10000,1);
    S(:,i)   = y;
end

信号の混合

ランダム混同行列を使用して信号を混合し、ランダム オフセットを追加します。

rng default % For reproducibility
mixdata = S*randn(6) + randn(1,6);

元の音声を聞くため、次のコードを実行します。

   for i = 1:6
       disp(i);
       sound(S(:,i));
       pause;
   end

混合した音声を聞くため、次のコードを実行します。

   for i = 1:6
       disp(i);
       sound(mixdata(:,i));
       pause;
   end

信号をプロットします。

figure
for i = 1:6
    subplot(2,6,i)
    plot(S(:,i))
    title(['Sound ',num2str(i)])
    subplot(2,6,i+6)
    plot(mixdata(:,i))
    title(['Mix ',num2str(i)])
end

元の信号には明確な構造があります。混合信号にはほとんど構造がありません。

混合信号の前置白色化

信号を効果的に分離するため、この例の終わりにある関数 prewhiten を使用して信号を "事前に白色化" します。この関数は、ゼロ平均と単位共分散をもつように mixdata を変換します。

考え方は以下のとおりです。 がゼロ平均ソースであり、その成分が統計的に独立している場合、次のようになります。

すると、 の平均と共分散は次のようになります。

は既知であるとします。実際には、標本の平均と の列の共分散からこれらの数量を推定します。次のように、 に関して を解くことができます。

が正方可逆行列ではない場合、後者の方程式は一定になります。

が半正定値行列 の左固有ベクトルによる 列の行列であり、 列の固有値の行列であるとします。その場合

また、以下とします。

この最後の方程式を満たす混同行列 は複数あります。 列の正規直交行列である場合、次のようになります。

についての方程式に代入します。

は事前に白色化されたデータです。rica は、 の成分が可能な限り独立であると仮定して、未知の行列 を計算します。

mixdata = prewhiten(mixdata);

すべての信号の分離

音声 1 のヒストグラムに示されているように、優ガウス性のソースではゼロの付近に鋭いピークがあります。

figure
histogram(S(:,1))

6 つの特徴量を要求して再構成 ICA を実行します。各ソースが優ガウス性であることを指定します。

q = 6;
Mdl = rica(mixdata,q,'NonGaussianityIndicator',ones(6,1));

特徴量を抽出します。分離手順に成功した場合、特徴量は元の信号に比例します。

unmixed = transform(Mdl,mixdata);

分離信号と元の信号の比較

元の信号と分離信号をプロットします。

figure
for i = 1:6
    subplot(2,6,i)
    plot(S(:,i))
    title(['Sound ',num2str(i)])
    subplot(2,6,i+6)
    plot(unmixed(:,i))
    title(['Unmix ',num2str(i)])
end

分離信号の順序は元の順序と異なります。列を並べ替えて、対応する元の信号と分離信号を一致させます。分離信号をスケーリングして、対応する元の信号と同じノルムにします (どのスケールでも同じ信号混合になるので、rica は元の信号のスケールを識別できません)。

unmixed = unmixed(:,[2,5,4,6,3,1]);
for i = 1:6
    unmixed(:,i) = unmixed(:,i)/norm(unmixed(:,i))*norm(S(:,i));
end

元の信号と分離信号をプロットします。

figure
for i = 1:6
    subplot(2,6,i)
    plot(S(:,i))
    ylim([-1,1])
    title(['Sound ',num2str(i)])
    subplot(2,6,i+6)
    plot(unmixed(:,i))
    ylim([-1,1])
    title(['Unmix ',num2str(i)])
end

分離信号は元の信号と同じように見えます。分離した音声を聞くため、次のコードを実行します。

   for i = 1:6
       disp(i);
       sound(unmixed(:,i));
       pause;
   end

以下は、関数 prewhiten のコードです。

function Z = prewhiten(X)
% X = N-by-P matrix for N observations and P predictors
% Z = N-by-P prewhitened matrix

    % 1. Size of X.
    [N,P] = size(X);
    assert(N >= P);

    % 2. SVD of covariance of X. We could also use svd(X) to proceed but N
    % can be large and so we sacrifice some accuracy for speed.
    [U,Sig] = svd(cov(X));
    Sig     = diag(Sig);
    Sig     = Sig(:)';

    % 3. Figure out which values of Sig are non-zero.
    tol = eps(class(X));
    idx = (Sig > max(Sig)*tol);
    assert(~all(idx == 0));

    % 4. Get the non-zero elements of Sig and corresponding columns of U.
    Sig = Sig(idx);
    U   = U(:,idx);

    % 5. Compute prewhitened data.
    mu = mean(X,1);
    Z = bsxfun(@minus,X,mu);
    Z = bsxfun(@times,Z*U,1./sqrt(Sig));
end

参考

| | |

関連する例

詳細