信号とイメージのノイズ除去

ソフトなしきい値処理またはハードなしきい値処理

しきい値処理は、入力信号のソフトなしきい値処理またはハードなしきい値処理を返す関数 wthresh を使用して実行できます。ハードなしきい値処理は最も簡潔な方法ですが、ソフトなしきい値処理には優れた数学的性質があります。thr がしきい値を表します。

thr = 0.4;

ハードなしきい値処理は、絶対値がしきい値未満である要素を 0 に設定する通常の設定処理であると説明できます。x > thr である場合のハードなしきい値信号は x であり、x $\le$ thr である場合は 0 です。

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

ソフトなしきい値処理はハードなしきい値処理を拡張したものです。最初に、絶対値がしきい値未満である要素を 0 に設定し、次に非ゼロの係数を 0 に向かって縮小していきます。x > thr である場合のソフトなしきい値信号は sign(x)(x-thr) であり、x $\le$ thr である場合は 0 です。

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 を呼び出すと、関数 thselect によってしきい値の選択が実行され、各レベルがしきい値処理されます。この 2 番目のステップは wthcoeff を使用して、元の信号のウェーブレット分解構造を直接処理することで実行できます。関数 thselect には 4 つのしきい値選択ルールが実装されています。通常は、入力信号がガウス ホワイト ノイズである場合に実施されているルールを示すと興味深いものになります。

rng default
y = randn(1,1000);

ルール 1: Stein の不偏リスク推定 (SURE) の原理を使用した選択

thr = thselect(y,'rigrsure')

thr = 2.0518

ルール 2: sqrt(2*log(length(y))) に等しい固定型のしきい値

thr = thselect(y,'sqtwolog')

thr = 3.7169

ルール 3: 最初の 2 つのオプションを混合させて使用した選択

thr = thselect(y,'heursure')

thr = 3.7169

ルール 4: ミニマックスの原理を使用した選択

thr = thselect(y,'minimaxi')

thr = 2.2163

ミニマックスと SURE のしきい値選択ルールはより保守的であり、信号の小さな Detail がノイズ範囲の付近にある場合に便利です。他の 2 つのルールはノイズをより効率的に除去します。

Donoho と Johnstone の功績である "blocks" テスト データを最初の例として使用します。元の信号 xref と、標準的なガウス ホワイト ノイズを追加するノイズが含まれるバージョン 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')

ノイズ信号の 2 回目のノイズ除去を行います。今回は、sym8 ウェーブレットによりレベル 3 で x の分解から取得した Detail 係数を経験的な SURE のソフトなしきい値処理を使用します。以前のノイズ除去後の信号と比較します。

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

少数の大きい係数のみが元の信号を特徴付けるため、ノイズ除去後の両方の信号を元の信号とよく比較します。ウェーブレット信号デノイザーを使用して、他のノイズ除去パラメーターがノイズ信号に与える影響を調べることができます。

非ホワイト ノイズの処理

非ホワイト ノイズが疑わしい場合、そのレベルのノイズのレベル依存推定によってしきい値を再スケーリングしなければなりません。2 番目の例として、電気信号の強く摂動を受ける部分に対してメソッドを試します。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 あたりのセンサー障害が始まる前後で、ノイズの性質に時間の異質性があるにもかかわらず、結果は非常に良好です。

イメージのノイズ除去

1 次元のケースに対して説明したノイズ除去方法はイメージにも適用され、幾何学的イメージに対してもうまく当てはまります。2 次元のノイズ除去手順も同様に 3 ステップで構成されていて、1 次元ウェーブレット ツールの代わりに 2 次元ウェーブレット ツールを使用します。しきい値の選択については、固定型のしきい値が使用されている場合は、length(y) の代わりに prod(size(y)) が使用されます。

ノイズを含むイメージを生成します。

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

この場合、固定型のしきい値がレベル ノイズの推定で使用され、しきい値処理モードはソフトで、Approximation 係数は維持されます。

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

まとめ

ウェーブレット分解に基づくノイズ除去方法は、ウェーブレットの最も重要な用途の 1 つです。詳細については、wdenoiseおよびウェーブレット信号デノイザーを参照してください。