Загрузите и постройте мультисигнал

Чтобы проиллюстрировать эти возможности, давайте рассмотрим реальный мультисигнал, представляющий 35 дней потребления электрической нагрузки, сосредоточенного и стандартизированного.

load elec35_nor;       % Load normalized original data.
X = signals;
[nbSIG,nbVAL] = size(X);
plot(X','r')
axis tight
title('Original Data: 35 days of electrical consumption');
xlabel('Minutes from 1 to 1440');

Мы видим, что сигналы являются локально неправильными и шумными, но тем не менее три различных общих формы можно отличить.

Поверхностное представление данных

Для того, чтобы подсветить периодичность мультисигнала, давайте теперь исследуем данные с помощью 3-D представления.

surf(X)
shading interp
axis tight
title('Original Data: 35 days of electrical consumption')
xlabel('Minutes from 1 to 1440','Rotation',4)
ylabel('Days from 1 to 35','Rotation',-60)
ax = gca;
ax.View = [-13.5 48];

Эти пять представленных недель могут теперь быть замечены более ясно.

Разложение строки мультисигнала

Выполните разложение вейвлета на уровне 7 с помощью 'sym4' вейвлета.

dirDec = 'r';         % Direction of decomposition
level  = 7;           % Level of decomposition
wname  = 'sym4';      % Near symmetric wavelet
decROW = mdwtdec(dirDec,X,level,wname);

Сгенерированная структура разложения следующие:

decROW

decROW = 

  struct with fields:

        dirDec: 'r'
         level: 7
         wname: 'sym4'
    dwtFilters: [1x1 struct]
       dwtEXTM: 'sym'
      dwtShift: 0
      dataSize: [35 1440]
            ca: [35x18 double]
            cd: {1x7 cell}

Сигналы и приближения на уровне 7

Давайте восстановим приближения на уровне 7 для каждого сигнала строки. Затем чтобы сравнить приближения с исходными сигналами, давайте отобразим два графика (см. рисунок ниже). Первый показывает все исходные сигналы, и второй показывает все соответствующие приближения на уровне 7.

A7_ROW  = mdwtrec(decROW,'a',7);

subplot(2,1,1)
plot(X(:,:)','r')
title('Original Data')
axis tight
subplot(2,1,2)
plot(A7_ROW(:,:)','b')
axis tight
title('Corresponding approximations at level 7')

Как это видно, общая форма получена приближениями на уровне 7, но некоторые интересные функции потеряны. Например, удары вначале и в конце сигналов исчезли.

Наложенные сигналы и приближения

Чтобы сравнить более тесно исходные сигналы с их соответствующими приближениями на уровне 7, давайте отобразим два графика (см. рисунок ниже). Первый касается четырех первых сигналов и второго пять последних сигналов. Каждый из этих графиков представляет исходные сигналы, наложенные их соответствующим приближением на уровне 7.

subplot(2,1,1)
idxDays = 1:4;
plot(X(idxDays,:)','r')
hold on
plot(A7_ROW(idxDays,:)','b')
axis tight
title(['Days ' int2str(idxDays), ' - Signals and Approximations at level 7'])
subplot(2,1,2)
idxDays = 31:35;
plot(X(idxDays,:)','r')
hold on
plot(A7_ROW(idxDays,:)','b')
axis tight
title(['Days ' int2str(idxDays), ' - Signals and Approximations at level 7'])

Как это видно, приближение исходного сигнала визуально точно в терминах общей формы.

Шумоподавление мультисигнала

Чтобы выполнить более тонкое упрощение мультисигнала, сохраняющего эти удары, которые явно относятся к электрическому сигналу, позвольте нам denoise мультисигнал.

Процедура шумоподавления выполняется с помощью трех шагов:

1) Разложение: Выберите вейвлет и уровень разложения N, и затем вычислите разложения вейвлета сигналов на уровне N.

2) Пороговая обработка: Для каждого уровня от 1 до N и для каждого сигнала, выбран порог, и пороговая обработка применяется к коэффициентам детали.

3) Реконструкция: Вычислите реконструкции вейвлета с помощью исходных коэффициентов приближения уровня N и модифицированных коэффициентов детали уровней от 1 до N.

Давайте выберем теперь уровень разложения N = 5 вместо N = 7 используемых ранее.

dirDec = 'r';         % Direction of decomposition
level  = 5;           % Level of decomposition
wname  = 'sym4';      % Near symmetric wavelet
decROW = mdwtdec(dirDec,X,level,wname);

[XD,decDEN] = mswden('den',decROW,'sqtwolog','mln');
Residuals = X-XD;

subplot(3,1,1)
plot(X','r')
axis tight
title('Original Data: 35 days of electrical consumption')
subplot(3,1,2)
plot(XD','b')
axis tight
title('Denoised Data: 35 days of electrical consumption')
subplot(3,1,3)
plot(Residuals','k')
axis tight
title('Residuals')
xlabel('Minutes from 1 to 1440')

Качество результатов хорошо. Удары вначале и в конце сигналов хорошо восстанавливаются и, с другой стороны, остаточные значения похожи на шум за исключением некоторых остающихся ударов из-за сигналов. Кроме того, величина этих остающихся ударов имеет маленький порядок.

Мультисигнал, сжимающий сигналы строки

Как шумоподавление, процедура сжатия выполняется с помощью трех шагов (см. выше).

Различие для процедуры шумоподавления найдено на шаге 2. Существует два доступные подхода сжатия:

Чтобы упростить представление данных и сделать более эффективным сжатие, давайте вернемся к разложению на уровне 7 и давайте сожмем каждую строку мультисигнала путем применения глобального порогового продвижения, чтобы восстановить 99% энергии.

dirDec = 'r';         % Direction of decomposition
level  = 7;           % Level of decomposition
wname  = 'sym4';      % Near symmetric wavelet
decROW = mdwtdec(dirDec,X,level,wname);

[XC,decCMP,THRESH] = mswcmp('cmp',decROW,'L2_perf',99);

subplot(3,1,1)
plot(X','r')
axis tight
title('Original Data: 35 days of electrical consumption')
subplot(3,1,2)
plot(XC','b')
axis tight
title('Compressed Data: 35 days of electrical consumption')
subplot(3,1,3)
plot((X-XC)','k')
axis tight
title('Residuals')
xlabel('Minutes from 1 to 1440')

Как это видно, общая форма сохраняется, в то время как локальными неисправностями пропускают. Остаточные значения содержат шум, а также компоненты из-за мелких масштабов по существу.

Эффективность сжатия

Давайте теперь вычислим соответствующую плотность ненулевых элементов.

cfs = cat(2,[decCMP.cd{:},decCMP.ca]);
cfs = sparse(cfs);
perf = zeros(1,nbSIG);
for k = 1:nbSIG
    perf(k) = 100*nnz(cfs(k,:))/nbVAL;
end

figure
plot(perf,'r-*')
title('Percentages of nonzero coefficients for the 35 days')
xlabel('Signal indices')
ylabel('% of nonzero coefficients')

Для каждого сигнала процент необходимых коэффициентов, чтобы восстановить 99% энергии находится между 1,25% и 1,75%. Это иллюстрирует возможность вейвлетов сконцентрировать энергию сигнала в немногих коэффициентах.

Кластеризация сигналов строки

Кластеризация предлагает удобную процедуру, чтобы обобщить большой набор сигналов с помощью разреженных представлений вейвлета. Можно реализовать иерархическую кластеризацию с помощью mdwtcluster. У вас должен быть Statistics and Machine Learning Toolbox™, чтобы использовать mdwtcluster.

Сравните три различного объединения в кластеры 35-дневных данных. Первый основан на исходном мультисигнале, втором на коэффициентах приближения на уровне 7, и последний основан на мультисигнале denoised.

Установите на количество кластеров к 3. Вычислите структуры P1 и P2, которые содержат соответственно два первых раздела и третий.

P1 = mdwtcluster(decROW,'lst2clu',{'s','ca7'},'maxclust',3);
P2 = mdwtcluster(decDEN,'lst2clu',{'s'},'maxclust',3);
Clusters = [P1.IdxCLU P2.IdxCLU];

Мы можем теперь протестировать равенство этих трех разделов.

EqualPART = isequal(max(diff(Clusters,[],2)),[0 0])

EqualPART =

  logical

   1

Таким образом, эти три раздела являются тем же самым. Давайте теперь построим и исследовать кластеры.

figure
stem(Clusters,'filled','b:')
title('The three clusters of the original 35 days')
xlabel('Signal indices')
ylabel('Index of cluster')
ax = gca;
xlim([1 35])
ylim([0.5 3.1])
ax.YTick = 1:3;

Первый кластер (пометил 3) содержит 25 дней в середине недели, и оба другие (пометил 2, и 1) содержат пять суббот и пять воскресений соответственно. Это иллюстрирует снова периодичность базовых временных рядов и трех различных общих форм, замеченных в двух первых графиках, отображенных в начале этого примера.

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

CA7 = mdwtrec(decROW,'ca');
IdxInCluster = cell(1,3);
for k = 1:3
    IdxInCluster{k} = find(P2.IdxCLU==k);
end
figure('Units','normalized','Position',[0.2 0.2 0.6 0.6])
for k = 1:3
    idxK = IdxInCluster{k};
    subplot(2,3,k)
    plot(X(idxK,:)','r')
    axis tight
    ax = gca;
    ax.XTick = [200 800 1400];
    if k==2
        title('Original signals')
    end
    xlabel(['Cluster: ' int2str(k) ' (' int2str(length(idxK)) ')'])
    subplot(2,3,k+3)
    plot(CA7(idxK,:)','b')
    axis tight
    if k==2
        title('Coefficients of approximations at level 7')
    end
    xlabel(['Cluster: ' int2str(k) ' (' int2str(length(idxK)) ')'])
end

Те же разделы получены из исходных сигналов (1 440 выборок для каждого сигнала) и от коэффициентов приближений на уровне 7 (18 выборок для каждого сигнала). Это иллюстрирует, что использования меньше чем 2% коэффициентов достаточно, чтобы получить те же разделы кластеризации этих 35 дней.

Сводные данные

Шумоподавление, сжатие и кластеризирующий использование вейвлетов являются очень эффективными инструментами. Возможность представлений вейвлета сконцентрировать энергию сигнала в немногих коэффициентах является ключом КПД. Кроме того, кластеризация предлагает удобную процедуру, чтобы обобщить большой набор сигналов с помощью разреженных представлений вейвлета.