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

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

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

Используя функцию rica, можно выполнить это восстановление даже, не зная смесительной матрицы или среднего значения. Учитывая набор нескольких наблюдений..., rica извлекает исходные сигналы....

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

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

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

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

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

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

Предположим, что это - матрица левых собственных векторов положительной полуопределенной матрицы и является матрицей собственных значений. Затем

Затем

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

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

предварительно побеленные данные. rica вычисляет неизвестную матрицу под предположением, что компоненты максимально независимы.

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

Документация