1D мультисигнал является набором 1D сигналов той же длины, сохраненной как матрица, организованная построчный (или по столбцам).
Цель этого примера состоит в том, чтобы показать, как анализировать, denoise или сжать мультисигнал, и затем кластеризировать различные представления или упрощенные версии сигналов, составляющих мультисигнал.
В этом примере сигналы сначала анализируются и различные представления, или производятся упрощенные версии:
Восстановленные приближения на данных уровнях,
Версии Denoised,
Сжатые версии.
Шумоподавление и сжатие являются двумя из главных приложений вейвлетов, часто используемых в качестве шага предварительной обработки перед кластеризацией.
Последний шаг этого примера выполняет несколько кластеризирующихся стратегий, и сравните их. Это позволяет обобщать большой набор сигналов с помощью разреженных представлений вейвлета.
Чтобы проиллюстрировать эти возможности, давайте рассмотрим реальный мультисигнал, представляющий 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.
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('Color','w') 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('Color','w') 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('Color','w','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 дней.
Шумоподавление, сжатие и кластеризирующий использование вейвлетов являются очень эффективными инструментами. Возможность представлений вейвлета сконцентрировать энергию сигнала в немногих коэффициентах является ключом КПД. Кроме того, кластеризация предлагает удобную процедуру, чтобы обобщить большой набор сигналов с помощью разреженных представлений вейвлета.