Мягкий или жесткий порог?

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

thr = 0.4;

Жесткое пороговое значение может быть описано как обычный процесс настройки для нуля элементов, абсолютные значения которых ниже порога. Жесткий пороговый сигнал x, если x > thr, и равен 0, если x $\le$ thr.

y = linspace(-1,1,100); 
ythard = wthresh(y,'h',thr);

Мягкое пороговое значение является расширением жесткого порога, сначала устанавливая, чтобы нуль элементы, абсолютные значения которых ниже порога, а затем сжимая ненулевые коэффициенты к 0. Мягкий пороговый сигнал является sign (x) (x - thr), если x > thr и равен 0, если x $\le$ thr.

ytsoft = wthresh(y,'s',thr);
subplot(1,3,1)
plot(y)
title('Original')
subplot(1,3,2)
plot(ythard)
title('Hard Thresholding')
subplot(1,3,3)
plot(ytsoft)
title('Soft Thresholding')

Как видно на рисунке выше, жесткая процедура создает разрывы в x = $$\pm$$ t, в то время как мягкая процедура не делает.

Правила выбора порогов

Ссылаясь на шаг 2 процедуры denoise, функция thselect выполняет выбор порога и затем каждый уровень устанавливается. Этот второй шаг может быть выполнен с помощью wthcoeff, непосредственно обрабатывая структуру вейвлета разложения исходного сигнала. В функции реализованы четыре правила выбора пороговых значений thselect. Обычно интересно показать их в действии, когда входной сигнал является Гауссовым белым шумом.

rng default
y = randn(1,1000);

Правило 1: Выбор с использованием принципа объективной оценки риска Штейна (SURE)

thr = thselect(y,'rigrsure')

thr = 2.0518

Правило 2: Фиксированный порог формы, равный sqrt (2 * журнал (length (y)))

thr = thselect(y,'sqtwolog')

thr = 3.7169

Правило 3: Выбор с использованием смеси первых двух опций

thr = thselect(y,'heursure')

thr = 3.7169

Правило 4: Выбор по минимаксному принципу

thr = thselect(y,'minimaxi')

thr = 2.2163

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

Давайте будем использовать «блоки» тестовых данных, зачисленных в Donoho и Johnstone в качестве первого примера. Сгенерируйте исходный сигнал внешней ссылки и шумную версию x, добавляя стандартный Гауссов белый шум.

sqrt_snr = 4;      % Set signal to noise ratio
init = 2055615866; % Set rand seed
[xref,x] = wnoise(1,11,sqrt_snr,init);

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

xd = wdenoise(x);
subplot(3,1,1)
plot(xref)
axis tight
title('Original Signal')
subplot(3,1,2)
plot(x)
axis tight
title('Noisy Signal')
subplot(3,1,3)
plot(xd)
axis tight
title('Denoised Signal - Signal to Noise Ratio = 4')

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

xd2 = wdenoise(x,3,'Wavelet','sym8',...
    'DenoisingMethod','SURE',...
    'ThresholdRule','Soft');
figure
subplot(3,1,1)
plot(xref)
axis tight
title('Original Signal')
subplot(3,1,2)
plot(xd)
axis tight
title('First Denoised Signal: Default Settings')
subplot(3,1,3)
plot(xd2)
axis tight
title('Second Denoised Signal')

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

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

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

load leleccum
indx = 2000:3450; 
x = leleccum(indx); % Load electrical signal and select part of it.

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

xd = wdenoise(x,3,'Wavelet','db3',...
    'DenoisingMethod','UniversalThreshold',...
    'ThresholdRule','Soft',...
    'NoiseEstimate','LevelDependent');
Nx = length(x);
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')

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

Шумоподавление изображений

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

Сгенерируйте шумное изображение.

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

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

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

thr = 107.9838

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

xd = wdencmp('gbl',x,'sym4',2,thr,sorh,keepapp);
figure('Color','white')
colormap(pink(255))
sm = size(map,1);
image(wcodemat(X,sm))
title('Original Image')

figure('Color','white')
colormap(pink(255))
image(wcodemat(x,sm))
title('Noisy Image')

image(wcodemat(xd,sm))
title('Denoised Image')

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

Шумоподавление, основанные на вейвлет-разложении, являются одним из наиболее значимых применений вейвлеты. См. wdenoise и Wavelet Signal Denoiser для получения дополнительной информации.