Вейвлет шумоподавление

Пример показывает, как обесценить сигнал, используя интервально-зависимые пороги.

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

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

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

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

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

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

load nelec.mat
sig = nelec;

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

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

Денуризируйте сигнал с помощью функции 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. Axes 1 with title Original Signal contains an object of type line. Axes 2 with title Denoised Signal contains an object of type line. Axes 3 with title Residual contains an object of type line.

perf0,perfl2

perf0 = 66.6995

perfl2 = 99.9756

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

Теперь мы обесцениваем сигнал с максимальным значением, предоставленным инструментами 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

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

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

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

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

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

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

I1 = [1 94] с порогом tr1 = 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. Axes 1 with title Original Signal and Detail Coefficients from 1 to 5 contains an object of type line. Axes 2 contains 9 objects of type line. Axes 3 contains 9 objects of type line. Axes 4 contains 9 objects of type line. Axes 5 contains 9 objects of type line. Axes 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')

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

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

Отображение исходных, деноизированных и остаточных сигналов.

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. Axes 1 with title Original and Denoised Signals contains 2 objects of type line. Axes 2 with title Denoised Signal contains an object of type line. Axes 3 with title Residual contains an object of type line.

Сравните три деноизированные версии сигнала.

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. The axes 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. The axes contains 3 objects of type line. These objects represent Denoised Min, Denoised Max, Denoised IDT.

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

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

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

Вместо того, чтобы вручную устанавливать интервалы и пороги для каждого уровня, мы можем использовать функцию 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. Axes 1 with title Original and Denoised Signals contains 2 objects of type line. Axes 2 with title Denoised Signal contains an object of type line.

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

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

[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. Axes 1 with title Original and Denoised Signals contains 2 objects of type line. Axes 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 для шумоподавления, в то же время предоставляя вам больше контроля над конкретными значениями параметров, чтобы получить лучшие результаты.

См. также

Wavelet Signal Denoiser | wdenoise | wdenoise2

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