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

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