Оценка вейвлет- Шумоподавление и непараметрическая функция

Wavelet Toolbox™ обеспечивает ряд функций для оценки неизвестной функции (сигнала или изображения) в шуме. Можно использовать эти функции для денуизации сигналов и как метод для непараметрической оценки функции.

Самая общая модель 1-D для этого является

s (n) = f (n) +, e (n)

где n = 0,1,2,...N-1. e (n) являются Гауссовыми случайными переменными, распределенными как N (0,1). Отклонение и e (n)².

На практике s (n) часто является сигналом дискретного времени с равными временными шагами, поврежденными аддитивным шумом, и вы пытаетесь восстановить этот сигнал.

В более общем случае можно просмотреть s (n) как N-мерный случайный вектор

$(\begin{matrix} f (0) + σ e (0) \\ f (1) + σ e (1) \\ f (2) + σ e (2) \\ . \\ . \\ . \\ f (N - 1) + σ e (N - 1) \end{matrix}) = (\begin{matrix} f (0) \\ f (1) \\ f (2) \\ . \\ . \\ . \\ f (N - 1) \end{matrix}) + (\begin{matrix} σ e (0) \\ σ e (1) \\ σ e (2) \\ . \\ . \\ . \\ σ e (N - 1) \end{matrix})$

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

Можно заменить N-by-1 случайный вектор на N на M случайных матриц, чтобы получить задачу восстановления изображения, поврежденного аддитивным шумом.

Вы можете получить 1-D пример этой модели со следующим кодом.

load cuspamax;
y = cuspamax+0.5*randn(size(cuspamax));
plot(y); hold on;
plot(cuspamax,'r','linewidth',2);
axis tight;
legend('f(n)+\sigma e(n)','f(n)', 'Location', 'NorthWest');

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

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

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

Вы можете результировать следующие шаги как:

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

Шумоподавление

Wavelet Toolbox поддерживает ряд методов шумоподавления. Четыре метода шумоподавления реализованы в thselect. Каждый метод соответствует tptr опция в команде

thr = thselect(y,tptr)

который возвращает пороговое значение.

Опция	Метод шумоподавления
`'rigrsure'`	Выбор с использованием принципа объективной оценки риска Штейна (SURE)
`'sqtwolog'`	Фиксированный порог формы (универсальный), равный $\sqrt{2 \ln (N)}$ с N длины сигнала.
`'heursure'`	Выбор с использованием смеси первых двух опций
`'minimaxi'`	Выбор по принципу minimax

Опция 'rigrsure' использует для оценщика мягкого порога правило выбора порога, основанное на объективной оценке риска Штейна (функция квадратичных потерь). Вы получаете оценку риска для конкретного порога значения t. Минимизация рисков в t дает выбор порогового значения.
Опция 'sqtwolog' использует фиксированный порог формы, дающий минимаксную эффективность, умноженную на небольшой коэффициент, пропорциональный логарифмической (длина (s)).
Опция 'heursure' представляет собой смесь двух предыдущих опций. В результате, если отношение сигнал/шум очень мало, оценка SURE очень шумна. Поэтому, если такая ситуация обнаружена, используется фиксированный порог формы.
Опция 'minimaxi' использует фиксированный порог, выбранный для получения минимаксной эффективности для средней квадратной ошибки против идеальной процедуры. Принцип минимакса используется в статистике для проектировщиков. Поскольку деноизированный сигнал может быть ассимилирован в оценщик неизвестной регрессионой функции, минимаксная оценка является опцией, которая реализует минимум, по заданному набору функций, максимальной средней квадратной ошибки.

Следующий пример показывает методы шумоподавления для вектора 1000 на 1 N (0,1). Сигнал здесь

$f (n) + e (n) e (n) ~ N (0, 1)$

с f(n) = 0.

rng default;
sig = randn(1e3,1);
thr_rigrsure = thselect(sig,'rigrsure')
thr_univthresh = thselect(sig,'sqtwolog')
thr_heursure = thselect(sig,'heursure')
thr_minimaxi = thselect(sig,'minimaxi')
histogram(sig);
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[0.7 0.7 0.7],'EdgeColor','w');
hold on;
plot([thr_rigrsure thr_rigrsure], [0 300],'linewidth',2);
plot([thr_univthresh thr_univthresh], [0 300],'r','linewidth',2);
plot([thr_minimaxi thr_minimaxi], [0 300],'k','linewidth',2);
plot([-thr_rigrsure -thr_rigrsure], [0 300],'linewidth',2);
plot([-thr_univthresh -thr_univthresh], [0 300],'r','linewidth',2);
plot([-thr_minimaxi -thr_minimaxi], [0 300],'k','linewidth',2);

Для объективной оценки риска (SURE) и минимаксных порогов Штейна сохраняются приблизительно 3% коэффициентов. В случае универсального порога все значения отклоняются.

Мы знаем, что вектор коэффициентов детализации является суперпозицией коэффициентов f и коэффициентов e, и что разложение e приводит к коэффициентам детализации, которые являются стандартными Гауссовыми белыми шумами.

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

Мягкое или жесткое пороговое значение

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

Самая простая схема - жесткое пороговое значение. Позвольте T обозначить порог и x свои данные. Жесткое пороговое значение

$η (x) = {\begin{array}{l} x & | x | \geq T \\ 0 & | x | < T \end{array}$

Мягкое пороговое значение

$η (x) = {\begin{array}{l} x - T & x > T \\ 0 & | x | \leq T \\ x + T & x < - T \end{array}$

Вы можете применить свой порог, используя жесткое или мягкое правило с wthresh.

y = linspace(-1,1,100);
thr = 0.4;
ythard = wthresh(y,'h',thr);
ytsoft = wthresh(y,'s',thr);
subplot(131);
plot(y); title('Original Data');
subplot(132);
plot(ythard,'*'); title('Hard Thresholding');
subplot(133);
plot(ytsoft,'*'); title('Soft Thresholding');

Борьба с немасштабированным шумом и небелым шумом

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

Самое простое использование wdenoise является

sd = wdenoise(s)

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

sd = wdenoise(s,n,'DenoisingMethod',tptr,'Wavelet',wav,...
     'ThresholdRule',sorh,'NoiseEstimate',scal)

который возвращает деноизированную версию sd исходного сигнала s полученный с использованием tptr метод шумоподавления. Другие необходимые параметры sorh, scal, и wname. Значение параметра sorh задает пороговое значение коэффициентов детализации разложения на уровне n от s по вейвлету, называемому wav. Оставшийся параметр scal должно быть задано. Это соответствует способу оценки отклонения шума в данных.

Опция	Метод оценки шума
`'LevelIndependent'`	`'LevelIndependent'` оценивает отклонение шума на основе коэффициентов вейвлет самой мелкой шкалы (с самым высоким разрешением).
`'LevelDependent'`	`'LevelDependent'` оценивает отклонение шума на основе коэффициентов вейвлета на каждом уровне разрешения.

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

 xd = wdencmp(opt,x,wav,n,thr,sorh,keepapp);

где

opt = 'gbl' и thr является положительным вещественным числом для равномерного порога.
opt = 'lvd' и thr является вектором для зависящего от уровня порога.
keepapp = 1 сохранить аппроксимационные коэффициенты, как ранее и
keepapp = 0 чтобы разрешить пороговое значение коэффициентов приближения.
x - сигнал, который будет денофицирован и wav, n, sorh являются такими же, как и выше.

Вейвлет Шумоподавления в действии

Открыть Live Script

Мы начинаем примеры 1-D методов шумоподавления с первого примера, приписываемого Донохо и Джонстоуну.

Блокирует пороговое значение сигнала

Сначала установите отношение сигнал/шум (ОСШ) и установите случайный seed.

sqrt_snr = 4;
init = 2055615866;

Сгенерируйте исходный сигнал xref и шумную версию x путем добавления стандартного Гауссова белого шума. Постройте график обоих сигналов.

[xref,x] = wnoise(1,11,sqrt_snr,init);
subplot(2,1,1)
plot(xref)
axis tight
title('Original Signal')
subplot(2,1,2)
plot(x)
axis tight
title('Noisy Signal')

Figure contains 2 axes. Axes 1 with title Original Signal contains an object of type line. Axes 2 with title Noisy Signal contains an object of type line.

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

xd = wdenoise(x,'Wavelet','sym8','DenoisingMethod','SURE','ThresholdRule','Soft');
figure
subplot(2,1,1)
plot(xref)
axis tight
title('Original Signal')
subplot(2,1,2)
plot(xd)
axis tight
title('Denoised Signal')

Figure contains 2 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.

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

Шумоподавление электрического сигнала

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

Сначала загрузите электрический сигнал и выберите из него сегмент. Постройте график сегмента.

load leleccum
indx = 2000:3450;
x = leleccum(indx);
figure
plot(indx,x)
axis tight
title('Original Signal')

Figure contains an axes. The axes with title Original Signal contains an object of type line.

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

xd = wdenoise(x,3,'Wavelet','db3',...
    'DenoisingMethod','UniversalThreshold',...
    'ThresholdRule','Soft',...
    'NoiseEstimate','LevelDependent');
figure
subplot(2,1,1)
plot(indx,x)
axis tight
title('Original Signal')
subplot(2,1,2)
plot(indx,xd)
axis tight
title('Denoised Signal')

Figure contains 2 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.

Результат довольно хорош, несмотря на временную неоднородность характера шума после и до начала отказа датчика около времени 2410.

Расширение для шумоподавления изображений

Открыть Live Script

Способ шумоподавления, описанный для 1-D случая, применяется также к изображениям и хорошо применяется к геометрическим изображениям. Прямым переводом модели 1-D является

$s (i, j) = f (i, j) + σ e (i, j)$

где $e$ является белым Гауссовым шумом с единичным отклонением.

Процедура 2-D шумоподавления имеет те же три шага и использует 2-D вейвлеты инструменты вместо 1-D инструментов. Для выбора порога prod(size(s)) используется вместо length(s) если используется фиксированный порог формы.

Обратите внимание, что кроме «автоматического» случая 1-D деноизации, 2-D деноизации и сжатия выполняются с помощью wdencmp. Чтобы проиллюстрировать 2-D шумоподавление, загрузите изображение и создайте его шумную версию. В целях воспроизводимости установите случайный seed.

init = 2055615866;
rng(init);
load woman
img = X;
imgNoisy = img + 15*randn(size(img));

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

[thr,sorh,keepapp] = ddencmp('den','wv',imgNoisy);
thr

thr = 107.9838

thr равно estimated_sigma*sqrt(log(prod(size(img)))).

Обесценивайте шумное изображение с помощью опции глобального порога. Отображение результатов.

imgDenoised = wdencmp('gbl',imgNoisy,'sym4',2,thr,sorh,keepapp);
figure
colormap(pink(255))
sm = size(map,1);
subplot(2,2,1)
image(wcodemat(img,sm))
title('Original Image')
subplot(2,2,2)
image(wcodemat(imgNoisy,sm))
title('Noisy Image')
subplot(2,2,3)
image(wcodemat(imgDenoised,sm))
title('Denoised Image')

Figure contains 3 axes. Axes 1 with title Original Image contains an object of type image. Axes 2 with title Noisy Image contains an object of type image. Axes 3 with title Denoised Image contains an object of type image.

Облицованное изображение хорошо сравнивается с оригинальным изображением.

1-D адаптивного порога отклонения вейвлета

Открыть Live Script

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

Точнее, модель принимает (как и ранее), что наблюдение равно интересному сигналу, наложенному на шум

$s (n) = f (n) + σ e (n)$

Но отклонение шума может варьироваться со временем. Существует несколько различных значений отклонения на нескольких временных интервалах. Значения, а также интервалы неизвестны.

Давайте сосредоточимся на проблеме оценки точек изменения или эквивалентно интервалам. Используемый алгоритм основан на исходной работе Марка Лавилля об обнаружении точек изменения с помощью динамического программирования (см. [Lav99] в ссылки).

Давайте сгенерируем сигнал от регрессионной модели с фиксированной конструкцией с двумя точками изменения отклонения шума, расположенными в положениях 200 и 600. В целях воспроизводимости установите случайный seed.

init = 2055615866;
rng(init);

x = wnoise(1,10);
bb = randn(1,length(x));
cp1 = 200;
cp2 = 600;
x = x+[bb(1:cp1),bb(cp1+1:cp2)/4,bb(cp2+1:end)];
plot(x)
title('Noisy Signal')

Figure contains an axes. The axes with title Noisy Signal contains an object of type line.

Цель этого примера состоит в том, чтобы восстановить две точки изменения из сигнала x.

Шаг 1. Восстановите сигнал с шумом путем подавления приближения. Сначала выполните одноуровневое разложение вейвлета, используя db3 вейвлет. Затем восстановите деталь на уровне 1.

wname = 'db3';
lev = 1;
[c,l] = wavedec(x,lev,wname);
det = wrcoef('d',c,l,wname,1);
figure
plot(det)
title('Level 1 Detail')

Figure contains an axes. The axes with title Level 1 Detail contains an object of type line.

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

Шаг 2. Чтобы удалить почти весь сигнал, замените 2% самых больших значений средним значением.

x = sort(abs(det));
v2p100 = x(fix(length(x)*0.98));
ind = find(abs(det)>v2p100);
det(ind) = mean(det);

Шаг 3. Используйте wvarchg функция для оценки точек изменения со следующими параметрами:

Минимальная задержка между двумя точками изменения d = 10.
Максимальное количество точек изменения - 5.

[cp_est,kopt,t_est] = wvarchg(det,5)

cp_est = 1×2

   259   611

kopt = 2

t_est = 6×6

        1024           0           0           0           0           0
         612        1024           0           0           0           0
         259         611        1024           0           0           0
         198         259         611        1024           0           0
         198         235         259         611        1024           0
         198         235         260         346         611        1024

Предлагаются две точки изменения и три интервала. Поскольку три отклонения интервала для шума очень различны, программа оптимизации легко обнаруживает правильную структуру. Предполагаемые точки изменения близки к истинным точкам изменения: 200 и 600.

Шаг 4. (Необязательно) Замените предполагаемые точки изменения.

Для 2 $\leq$ i $\leq$ 6, t_est(i,1:i-1) содержит i-1 моменты точек изменения отклонения, и с kopt предлагаемое количество точек изменения, затем

cp_est = t_est(kopt+1,1:kopt);

Можно заменить предполагаемые точки изменения вычислением:

for k=1:5
    cp_New = t_est(k+1,1:k)
end

cp_New = 612

cp_New = 1×2

   259   611

cp_New = 1×3

   198   259   611

cp_New = 1×4

   198   235   259   611

cp_New = 1×5

   198   235   260   346   611

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

Следующие измерения и настройки полезны для анализа вейвлет и изображений:

M S E - Средняя квадратная ошибка (MSE) является квадратной нормой различия между данными и приближением сигнала или изображения, разделенным на количество элементов.
Max Error - Максимальное абсолютное квадратное отклонение в приближении сигнала или изображения.
L2-Norm Ratio - Отношение квадратов L2-norm сигнала или приближения изображения к входному сигналу или изображению. Для изображений изображение изменяется как вектор-столбец перед принятием L2-norm
P S N R - отношение пикового сигнала к шуму (PSNR) в децибелах. PSNR имеет значение только для данных, закодированных с точки зрения бит на выборку или биты на пиксель.
B P P - Биты на пиксель (BPP), которое является степенью сжатия (Comp. Ratio), умноженной на 8, принимая один байт на пиксель (8 биты).
Comp Ratio - Коэффициент сжатия, который является количеством элементов в сжатом изображении, разделенным на количество элементов в оригинальном изображении, выраженное в процентах.

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