小波去噪和非参数函数估计- MATLAB和Simulink MathWorks西班牙万博1manbetx

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

小波工具箱™提供了一些函数的估计一个未知函数(信号或图像)的噪声。您可以使用这些函数来消除干扰信号,作为非参数函数估计的一种方法。

这是最一般的一维模型

年代(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×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只有少数几个非零小波系数。这些功能有稀疏的小波表示。举个例子,一个光滑函数几乎无处不在,只有一些突然的变化,有这样一个属性。

一般的小波去噪方法和非参数函数估计是将数据变换到小波域,阈值的小波系数和反变换。

你可以总结一下这些步骤为:

分解
选择小波和水平N。计算信号的小波分解年代到水平N。
阈值的细节系数
为每一个从1到水平N,细节系数阈值。
重建
计算小波重建使用原始的近似系数的水平N和修改后的细节系数从1到水平N。

去噪方法

小波工具箱支持多种去噪方法。万博1manbetx四种去噪方法的实现 thselect。每个方法对应于一个tptr选项在命令

用力推= thselect (y, tptr)

它返回阈值。

选项	去噪方法
`“rigrsure”`	选择使用斯坦的无偏风险估计原理(确定)
`“sqtwolog”`	固定形式(通用)阈值等于 $\sqrt{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(团体,<年代pan style="color:#A020F0">“sqtwolog”)thr_heursure = thselect(团体,<年代pan style="color:#A020F0">“heursure”)thr_minimaxi = thselect(团体,<年代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">在;情节([thr_rigrsure thr_rigrsure], [0 300),<年代pan style="color:#A020F0">“线宽”2);情节([thr_univthresh thr_univthresh], [0 300),<年代pan style="color:#A020F0">“r”,<年代pan style="color:#A020F0">“线宽”2);情节([thr_minimaxi thr_minimaxi], [0 300),<年代pan style="color:#A020F0">“k”,<年代pan style="color:#A020F0">“线宽”2);情节([-thr_rigrsure -thr_rigrsure], [0 300),<年代pan style="color:#A020F0">“线宽”2);情节([-thr_univthresh -thr_univthresh], [0 300),<年代pan style="color:#A020F0">“r”,<年代pan style="color:#A020F0">“线宽”2);情节([-thr_minimaxi -thr_minimaxi], [0 300),<年代pan style="color:#A020F0">“k”,<年代pan style="color:#A020F0">“线宽”2);

斯坦的无偏风险估计和极大极小阈值(确定),大约有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 (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函数执行小波系数阈值去噪和压缩的目的,而直接处理一维,二维数据。它允许您定义自己的阈值策略选择

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

在哪里

选择= ' gbl ('和用力推是一个正实数为统一的阈值。
选择= ' lvd '和用力推是一个向量的依赖阈值水平。
keepapp = 1保持近似系数,正如前面
keepapp = 0允许近似系数的阈值。
x信号去噪和吗wav,n,sorh上面是一样的。

小波去噪在行动

打开生活的脚本

我们开始第一个例子的一维去噪方法的示例归功于Donoho和约翰斯通。

块信号阈值

第一组一组信噪比(信噪比)和一个随机种子。

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标题原始信号包含一个类型的对象。坐标轴对象2标题嘈杂的信号包含一个类型的对象。

消除干扰噪声信号使用软启发式确定阈值对小波分解的细节系数获得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">“软”);图次要情节(2,1,1)情节(xref)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">原始信号的次要情节(2,1,2)情节(xd)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">的去噪信号)

图包含2轴对象。坐标轴对象1标题原始信号包含一个类型的对象。坐标轴对象2标题去噪信号包含一个类型的对象。

因为只有少量的大系数描述原始信号,该方法执行得很好。

电信号去噪

当你怀疑一个非白人噪音,阈值必须新等级相关估计的噪声水平。第二个例子,让我们试着高度摄动的方法一个电信号的一部分。

第一次加载电信号并选择一段。绘制部分。

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

图包含一个坐标轴对象。标题为原始信号的坐标轴对象包含一个类型的对象。

消除干扰信号使用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”);图次要情节(2,1,1)情节indx x轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">原始信号的次要情节(2,1,2)情节(indx xd)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">的去噪信号)

图包含2轴对象。坐标轴对象1标题原始信号包含一个类型的对象。坐标轴对象2标题去噪信号包含一个类型的对象。

结果是很好尽管时间异质性后噪声的性质和周围的传感器故障开始时间前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(大小(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 =大小(图1);次要情节(2,2,1)图像(wcodemat (img, sm))标题(<年代pan style="color:#A020F0">原始图像的次要情节(2,2,2)图像(wcodemat (imgNoisy, sm))标题(<年代pan style="color:#A020F0">“嘈杂的图像”次要情节(2,2,3)图像(wcodemat (imgDenoised, sm))标题(<年代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">噪声信号的)

图包含一个坐标轴对象。坐标轴对象与标题嘈杂的信号包含一个类型的对象。

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

步骤1。恢复的信号通过抑制一种近似。首先执行一个单层小波分解使用db4小波。然后在一级重建的细节。

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级细节”)

图包含一个坐标轴对象。坐标轴对象与标题1级细节包含一个类型的对象。

在1级恢复重建的细节信号在这个阶段几乎是免费的。它捕获的噪声变化的主要特征点检测的观点如果有趣的部分信号稀疏的小波表示。把几乎所有的信号,我们最大的值替换的意思。

步骤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 0 612 1024 0 0 0 0 259 611 1024 0 0 0 0 0 198 198 259 611 1024 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);

你可以取代估计改变点通过计算:

为k = 1:5 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-norm信号或图像近似输入信号或图像。对于图像,图像是重塑L2-norm前一个列向量
P S N R峰值信噪比(PSNR)在分贝。PSNR值仅为有意义的数据编码的比特每样或比特每像素。
B P P每像素比特率(BPP),这是压缩比(Comp。比)乘以8,假设每像素一个字节(8位)。
薪酬比率——压缩比,压缩图像的元素数量除以原始图像中的元素表示为一个百分比。