Загрузить данные

Загрузите набор из шести аудиофайлов, поставляемых с MATLAB ®. Обрезать каждый файл до 10000 образцов.

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

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

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

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

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

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

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

Тогда

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

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

- предварительно выделенные данные. 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