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

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

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

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

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

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

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

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

load nelec.mat
sig = nelec;

Теперь мы используем дискретный анализ вейвлета на уровне 5 с sym4 вейвлетом. Мы устанавливаем тип пороговой обработки на 's' (мягкий).

wname = 'sym4';
level = 5;
sorh  = 's';

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')

Figure contains 3 axes objects. Axes object 1 with title Original Signal contains an object of type line. Axes object 2 with title Denoised Signal contains an object of type line. Axes object 3 with title Residual contains an object of type line.

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 шумоподавления.

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

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

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

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

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

denPAR = {[1 94 5.9 ; 94 1110 19.5 ; 1110 2000 4.5]};
thrParams = cell(1,level);
thrParams(:) = denPAR;

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

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

% Display the coefficients of the decomposition
figure
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])
    %hold off
end
subplot(6,1,level+1)
xlabel('Time or Space')

Figure contains 6 axes objects. Axes object 1 with title Original Signal and Detail Coefficients from 1 to 5 contains an object of type line. Axes object 2 contains 9 objects of type line. Axes object 3 contains 9 objects of type line. Axes object 4 contains 9 objects of type line. Axes object 5 contains 9 objects of type line. Axes object 6 contains 9 objects of type line.

Для каждого уровня 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.
figure
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])
    hold off
end
subplot(6,1,level+1)
xlabel('Time or Space')

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

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

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

figure
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')

Figure contains 3 axes objects. Axes object 1 with title Original and Denoised Signals contains 2 objects of type line. Axes object 2 with title Denoised Signal contains an object of type line. Axes object 3 with title Residual contains an object of type line.

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

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

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent Denoised Min, Denoised Max, Denoised IDT.

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

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

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent Denoised Min, Denoised Max, Denoised IDT.

Мы видим, что, когда максимальное пороговое значение используется, сигнал 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);
figure
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
hold off
title('Denoised Signal')

Figure contains 2 axes objects. Axes object 1 with title Original and Denoised Signals contains 2 objects of type line. Axes object 2 with title Denoised Signal contains an object of type line.

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

В режиме командной строки мы можем использовать функциональный 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.

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

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

Figure contains 2 axes objects. Axes object 1 with title Original and Denoised Signals contains 2 objects of type line. Axes object 2 with title Denoised Signal contains an object of type line.

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

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

Вместо того, чтобы использовать значения по умолчанию для входных параметров, мы можем задать их при вызывании функции. Здесь, тип порога выбран в качестве 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=1×6 cell array
  Columns 1 through 4

    {[1 2000 8.3611]}    {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 для шумоподавления при давании вам большего количества контроля конкретными значениями параметров, чтобы получить лучшие результаты.

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

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

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

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

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

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

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

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

Смотрите также

Похожие темы

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

Поддержка