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

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

Исходные сигналы имеют четкую структуру. Смешанные сигналы имеют гораздо меньшую структуру.

Предварительные смешанные сигналы

Чтобы эффективно разделять сигналы, «предварительно whiten» сигналы при помощи prewhiten функция, которая появляется в конце этого примера. Эта функция преобразует mixdata так, что это имеет нулевое среднее и тождества ковариации.

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

Тогда среднее и ковариация

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

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

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

Тогда

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

Подстановка в уравнение на,

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

mixdata = prewhiten(mixdata);

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

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

figure
histogram(S(:,1))

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