Извлеките смешанные сигналы

В этом примере показано, как использовать rica распутывать смешанные звуковые сигналы. Можно использовать rica выполнять независимый анализ компонента (ICA), когда предварительное отбеливание включено как шаг предварительной обработки. Модель ICA

$$x = \mu + As.$$

Здесь,$x$$p$-by-1 вектор из смешанных сигналов,$\mu$$p$-by-1 вектор из значений смещения,$A$ $p$$q$смесительная матрица и$s$$q$-by-1 вектор из исходных сигналов. Предположим сначала, что это$A$ - квадратная матрица. Если вы знаете$\mu$ и$A$, можно восстановить исходный сигнал$s$ с данных$x$:

$$s = A^{-1}(x - \mu).$$

Используя rica функция, можно выполнить это восстановление даже, не зная смесительной матрицы$A$ или среднего значения$\mu$. Учитывая набор нескольких наблюдений$x(1)$$x(2)$..., rica извлекает исходные сигналы$s(1)$$s(2)$....

Загрузка данных

Загрузите набор шести звуковых файлов, которые поставляются с MATLAB®. Обрежьте каждый файл к 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 так, чтобы это имело нулевое среднее значение и единичную ковариацию.

Идея следующая. Если$s$ нулевой средний источник со статистически независимыми компонентами, то

$$E(s) = 0$$

$$E(ss^T) = I.$$

Затем среднее значение и ковариация$x$

$$ E(x) = \mu$$

$${\rm Cov}(x) = AA^T = C.$$

Предположим, что вы знаете$\mu$ и$C$. На практике вы оценили бы эти количества от демонстрационного среднего значения и ковариации столбцов$x$. Можно решить для$s$ в терминах$x$

$$ s = A^{-1}(x-\mu) = (A^TA)^{-1}A^T(x-\mu).$$

Последнее уравнение содержит, даже когда$A$ не квадратная обратимая матрица.

Предположим, что это$U$ - $p$$q$матрица левых собственных векторов положительной полуопределенной матрицы $C$и$\Sigma$ является $q$$q$матрицей собственных значений. То

$$ C = U\Sigma U^T$$

$$U^TU = I.$$

То

$$ AA^T = U\Sigma U^T.$$

Существует много смесительных матриц$A$, которые удовлетворяют этому последнему уравнению. Если$W$ $q$$q$ортонормированной матрицей, то

$$ W^TW = W W^T = I$$

$$ A = U\Sigma^{1/2} W.$$

Замена в уравнение для$s$,

$$s = W^T\tilde{x},\ \rm{where}$$

$$ \tilde{x} = \Sigma^{-1/2}U^T(x-\mu).$$

$\tilde{x}$ предварительно побеленные данные. rica вычисляет неизвестную матрицу$W$ под предположением, что компоненты$s$ максимально независимы.

mixdata = prewhiten(mixdata);

Разделите все сигналы

Супергауссов источник имеет резкий пик около нуля, такого как гистограмма звукового 1 показывает.

figure
histogram(S(:,1))

Выполните 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

Смотрите также

| | |

Связанные примеры

Больше о

Для просмотра документации необходимо авторизоваться на сайте