小波去噪与非参数函数估计- MATLAB & Simulink - MathWorks España金宝app

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

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

最通用的一维模型是

年代（n) =f（n) + σe（n）

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

在实践中,年代（n)通常是一个离散时间信号，具有相等的时间步长，被附加噪声损坏，您正在尝试恢复该信号。

更一般地，您可以查看年代（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 × 1随机向量替换为n × m随机矩阵，得到被加性噪声损坏的图像的恢复问题。

您可以使用以下代码获得该模型的1-D示例。

负载<年代pan style="color:#A020F0">cuspamax；Y = cuspamax+0.5*randn(size(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只有几个非零小波系数。这些函数具有稀疏小波表示。例如，一个平滑函数几乎无处不在，只有很少的突变，具有这样的性质。

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

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

分解
选择一个小波和一个电平N．计算信号的小波分解年代降至水平N．
阈值详细系数
对于从1到的每个级别N，对细节系数进行阈值。
重建
利用水平的原始近似系数计算小波重构N将层次细节系数从1修改为N．

去噪方法

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

THR = thselect(y,tptr)

它返回阈值。

选项	去噪方法
`“rigrsure”`	运用Stein无偏风险估计(SURE)原则进行选择
`“sqtwolog”`	固定形式(通用)阈值等于 $\sqrt{2 \ln （ N ）}$ 与N信号的长度。
`“heursure”`	混合使用前两个选项的选择
`“minimaxi”`	选择使用极大极小原理

选项“rigrsure”软性阈值估计器采用基于Stein无偏风险估计(二次损失函数)的阈值选择规则。你会得到一个特定阈值的风险估计值t．最大限度地降低t给出阈值的选择。
选项“sqtwolog”使用一个固定形式的阈值产生的最小最大性能乘以小因子成正比<年代pan class="emphasis">日志（<年代pan class="emphasis">长度（年代))。
选项“heursure”是前面两个选项的混合。因此，如果信噪比非常小，就会产生SURE估计是非常嘈杂的。因此，如果检测到这种情况，就使用固定的表单阈值。
选项“minimaxi”使用选定的固定阈值来yield对于理想程序的均方误差的最小最大性能。在统计学中使用极大极小原理来设计估计器。由于去噪后的信号可以被同化为未知回归函数的估计量，因此极小极大估计量是在给定函数集上实现最大均方误差的最小值的选项。

下面的例子展示了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(1e3,1);Thr_rigrsure = thselect(sig，<年代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);

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

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

使用后thselect为确定阈值，可以对属性的各个级别设置阈值。第二步可以使用wthcoef，直接对原始信号的小波分解结构进行处理年代．

软阈值或硬阈值

硬阈值和软阈值就是其中的例子<年代pan class="emphasis">收缩规则。在确定了阈值之后，必须决定如何将该阈值应用于数据。

最简单的方案是<年代pan class="emphasis">硬阈值。让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，<年代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，'Wavelet'，wav，…ThresholdRule, sorh NoiseEstimate,宏大)

哪个返回去噪后的版本sd原始信号的年代使用tptr去噪方法。所需的其他参数为sorh，公司拥有,wname．的参数sorh在级别上指定分解细节系数的阈值n的年代通过小波变换wav．剩下的参数公司拥有是要指定的。它对应于估计数据中噪声方差的方法。

选项	噪声估计方法
`“LevelIndependent”`	`“LevelIndependent”`基于最佳尺度(最高分辨率)小波系数估计噪声的方差。
`“LevelDependent”`	`“LevelDependent”`基于小波系数在每个分辨率水平上估计噪声的方差。

对于更一般的过程，请使用wdencmp函数在直接处理一维和二维数据的同时，对小波系数进行阈值处理以达到去噪和压缩的目的。它允许您定义自己的阈值策略选择

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

在哪里

Opt = 'gbl'而且用力推为均匀阈值的正实数。
Opt = 'lvd'而且用力推表示与级别相关的阈值的向量。
Keepapp = 1保持近似系数，如前面和
Keepapp = 0允许近似系数阈值。
x信号要去噪吗wav，n，sorh同上。

小波去噪的应用

打开实时脚本

我们开始1-D去噪方法的例子，第一个例子归功于Donoho和Johnstone。

信号阈值

首先设置信噪比(SNR)并设置随机种子。

Sqrt_snr = 4;Init = 2055615866;

生成原始信号xref还有一个有噪声的版本x加入标准高斯白噪声。把两个信号都画出来。

[xref,x] = wnoise(1,11,sqrt_snr,init);Subplot (2,1,1) plot(xref)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">原始信号的) subplot(2,1,2) plot(x)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">噪声信号的）

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

对小波分解得到的细节系数采用软启发式SURE阈值对噪声信号进行降噪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">原始信号的) subplot(2,1,2) plot(xd)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">的去噪信号）

图中包含2个轴对象。标题为Original Signal的Axes对象1包含一个类型为line的对象。标题为降噪信号的Axes对象2包含一个类型为line的对象。

由于只有少量的大系数表征原始信号，该方法性能非常好。

电信号去噪

当您怀疑是非白噪声时，必须通过对电平噪声的电平依赖估计来重新调整阈值。作为第二个例子，让我们在电信号的高度扰动部分尝试这种方法。

首先加载电信号并从中选择一段。绘制线段。

负载<年代pan style="color:#A020F0">leleccum指数= 2000:3450;X = leceleccum (indx);图示(indx,x)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">原始信号的）

图中包含一个轴对象。标题为Original Signal的axes对象包含一个类型为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">原始信号的) subplot(2,1,2) plot(indx,xd)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">的去噪信号）

图中包含2个轴对象。标题为Original Signal的Axes对象1包含一个类型为line的对象。标题为降噪信号的Axes对象2包含一个类型为line的对象。

结果是相当好的，尽管在2410左右传感器故障开始之前和之后的噪声性质的时间异质性。

图像去噪的扩展

打开实时脚本

对一维情况所描述的去噪方法也适用于图像，并且很好地适用于几何图像。1-D模型的直接翻译是

$年代（我， j ）＝ f （我， j ） + σ e （我， j ）$

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

二维去噪过程同样有三个步骤，并且使用二维小波工具而不是一维工具。对于阈值的选择，刺激(大小(s))用来代替长度(年代)如果使用固定的表单阈值。

注意，除了“自动”1-D去噪情况外，2-D去噪和压缩是使用wdencmp．为了说明2-D去噪，加载一个图像并创建一个有噪声的版本。为了重现性的目的，设置随机种子。

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

使用ddencmp找到去噪值。在这种情况下，固定形式的阈值与估计水平噪声，阈值是软的，并保持近似系数。

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

THR = 107.9838

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

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

imgDenoised = wdencmp(<年代pan style="color:#A020F0">gbl(的imgNoisy,<年代pan style="color:#A020F0">“sym4”2用力推sorh keepapp);图colormap(pink(255)) sm = size(map,1);Subplot (2,2,1) image(wcodemat(img,sm)) title(<年代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个轴对象。标题为Original Image的坐标轴对象1包含一个Image类型的对象。标题为“噪声图像”的坐标轴对象2包含一个图像类型的对象。标题为降噪图像的坐标轴对象3包含一个图像类型的对象。

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

一维小波方差自适应阈值

打开实时脚本

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

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

$年代（ n ）＝ f （ n ） + σ e （ n ）$

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

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

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

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)];情节(x)标题(<年代pan style="color:#A020F0">噪声信号的）

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

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

步骤1．通过抑制近似来恢复一个有噪声的信号。方法执行单级小波分解db4小波。然后在第1级重构细节。

wname =<年代pan style="color:#A020F0">“db4”；Lev = 1;[c,l] = wavedec(x,lev,wname);Det = wrcoef(<年代pan style="color:#A020F0">' d 'c l wname, 1);图示(det)标题(<年代pan style="color:#A020F0">“一级细节”）

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

在此阶段恢复的1级重建细节几乎是无信号的。如果信号中感兴趣的部分具有稀疏小波表示，则从变化点检测的角度捕获噪声的主要特征。为了去除几乎所有的信号，我们用平均值代替最大值。

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

X = sort(abs(det));V2p100 = x(fix(length(x)*0.98));Ind = find(abs(det)>v2p100);Det (ind) =平均值(Det);

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

两个更改点之间的最小延迟是d = 10。
变更点的最大数量是5个。

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

cp_est =<年代pan class="emphasis">1×2259 611

Kopt = 2

t_est =<年代pan class="emphasis">6×61024 0 0 0 0 0 612 1024 0 0 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"> $\leq$ 我 $\leq$ 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位)。
薪酬比率-压缩比，即压缩图像中的元素数量除以原始图像中以百分比表示的元素数量。