Main Content

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

この例では、ノイズの多いデータからの信号再生に関する問題について説明します。一般的なノイズ除去手順には 3 つのステップがあります。基本的なバージョンの手順では、以下で説明するステップに従います。

  • 分解: ウェーブレットを選択し、レベル N を選択します。レベル N の信号のウェーブレット分解を計算します。

  • Detail 係数のしきい値処理: 1 から N までの各レベルで、しきい値を選択し、Detail 係数にソフトなしきい値処理を適用します。

  • 再構成: レベル N の元の Approximation 係数と、レベル 1 から N までの変更済みの Detail 係数を使用して、ウェーブレットの再構成を計算します。

特に、次の 2 点に対処しなければなりません。

  • しきい値の選択方法

  • しきい値処理の実行方法

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

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

thr = 0.4;

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

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

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

Figure contains 3 axes objects. Axes object 1 with title Original contains an object of type line. Axes object 2 with title Hard Thresholding contains an object of type line. Axes object 3 with title Soft Thresholding contains an object of type line.

上の図でわかるように、ハードな手順では x = ± 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')

Figure contains 3 axes objects. Axes object 1 with title Original Signal contains an object of type line. Axes object 2 with title Noisy Signal contains an object of type line. Axes object 3 with title Denoised Signal - Signal to Noise Ratio = 4 contains an object of type line.

ノイズ信号の 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')

Figure contains 3 axes objects. Axes object 1 with title Original Signal contains an object of type line. Axes object 2 with title First Denoised Signal: Default Settings contains an object of type line. Axes object 3 with title Second Denoised Signal contains an object of type line.

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

非ホワイト ノイズの処理

非ホワイト ノイズが疑わしい場合、そのレベルのノイズのレベル依存推定によってしきい値を再スケーリングしなければなりません。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')

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

時間 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 contains an axes object. The axes object with title Original Image contains an object of type image.

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

Figure contains an axes object. The axes object with title Noisy Image contains an object of type image.

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

Figure contains an axes object. The axes object with title Denoised Image contains an object of type image.

まとめ

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

参考

|

関連するトピック