主要内容

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

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

最常用的一维模型是

年代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到1N

去噪方法

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

thr=thselect(y,tptr)

它返回阈值。

选项

去噪方法

“rigrsure”

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

“sqtwolog”

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

2 自然对数 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(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);绘图([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 x | x | T 0 | x | < T

软阈值是

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

可以使用硬规则或软规则应用阈值威瑟什

y = linspace (1100);用力推= 0.4;ythard = wthresh (y,<年代pan style="color:#A020F0">“h”,thr);Ytsoft = wthresh(y,<年代pan style="color:#A020F0">“年代”第(131)子地块;地块(y);标题(<年代pan style="color:#A020F0">'原始数据'); 小批(132);阴谋,<年代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去噪方法。需要的其他参数包括索尔scalwname.参数索尔指定层次上分解的细节系数的阈值n属于年代被称为.剩下的参数scal是要指定的。它对应于估计数据中噪声方差的方法。

选项

噪声估计方法

“LevelIndependent”

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

'leveldeDependent'

'leveldeDependent'基于每个分辨率的小波系数估计噪声的方差。

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

xd = wdencmp(opt,x,wav,n,thr,sorh,legpapp);

在哪里

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

  • 选择= ' lvd '用力推是级别相关阈值的矢量。

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

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

  • x信号要去噪和吗n索尔和上面一样。

小波去噪的作用

我们开始使用第一个例子归功于Donoho和Johnstone的第一个例子的例子。

阻止信号阈值

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

sqrt_snr = 4;init = 2055615866;

生成原始信号xref还有一个嘈杂的版本x通过添加标准高斯白噪声。绘制两个信号。

[xref,x]=wnoise(1,11,sqrt_snr,init);子地块(2,1,1)打印(外部参照)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">原始信号的)子图(2,1,2)图(x)轴<年代pan style="color:#A020F0">紧标题(<年代pan style="color:#A020F0">噪声信号的)

使用从小波分解获得的细节系数上使用软启发式肯定阈值阈值来表示噪声信号x使用符号小波。使用默认设置wdenoise对于其余参数,请与原始信号进行比较。

xd = wddoise(x,<年代pan style="color:#A020F0">“小波”,<年代pan style="color:#A020F0">“sym8”,<年代pan style="color:#A020F0">“DenoisingMethod”,<年代pan style="color:#A020F0">“当然”,<年代pan style="color:#A020F0">'阈值',<年代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">的去噪信号)

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

电信号去噪

当您怀疑存在非白噪声时,必须通过对电平噪声的电平相关估计来重新调整阈值。作为第二个示例,让我们在电信号的高度扰动部分尝试该方法。

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

负载<年代pan style="color:#A020F0">勒勒库姆indx = 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">‘db3’,<年代pan style="color:#0000FF">......“DenoisingMethod”,<年代pan style="color:#A020F0">“大学三年级”,<年代pan style="color:#0000FF">......'阈值',<年代pan style="color:#A020F0">“软的”,<年代pan style="color:#0000FF">......“NoiseEstimate”,<年代pan style="color:#A020F0">'leveldeDependent');图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">的去噪信号)

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

扩展到图像去噪

1-DICACE描述的去噪方法也适用于图像并适用于几何图像。1-D模型的直接翻译是

年代 j f j + σ e j

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

二维去噪过程具有相同的三个步骤,使用二维小波工具代替一维工具。对于阈值选择,产品(尺寸)被用来代替长度如果使用固定的表单阈值。

注意,除了“自动”的一维去噪情况,二维去噪和压缩使用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">西弗吉尼亚州的,不认识);苏氨酸
thr = 107.9838.

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

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

imgDenoised = wdencmp (<年代pan style="color:#A020F0">“gbl”,imgnoisy,<年代pan style="color:#A020F0">'符号4'2用力推sorh keepapp);图colormap(粉色(255))sm = size(map,1);次要情节(2,2,1)图像(wcodemat (img, sm))标题(<年代pan style="color:#A020F0">'原始图像​​')子图(2,2,2)图像(WCODEMAT(IMGNoisy,SM))标题(<年代pan style="color:#A020F0">'吵闹的形象') subplot(2,2,3) image(wcodemat(imgDenoised,sm)) title(<年代pan style="color:#A020F0">'去世图像')

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

一维小波方差自适应阈值

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

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

年代 n f n + σ e n

但是噪声方差可以随时间变化。在几个时间间隔上有几个不同的方差值。这些值以及间隔都是未知的。

让我们专注于估计变化点或等效的时间间隔。所使用的算法基于Marc Lavielle的原始工作,关于使用动态编程检测变化点(见[LAV99]参考文献).

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

初始值=205565866;rng(初始值);x=wnoise(1,10);bb=randn(1,长度(x));cp1=200;cp2=600;x=x+[bb(1:cp1),bb(cp1+1:cp2)/4,bb(cp2+1:end)];绘图(x)标题(<年代pan style="color:#A020F0">噪声信号的)

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

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

wname =<年代pan style="color:#A020F0">‘db3’; lev=1;[c,l]=wavedec(x,lev,wname);det=wrcoef(<年代pan style="color:#A020F0">'D',c,l,wname,1);图情节(DET)标题(<年代pan style="color:#A020F0">“1级细节”)

在这个阶段恢复的1级的重建细节几乎是无信号。如果信号的有趣部分具有稀疏小波表示,它捕获来自改变点检测视点的噪声的主要特征。要删除几乎所有信号,我们通过平均值替换最大的值。

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

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

第3步。使用WVARCHG.函数用以下参数估计变化点:

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

  • 最大更改点数为5。

[cp_est, kopt t_est] = wvarchg(侦破,5)
美国东部=<年代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 259 611 1024 0 0 0 0 198 235 23124 611 1024 0 198 235 235 260 328 235 235 260 346 61235 235 260 324 235 235 260 324 235 235 260 326 235 24046 235 235 240 346 61235 231240 346 61312524 611 1024 0 198 235 260 324 235 235 240 324 235 235 260 326 235 246 246 244124

提出了两个变化点和三个间隔。由于噪声的三个区间方差是非常不同的,优化程序很容易检测到正确的结构。估计的变化点与真实的变化点很接近:200和600。

步骤4。(可选)替换估计的变更点。

2.<年代pan class="inlineequation"> ≤. ≤. 6,测试(i,1:i-1)包含瞬间的方差变化点,和自kopt.那么,提议的改变点数目是多少呢

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

您可以通过计算替换估计的变化点:

k=1:5 cp\u New=t\u 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-NOM比率-信号或图像近似值的平方L2范数与输入信号或图像的比值。对于图像,在采用L2范数之前,将图像重塑为列向量

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

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

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