主要内容

小波去噪与非参数函数估计

小波工具箱™提供了许多函数,用于在噪声中估计未知函数(信号或图像)。您可以使用这些函数去噪信号,并作为非参数函数估计的方法。

最常用的一维模型是

年代n) =fn) +σen

在哪里<年代pan class="inlineequation">n= 0, 1, 2,…n - 1。的en)为高斯随机变量分布为N(0, 1)。σ的方差en)是σ<年代up>2.

在实践中,年代n)通常是一个等时间步长的离散信号,被加性噪声破坏,而你正试图恢复该信号。

更一般地说,你可以查看年代n)作为n维随机向量

f 0 + σ e 0 f 1 + σ e 1 f 2 + σ e 2 f N 1 + σ e N 1 f 0 f 1 f 2 f N 1 + σ e 0 σ e 1 σ e 2 σ e N 1

在一般情况下,去噪和回归之间的关系是清楚的。

你可以用n × m的随机矩阵替换n × 1的随机向量来获得被加性噪声损坏的图像的恢复问题。

您可以使用以下代码获得该模型的一维示例。

负载<年代pan style="color:#A020F0">cuspamax;y = cuspamax + 0.5 * randn(大小(cuspamax));情节(y);持有<年代pan style="color:#A020F0">在;情节(cuspamax<年代pan style="color:#A020F0">“r”,<年代pan style="color:#A020F0">“线宽”2);轴<年代pan style="color:#A020F0">紧;传奇(<年代pan style="color:#A020F0">f (n) + \σe (n) ',<年代pan style="color:#A020F0">“f (n)”,<年代pan style="color:#A020F0">“位置”,<年代pan style="color:#A020F0">“西北”);

对于一类具有一定平滑特性的函数(信号、图像),小波技术对于函数恢复是最优或接近最优的。

具体来说,该方法对函数族是有效的f只有几个非零小波系数。这些函数具有稀疏小波表示。例如,一个几乎处处平滑的函数,只有少数突变,就具有这样的性质。

一般的基于小波的去噪和非参数函数估计方法是将数据变换到小波域,对小波系数进行阈值化,然后进行逆变换。

你可以把这些步骤总结为:

  1. 分解

    选择小波和电平N.计算信号的小波分解年代到水平N

  2. 阈值的细节系数

    从1到N,阈值细节系数。

  3. 重建

    利用原水平近似系数计算小波重构N以及从1到N

去噪方法

小波工具箱支持多种去噪方法。金宝app实现了四种去噪方法thselect.每个方法对应于atptr命令中的选项

用力推= thselect (y, tptr)

返回阈值。

选项

去噪方法

“rigrsure”

运用斯坦的无偏风险估计原则进行选择

“sqtwolog”

固定形式(通用)阈值等于

2 ln N

N信号的长度。

“heursure”

使用前两个选项的混合选择

“minimaxi”

选择使用极大极小原理

  • 选项“rigrsure”软阈值估计使用基于斯坦的无偏风险估计(二次损失函数)的阈值选择规则。你会得到一个特定阈值的风险估计t.将风险降到最低t给出阈值的选择。

  • 选项“sqtwolog”使用一个固定形式的阈值产生的极小极大的性能乘以一个小的因素成正比<年代pan class="emphasis">日志(<年代pan class="emphasis">长度年代))。

  • 选项“heursure”是前两种选择的混合。因此,在信噪比很小的情况下肯定的估计是非常嘈杂的。因此,如果检测到这种情况,则使用固定形式阈值。

  • 选项“minimaxi”使用选择屈服的固定阈值对均方误差的极小极大性能的理想程序。在统计学中,极小极大原理被用来设计估计量。由于去噪信号可以被同化为未知回归函数的估计器,极小极大估计器是实现在给定函数集上的最大均方误差最小的选项。

下面的例子展示了一个1000 × 1 N(0,1)向量的去噪方法。这里的信号是

f n + e n e n N 0 1

与<年代pan class="inlineequation">f (n)= 0.

rng<年代pan style="color:#A020F0">默认的;sig = randn (e3, 1);thr_rigrsure = thselect(团体,<年代pan style="color:#A020F0">“rigrsure”thr_univthresh = thselect(sig,<年代pan style="color:#A020F0">“sqtwolog”) thr_heursure = thselect(sig,<年代pan style="color:#A020F0">“heursure”) thr_minimaxi = thselect(sig,<年代pan style="color:#A020F0">“minimaxi”)直方图(团体);h = findobj (gca),<年代pan style="color:#A020F0">“类型”,<年代pan style="color:#A020F0">“补丁”);集(h,<年代pan style="color:#A020F0">“FaceColor”(0.7 0.7 0.7),<年代pan style="color:#A020F0">“EdgeColor”,<年代pan style="color:#A020F0">' w ');持有<年代pan style="color:#A020F0">在;Plot ([thr_rigrsure thr_rigrsure], [0 300],<年代pan style="color:#A020F0">“线宽”2);Plot ([thr_univthresh thr_univthresh], [0 300],<年代pan style="color:#A020F0">“r”,<年代pan style="color:#A020F0">“线宽”2);Plot ([thr_minimaxi thr_minimaxi], [0 300],<年代pan style="color:#A020F0">“k”,<年代pan style="color:#A020F0">“线宽”2);Plot ([-thr_rigrsure -thr_rigrsure], [0 300],<年代pan style="color:#A020F0">“线宽”2);Plot ([-thr_univthresh -thr_univthresh], [0 300],<年代pan style="color:#A020F0">“r”,<年代pan style="color:#A020F0">“线宽”2);Plot ([-thr_minimaxi -thr_minimaxi], [0 300],<年代pan style="color:#A020F0">“k”,<年代pan style="color:#A020F0">“线宽”2);

对于Stein的无偏风险估计(SURE)和极小极大阈值,大约保留了3%的系数。在通用阈值的情况下,所有值都将被拒绝。

我们知道细节系数向量是系数的叠加f系数e,即分解e导致细节系数,这是标准的高斯白噪声。

你使用后thselect要确定一个阈值,您可以为每个级别设置一个阈值。第二步可以使用wthcoef,直接处理原始信号的小波分解结构年代

软阈值或硬阈值

硬阈值和软阈值是<年代pan class="emphasis">收缩规则。确定阈值后,必须决定如何将该阈值应用于数据。

最简单的方案是<年代pan class="emphasis">阈值。让T表示阈值x您的数据。硬性阈值是

η x x | x | T 0 | x | < T

软阈值为

η x x T x > T 0 | x | T x + T x < T

可以使用硬规则或软规则应用阈值wthresh

y = linspace (1100);用力推= 0.4;ythard = wthresh (y,<年代pan style="color:#A020F0">“h”,用力推);ytsoft = wthresh (y,<年代pan style="color:#A020F0">“年代”,用力推);次要情节(131);情节(y);标题(<年代pan style="color:#A020F0">“原始数据”);次要情节(132);情节(ythard<年代pan style="color:#A020F0">‘*’);标题(<年代pan style="color:#A020F0">“硬阈值”);次要情节(133);情节(ytsoft<年代pan style="color:#A020F0">‘*’);标题(<年代pan style="color:#A020F0">“软阈值”);

处理无比例噪声和非白噪声

在实践中,通常不能直接使用基本模型。我们在这里检查可用的选项来处理主要去噪函数中的模型偏差wdenoise

最简单的用法wdenoise

sd = wdenoise (s)
哪个返回去噪版本sd原始信号的年代使用默认的参数设置,包括小波,去噪方法,软或硬阈值。任何默认设置都可以更改:
sd = wdenoise (s、n, DenoisingMethod, tptr,“小波”,wav,…ThresholdRule, sorh NoiseEstimate,宏大)
哪个返回去噪版本sd原始信号的年代获得使用tptr去噪方法。需要的其他参数有sorh公司拥有,wname.的参数sorh指定在级别上分解的细节系数的阈值n年代通过小波wav.其余的参数公司拥有是要指定的。它对应于估计数据中噪声方差的方法。

选项

噪声估计方法

“LevelIndependent”

“LevelIndependent”基于最小尺度(最高分辨率)小波系数估计噪声的方差。

“LevelDependent”

“LevelDependent”基于每个分辨率的小波系数估计噪声的方差。

对于更一般的程序,wdencmp函数执行小波系数阈值去噪和压缩目的,同时直接处理1-D和2-D数据。它允许您定义自己的阈值选择策略

xd = wdencmp(选择x, wav, n,用力推,sorh, keepapp);

在哪里

  • 选择= ' gbl ('用力推为统一阈值的正实数。

  • 选择= ' lvd '用力推是电平相关阈值的向量。

  • keepapp = 1保持近似系数,如前所述和

  • keepapp = 0允许近似系数阈值。

  • x信号要去噪和吗wavnsorh和上面一样。

小波去噪的作用

我们以Donoho和Johnstone的第一个例子开始一维去噪方法的例子。

块信号阈值

首先设置信噪比,并设置随机种子。

sqrt_snr = 4;init = 2055615866;

产生原始信号xref还有一个嘈杂的版本x加入标准高斯白噪声。同时绘制信号。

[xref x] = wnoise (sqrt_snr 1, 11日,init);次要情节(2,1,1)情节(xref)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">原始信号的)子图(2,1,2)图(x)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">噪声信号的)

图中包含2个轴对象。标题为“原始信号”的轴对象1包含一个类型为line的对象。标题为“噪声信号”的轴对象2包含一个类型为line的对象。

对小波分解得到的细节系数采用软启发式确定阈值对噪声信号去噪x使用sym8小波。使用默认设置wdenoise查看其余参数。与原始信号进行比较。

xd = wdenoise (x,<年代pan style="color:#A020F0">“小波”,<年代pan style="color:#A020F0">“sym8”,<年代pan style="color:#A020F0">“DenoisingMethod”,<年代pan style="color:#A020F0">“确定”,<年代pan style="color:#A020F0">“ThresholdRule”,<年代pan style="color:#A020F0">“软”);图subplot(2,1,1) plot(xref)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">原始信号的)子图(2,1,2)图(xd)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">的去噪信号)

图中包含2个轴对象。标题为“原始信号”的轴对象1包含一个类型为line的对象。标题为“去噪信号”的轴对象2包含一个类型为line的对象。

由于只有少量的大系数能表征原始信号,因此该方法具有很好的性能。

电信号去噪

当您怀疑非白噪声时,必须根据电平噪声的电平相关估计重新调整阈值。作为第二个例子,让我们在电信号的高扰动部分尝试这种方法。

首先加载电信号并从中选择一个片段。绘制部分。

负载<年代pan style="color:#A020F0">leleccumindx = 2000:3450;x = leleccum (indx);图绘制(indx x)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">原始信号的)

图中包含一个轴对象。标题为Original Signal的axis对象包含一个类型为line的对象。

去噪信号使用db4小波和三层小波分解和软固定形式阈值。为了处理非白噪声,使用电平相关噪声大小估计。与原始信号进行比较。

xd = wdenoise (x 3<年代pan style="color:#A020F0">“小波”,<年代pan style="color:#A020F0">“db4”,<年代pan style="color:#0000FF">...“DenoisingMethod”,<年代pan style="color:#A020F0">“UniversalThreshold”,<年代pan style="color:#0000FF">...“ThresholdRule”,<年代pan style="color:#A020F0">“软”,<年代pan style="color:#0000FF">...“NoiseEstimate”,<年代pan style="color:#A020F0">“LevelDependent”);图subplot(2,1,1) plot(indx,x)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">原始信号的)子plot(2,1,2) plot(indx,xd)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">的去噪信号)

图中包含2个轴对象。标题为“原始信号”的轴对象1包含一个类型为line的对象。标题为“去噪信号”的轴对象2包含一个类型为line的对象。

尽管在2410左右传感器故障开始前后的噪声性质存在时间不均一性,但结果相当好。

扩展到图像去噪

一维情况下的去噪方法同样适用于图像,也适用于几何图像。一维模型的直接平移是

年代 j f j + σ e j

在哪里<年代pan class="inlineequation"> e 是单位方差的高斯白噪声。

二维去噪过程采用相同的三个步骤,使用二维小波工具代替一维小波工具。对于阈值的选择,刺激(大小(s))被用来代替长度(年代)如果使用固定的表单阈值。

注意,除了“自动”的一维去噪情况,二维去噪和压缩使用wdencmp.为了说明二维去噪,加载一个图像并创建一个有噪声的版本。为了重现性,设置随机种子。

init = 2055615866;rng (init);负载<年代pan style="color:#A020F0">女人img = X;imgNoisy = img + 15*randn(size(img));

使用ddencmp求去噪值。在这种情况下,采用固定形式的阈值和电平噪声估计,阈值较软,并保留了近似系数。

(用力推,sorh keepapp] = ddencmp (<年代pan style="color:#A020F0">“窝”,<年代pan style="color:#A020F0">西弗吉尼亚州的, imgNoisy);用力推
用力推= 107.9838

用力推等于estimated_sigma *√日志(prod(大小(img))))

使用全局阈值选项去噪噪声图像。显示结果。

imgDenoised = wdencmp (<年代pan style="color:#A020F0">gbl(的imgNoisy,<年代pan style="color:#A020F0">“sym4”2用力推sorh keepapp);图colormap(粉色(255))sm = size(map,1);次要情节(2,2,1)图像(wcodemat (img, sm))标题(<年代pan style="color:#A020F0">原始图像的) subplot(2,2,2) image(wcodemat(imgNoisy,sm)) title(<年代pan style="color:#A020F0">“嘈杂的图像”) subplot(2,2,3) image(wcodemat(imgDenoised,sm)) title(<年代pan style="color:#A020F0">”“去噪图像)

图中包含3个轴对象。标题为“原始图像”的轴对象1包含一个类型为“图像”的对象。标题为“噪声图像”的轴对象2包含一个类型为图像的对象。标题为“去噪图像”的轴对象3包含一个类型为“图像”的对象。

去噪后的图像与原始图像对比良好。

一维小波方差自适应阈值

其思想是逐级定义时间相关阈值,然后提高去噪策略处理非平稳方差噪声模型的能力。

更准确地说,该模型假定(如前所述)观测到的信号等于叠加在噪声上的有趣信号

年代 n f n + σ e n

但噪声方差会随时间而变化。在几个时间间隔上有几个不同的方差值。值和间隔都是未知的。

让我们关注于估计变更点或等价的间隔的问题。使用的算法是基于Marc Lavielle关于使用动态规划检测变化点的原创著作(见[Lav99]参考文献).

让我们从一个固定设计的回归模型中生成一个信号,该模型有两个噪声方差变化点,分别位于位置200和600。为了重现性,设置随机种子。

init = 2055615866;rng (init);x = wnoise (10);bb = randn(1、长度(x));cp1 = 200;cp2 = 600;x = x + (bb (1: cp1), bb (cp1 + 1: cp2) / 4, bb (cp2 + 1:结束)];情节(x)标题(<年代pan style="color:#A020F0">噪声信号的)

图中包含一个轴对象。标题为“噪声信号”的axis对象包含一个类型为line的对象。

这个例子的目的是从信号中恢复两个变化点x

步骤1.通过抑制近似恢复噪声信号。首先执行单级小波分解使用db4小波。然后在第1级重建细节。

wname =<年代pan style="color:#A020F0">“db4”;列弗= 1;[c、l] = wavedec (x,列弗,wname);侦破= wrcoef (<年代pan style="color:#A020F0">' d 'c l wname, 1);图的阴谋(检波器)标题(<年代pan style="color:#A020F0">“1级细节”)

图中包含一个轴对象。标题为Level 1 Detail的axes对象包含一个类型为line的对象。

在这一阶段恢复的第一级的重建细节几乎是无信号的。如果信号的感兴趣部分有一个稀疏小波表示,它从变化点检测的观点捕获噪声的主要特征。为了去除几乎所有的信号,我们用均值替换最大的值。

步骤2。为了去除几乎所有的信号,用平均值替换2%的最大值。

x = (abs(检波器));v2p100 = x(修复(长度(x) * 0.98));印第安纳州=找到(abs(检波器)> v2p100);依据(印第安纳州)=意味着(依据);

步骤3。使用wvarchg函数用以下参数估计变化点:

  • 两个改变点之间的最小延迟为d = 10。

  • 最大更改点数为5。

[cp_est, kopt t_est] = wvarchg(侦破,5)
cp_est =<年代pan class="emphasis">1×2259 611
kopt = 2
t_est =<年代pan class="emphasis">6×61024 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<年代pan class="inlineequation"> 6,t_est(1张):包含了瞬间的方差变化点,和自kopt那么,提议的改变点数目是多少呢

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

您可以通过计算来替换估计的更改点:

cp_New = t_est(k+1,1:k)<年代pan style="color:#0000FF">结束
cp_New = 612
cp_New =<年代pan class="emphasis">1×2259 611
cp_New =<年代pan class="emphasis">1×3198 259 611
cp_New =<年代pan class="emphasis">1×4198 235 259 611
cp_New =<年代pan class="emphasis">1×5198 235 260 346 611

小波去噪分析测量

下面的测量和设置对于分析小波信号和图像是有用的:

  • M S E-均方误差(MSE)是数据与信号或图像近似之间的差除以元素数的平方范数。

  • 最大的错误-信号或图像近似中的最大绝对平方偏差。

  • L2-Norm比率-信号或图像近似的平方l2范数与输入信号或图像的比率。对于图像,图像在取l2范数之前被重塑为列向量

  • P s n r-峰值信噪比(PSNR),单位为分贝。PSNR仅对以每样本位或每像素位编码的数据有意义。

  • B P P-比特/像素比(BPP),这是压缩比(Comp. ratio)乘以8,假设每像素一个字节(8位)。

  • 薪酬比率—压缩比,即压缩图像中的元素个数除以原始图像中的元素个数,以百分比表示。