Шумоподавление зависимого интервала вейвлета

Пример показывает как denoise сигнал с помощью зависимых интервалом порогов.

Инструменты Wavelet GUI обеспечивают точный процесс шумоподавления, позволяя нам точно настроить параметры, требуемые к denoise сигнал. Затем мы можем сохранить сигнал denoised, разложение вейвлета и все параметры шумоподавления.

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

В этом примере мы выполняем шесть испытаний к denoise электрический сигнал nelec использующий эти процедуры:

Используя один интервал с минимальным пороговым значением: 4.5
Используя один интервал с максимальным пороговым значением: 19.5
Вручную выбор трех интервалов и трех порогов, и использование wthresh функционируют к порогу коэффициенты.
Используя utthrset_cmd функционируют, чтобы автоматически найти интервалы и пороги.
Используя cmddenoise функционируют, чтобы выполнить все процессы автоматически.
Используя cmddenoise функционируют с дополнительными параметрами.

Шумоподавление Используя один интервал

Загрузите электрический nelec сигнала потребления.

load nelec.mat; 
sig = nelec;
denPAR = {[1 94 5.9 ; 94 1110 19.5 ; 1110 2000 4.5]};

Используя инструменты GUI, чтобы выполнить зависимую интервалом пороговую обработку для сигнала nelec и определение номера интервалов в три, мы получаем содержимое переменной denPAR, которая может быть интерпретирована можно следующим образом:

I1 = [1 94] с порогом thr1 = 5.9
I2 = [94 1110] с порогом thr2 = 19.5
I3 = [1110 2000] с порогом thr3 = 4.5

Теперь мы используем дискретный анализ вейвлета на уровне 5 с sym4 вейвлетом.

wname = 'sym4';
level = 5;
sorh  = 's'; % type of thresholding

Denoise сигнал с помощью функционального wdencmp с пороговым набором в 4,5, который является минимальным значением, обеспеченным инструментами GUI.

thr = 4.5;
[sigden_1,~,~,perf0,perfl2] = wdencmp('gbl',sig,wname,level,thr,sorh,1);
res = sig-sigden_1;
subplot(3,1,1);
plot(sig,'r'); 
axis tight
title('Original Signal');
subplot(3,1,2);
plot(sigden_1,'b'); 
axis tight
title('Denoised Signal');
subplot(3,1,3);
plot(res,'k'); 
axis tight
title('Residual');

perf0,perfl2

perf0 = 66.6995

perfl2 = 99.9756

Полученный результат не хорош. Шумоподавление очень эффективно вначале и конец сигнала, но между 100 и 1100 не удален шум. Обратите внимание на то, что perf0 значение дает процент обнуленных коэффициентов, и perfl2 значение дает процент сохраненной энергии.

Теперь, мы denoise сигнал с максимальным значением, обеспеченным инструментами GUI для порога, 19.5

thr = 19.5;
[sigden_2,cxd,lxd,perf0,perfl2] = wdencmp('gbl',sig,wname,level,thr,sorh,1);
res = sig-sigden_2;
subplot(3,1,1);
plot(sig,'r');   
axis tight
title('Original Signal');
subplot(3,1,2);
plot(sigden_2,'b'); 
axis tight
title('Denoised Signal');
subplot(3,1,3);
plot(res,'k');   
axis tight
title('Residual');

perf0,perfl2

perf0 = 94.7860

perfl2 = 99.9564

Сигнал denoised очень сглажен. Кажется довольно хорошим, но если мы смотрим на невязку после положения 1100, мы видим, что отклонение базового шума не является постоянным. Некоторые компоненты сигнала, вероятно, остались в невязке, такой как, около положения 1300 и между положениями 1900 и 2000.

Шумоподавление Используя Зависимые пороги интервала (IDT)

Теперь мы будем использовать зависимую интервалом пороговую обработку, как в инструментах GUI шумоподавления.

Задайте зависимые интервалом пороги.

thrParams = cell(1,level);
thrParams(:) = denPAR;

Выполните дискретный анализ вейвлета сигнала.

[coefs,longs] = wavedec(sig,level,wname);

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

% Replicate the coefficients
cfs_beg = wrepcoef(coefs,longs);

% Display the coefficients of the decomposition
clf
subplot(6,1,1);
plot(sig,'r');
axis tight
title('Original Signal and Detail Coefficients from 1 to 5');
ylabel('S','Rotation',0);
for k = 1:level
    subplot(6,1,k+1); plot(cfs_beg(k,:),'Color',[0.5 0.8 0.5]);
    ylabel(['D' int2str(k)],'Rotation',0);
    axis tight 
    hold on
    maxi = max(abs(cfs_beg(k,:)));
    hold on
    par = thrParams{k};
    plotPar = {'Color','m','LineStyle','-.'};
    for j = 1:size(par,1)-1
        plot([par(j,2),par(j,2)],[-maxi maxi],plotPar{:});
    end
    for j = 1:size(par,1)
        plot([par(j,1),par(j,2)],[par(j,3) par(j,3)],plotPar{:});
        plot([par(j,1),par(j,2)],-[par(j,3) par(j,3)],plotPar{:});                
    end
    ylim([-maxi*1.05 maxi*1.05]);
end
subplot(6,1,level+1);
xlabel('Time or Space');

Для каждого уровня k переменная thrParams{k} содержит интервалы и соответствующие пороги для процедуры шумоподавления.

Порог содействующий уровень уровнем вейвлета и интервал "интервал", с помощью значений содержится в переменной thrParams.

Используя функциональный wthresh, мы порог содействующие значения вейвлета между горизонтальными строками, заменяя их на нули, в то время как значения других или уменьшаются, если sorh = 's' или остаются неизменными если sorh = 'h'.

first = cumsum(longs)+1;
first = first(end-2:-1:1);
tmp   = longs(end-1:-1:2);
last  = first+tmp-1;
for k = 1:level
    thr_par = thrParams{k};
    if ~isempty(thr_par)
        cfs = coefs(first(k):last(k));
        nbCFS = longs(end-k);
        NB_int = size(thr_par,1);
        x = [thr_par(:,1) ; thr_par(NB_int,2)];
        alf = (nbCFS-1)/(x(end)-x(1));
        bet = 1 - alf*x(1);
        x = round(alf*x+bet);
        x(x<1) = 1;
        x(x>nbCFS) = nbCFS;
        thr = thr_par(:,3);
        for j = 1:NB_int
            if j == 1
                d_beg = 0; 
            else
                d_beg = 1;
            end
            j_beg = x(j)+ d_beg;
            j_end = x(j+1);
            j_ind = (j_beg:j_end);
            cfs(j_ind) = wthresh(cfs(j_ind),sorh,thr(j));
        end
        coefs(first(k):last(k)) = cfs;
    end
end

Покажите пороговые коэффициенты вейвлета сигнала.

% Replicate the coefficients.
cfs_beg = wrepcoef(coefs,longs);

% Display the decomposition coefficients.
clf
subplot(6,1,1); 
plot(sig,'r'); 
axis tight
title('Original Signal and Detail Coefficients from 1 to 5');
ylabel('S','Rotation',0);
for k = 1:level
    subplot(6,1,k+1);
    plot(cfs_beg(k,:),'Color',[0.5 0.8 0.5]);
    ylabel(['D' int2str(k)],'Rotation',0);
    axis tight
    hold on
    maxi = max(abs(cfs_beg(k,:)));
    hold on
    par = thrParams{k};
    plotPar = {'Color','m','LineStyle','-.'};
    for j = 1:size(par,1)-1
        plot([par(j,2),par(j,2)],[-maxi maxi],plotPar{:});
    end
    for j = 1:size(par,1)
        plot([par(j,1),par(j,2)], [par(j,3) par(j,3)],plotPar{:});
        plot([par(j,1),par(j,2)],-[par(j,3) par(j,3)],plotPar{:});                
    end
    ylim([-maxi*1.05 maxi*1.05]);
end
subplot(6,1,level+1);
xlabel('Time or Space');

Восстановите сигнал denoised.

sigden = waverec(coefs,longs,wname);
res = sig - sigden;

Отобразите оригинал, denoised, и остаточные сигналы.

clf
subplot(3,1,1);
plot(sig,'r');
hold on 
plot(sigden,'b'); 
axis tight
title('Original and Denoised Signals');
subplot(3,1,2);
plot(sigden,'b');  
axis tight
title('Denoised Signal');
subplot(3,1,3);
plot(res,'k');   
axis tight
title('Residual');

Сравните три denoised версии сигнала.

clf
plot(sigden_1,'g');
hold on 
plot(sigden_2,'r');
plot(sigden,'k');
axis tight
legend('Denoised Min','Denoised Max','Denoised IDT','Location','North');

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

xlim([1200 2000]);
ylim([180 350]);

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

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

Автоматическое вычисление зависимых интервалом порогов

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

% Wavelet Analysis.
[coefs,longs] = wavedec(sig,level,wname);
siz = size(coefs);
thrParams = utthrset_cmd(coefs,longs);
first = cumsum(longs)+1;
first = first(end-2:-1:1);
tmp   = longs(end-1:-1:2);
last  = first+tmp-1;
for k = 1:level
    thr_par = thrParams{k};
    if ~isempty(thr_par)
        cfs = coefs(first(k):last(k));
        nbCFS = longs(end-k);
        NB_int = size(thr_par,1);
        x = [thr_par(:,1) ; thr_par(NB_int,2)];
        alf = (nbCFS-1)/(x(end)-x(1));
        bet = 1 - alf*x(1);
        x = round(alf*x+bet);
        x(x<1) = 1;
        x(x>nbCFS) = nbCFS;
        thr = thr_par(:,3);
        for j = 1:NB_int
            if j==1 , 
                d_beg = 0;
            else
                d_beg = 1;
            end
            j_beg = x(j)+d_beg;
            j_end = x(j+1);
            j_ind = (j_beg:j_end);
            cfs(j_ind) = wthresh(cfs(j_ind),sorh,thr(j));
        end
        coefs(first(k):last(k)) = cfs;
    end
end
sigden = waverec(coefs,longs,wname);

clf
subplot(2,1,1); 
plot(sig,'r');
axis tight
hold on
plot(sigden,'k');
title('Original and Denoised Signals');
subplot(2,1,2);
plot(sigden,'k');
axis tight
title('Denoised Signal');

Автоматическое зависимое интервалом шумоподавление

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

[sigden,~,thrParams] = cmddenoise(sig,wname,level);
thrParams{1}   % Denoising parameters for level 1.

ans = 2×3
10³ ×

    0.0010    1.1100    0.0176
    1.1100    2.0000    0.0045

Автоматическая процедура находит два интервала для шумоподавления:

I1 = [1 1110] с порогом thr1 = 17.6
I2 = [1110 2000] с порогом thr2 = 4.5.

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

clf
subplot(2,1,1);
plot(sig,'r');
axis tight
hold on
plot(sigden,'k');
title('Original and Denoised Signals');
subplot(2,1,2);
plot(sigden,'k');
axis tight
title('Denoised Signal');

Усовершенствованное автоматическое зависимое интервалом шумоподавление

Теперь, мы смотрим на более полный пример автоматического шумоподавления.

Вместо того, чтобы использовать значения по умолчанию для входных параметров, мы можем задать их при вызывании функции. Здесь, тип порога выбран в качестве (мягкого) s, и номер интервалов определяется в 3.

load nelec.mat;
sig = nelec;   % Signal to analyze.
wname  = 'sym4';               % Wavelet for analysis.
level  = 5;                    % Level for wavelet decomposition.
sorh   = 's';                  % Type of thresholding.
nb_Int = 3;                    % Number of intervals for thresholding.

[sigden,coefs,thrParams,int_DepThr_Cell,BestNbOfInt] = ...
            cmddenoise(sig,wname,level,sorh,nb_Int);

Для выходных параметров переменная thrParams{1} дает параметры шумоподавления для уровней от 1 до 5. Например, вот параметры шумоподавления для уровня 1.

thrParams{1}

ans = 3×3
10³ ×

    0.0010    0.0940    0.0059
    0.0940    1.1100    0.0195
    1.1100    2.0000    0.0045

Мы находим те же значения, которые были установлены ранее этот пример. Они соответствуют выбору, который мы сделали путем фиксации количества интервалов к три во входном параметре: nb_Int = 3.

Автоматическая процедура предлагает 2 как лучшее количество интервалов для шумоподавления. Это выходное значение BestNbOfInt = 2 эквивалентно используемый на предыдущем шаге этого примера.

BestNbOfInt

BestNbOfInt = 2

Переменная int_DepThr_Cell содержит местоположения интервала и пороговые значения для многих интервалов от 1 до 6.

int_DepThr_Cell

int_DepThr_Cell = 1x6 cell array
  Columns 1 through 4

    {1x3 double}    {2x3 double}    {3x3 double}    {4x3 double}

  Columns 5 through 6

    {5x3 double}    {6x3 double}

Наконец, мы смотрим на значения, соответствующие местоположениям и порогам для 5 интервалов.

int_DepThr_Cell{5}

ans = 5×3
10³ ×

    0.0010    0.0940    0.0059
    0.0940    0.5420    0.0184
    0.5420    0.5640    0.0056
    0.5640    1.1100    0.0240
    1.1100    2.0000    0.0045

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

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

Документация