2014-10-31 16:50:38 adong76 阅读数 3921


参考:

http://media.cs.tsinghua.edu.cn/~ahz/digitalimageprocess/chapter12/chapt12_ahz.htm


Matlab 小波变换

lean图像的行列应该满足2的幂次方

img  =  imread('lena.jpg');
img =rgb2gray(img);
img = double(img);

[ca1 ch1,cv1,cd1] = dwt2(img,'haar');
figure;
imshow([uint8(ca1),uint8(ch1);uint8(cv1),uint8(cd1)]);
title('first wavelet decomposition');

[ca2 ch2,cv2,cd2] = dwt2(ca1,'haar');
figure;
imshow([uint8(ca2),uint8(ch2);uint8(cv2),uint8(cd2)]);
title('second wavelet decomposition');

% restruct img
img_ca1 =  idwt2(ca2,ch2,cv2,cd2,'haar');
img_re = idwt2(img_ca1,ch1,cv1,cd1,'haar');
figure;
imshow( uint8(img_re));
title('wavelet restruct');

2017-10-24 16:15:56 HelloZEX 阅读数 9199

内容完全转载:

小波理论的基本概念及概述(第二版)

欢迎阅读此份关于小波变换的入门教程。小波变换是一个相对较新的概念(其出现大约是在20世纪80年代),但是有关于它的文章和书籍却不少。这其中大部分都是由数学专业人士写给其他同行看的,不过,仍然有大量数学专家不知道其他同行们讨论的是什么(我的一个数学教授就承认过)。换言之,大多数介绍小波变换的文献对那些小波新手们来说用处不大(此为个人观点)。

我刚开始接触小波变换的时候,曾经为了搞清楚小波变换这个这个神奇的世界到底发生了什么而苦苦挣扎,因为在这个领域的入门教材非常少。因而,我决定为新手们写一份教程。我自认为也是一个新手,必须承认,我也有很多理论细节没有弄清楚。不过,就工程应用而言,我认为弄清楚所有的理论细节大可不必。

这份教程将试着介绍一些小波理论的基本原理,并且不会给出这些原理和相关公式的证明,因为这份教程的目标读者暂时还不需要知道这些。不过,感兴趣的读者可以参阅引用的文献以便了解更深入的内容。

此篇文档假定你没有任何相关知识背景。要是有的话,请跳过以下内容,这些对你而言可能都是显然的。

要是你发现教程里有任何前后不协调或不正确的内容,请联系我。我很乐于收到关于教程的任何评论。

变换…啥?

首先,为什么需要变换,或者说到底什么是变换?

为了获取在原始信号中不易获得的信息,往往要对信号进行数学变换。以下篇幅均假定时域内信号为原始信号,经过数学变换后的信号为处理信号。

可用的变换有很多种,其中,傅立叶变换大概是目前最流行的。

实际中,多数信号的原始形式都是时域信号,也即不论如何测得的,信号总是关于时间的函数。换言之,绘制信号的图形时,一个轴代表时间(自变量),另一轴代表信号幅值(因变量)。在时域内作图,便可得到信号的时-幅表示。在多数信号处理有关的应用场景中,这种表示并不是最好的表示。很多时候,最易分辨的信息往往隐藏在信号的频率成分中。信号的频谱是指信号中的频率分量(或谱分量),其表示的是信号中存在哪些频率成分。

直觉上,我们都知道频率是跟事物的变化率有关的量。如果一样东西(专业术语应该为数学量或物理量)变化得很快,则它的频率就高;变换得慢,或者说变化得很平滑,则它的频率就低。如果该量保持不变,则其频率为零,或者说没有频率。例如,日报的频率就比月刊高(因为日报出版快)。

频率用“循环次数/秒”,或者用更常用的“赫兹”来衡量。例如,在美国,日常生活中所用交流电的频率是60Hz(世界上其他一些地区是50Hz)。这意味着,如果我们想要绘制电流变化曲线,得到的将是1秒内往复50次的正弦波。看下面几张图,第一幅图中是频率是3Hz的正弦信号,第二幅是频率10Hz的,第三幅则是频率50Hz,对比下吧。

那么怎样测量频率,或者说怎样得到一个信号中所含的频率成分呢?答案是傅立叶变换(FT)。对时域信号做傅立叶变换,就会得到信号的频谱。也就是说,此时我们绘制信号图形的话,一个轴是频率,另一个轴是频率分量的幅值。所得图像将告诉我们信号中包含的各种频率成分分别有多少。

频率轴从零开始,直至正无穷。每个频率都对应一个幅值。例如,如果我们对房间所用的电流信号做傅立叶变换,频谱图中在50Hz处会出现尖峰,其它频率对应的幅值则为零,因为信号中只包含了50Hz的频率分量。然而,很少有信号的傅立叶变换是如此简单的。实际中的信号大都包含多个频率分量。50Hz信号的傅立叶变换如下图所示:

图 1.4 50 Hz 信号的傅里叶变换

注意,图1.4给出了上下两张图,下图显示的其实是上图的前半部分。这是因为实值信号的频谱图是左右对称的,这点暂时不理解也无妨。上图能够看出这一特性。不过,由于后一半对称部分只不过是前一半图形的镜像,并未提供额外信息,因此,这部分经常不画出来。下文中出现的多数频谱图,我将只绘出前半部分。

为什么需要频率信息?

通常,一些在时域中不易看出的信息很容易在频域中观察到。

看一个生物信号的例子。设想我们正在观察一个心电信号。心脏专家一般都熟知典型的健康人心电图的形状。与这些典型形状存在显著偏差往往是疾病的征兆。

一些病征在时域表示的心电信号中并不明显。过去,心脏专家一般用记录在磁带上的时域心电图来分析心电信号。最近,新型的数字心电记录仪/分析仪可以利用心电图的频域信息来判断病征是否存在。对心电信号的频率成分进行分析能使他们更容易的诊断病情。

上面只是一个说明频率成分作用的简单例子。当前,傅立叶变换已经被用于不同的领域,涵盖了工程领域的各个分支。

尽管傅立叶变换可能是使用最多的(特别在电气工程领域),但它并不是唯一的变换。许多其他的变换也常为工程师和数学家们所用,如希尔伯特变换、短时傅立叶变换(下文会有更多介绍)、魏格纳分布和雷登变换,当然还有教程的主角——小波变换,而这些也仅是工程师和数学家们所用变换中的一小部分。每种变换都有其应用领域,也有其优缺点,小波变换也不例外。

为了更好地理解为什么需要小波变换,我们需要更深刻地认识傅立叶变换。傅立叶变换是一种可逆变换,即它允许原始信号和处理信号之间互相变换。但是,在任意时刻只有一种信号形式是可用的。也就是说,在时域信号中不包含频率信息,而经过傅里叶变换后的信号则不包含任何时间信息。说到这,头脑里很自然地会提出一个问题,为什么需要同时知道时间和频率信息呢?

我们马上就会知道,答案是具体问题具体分析。回想一下,傅立叶变换给出了信号中的频率信息,即它可以告诉我们原始信号包含各个频率成分到底有多少,但是并未告诉我们某个频率信号何时出现。对于所谓的平稳信号,这些信息并不需要。

让我们进一步探讨一下平稳的概念,因为它在信号分析中具有重要意义。如果信号中的频率分量不随时间变化,则称这类信号为平稳信号。平稳信号中的频率分量一直保持不变,那么,自然无需知道频率分量是何时出现的,因为所有的频率分量出现在信号的每一刻!!!

以如下信号为例:

\( x(t)=cos(2picdot10t)+cos(2picdot25t)+cos(2picdot50t)+cos(2picdot100t) \)

这是个平稳信号,因为任何时刻都包含10,25,50和100Hz的频率。信号的图形如下:

图 1.5

下图为它的傅立叶变换:

图 1.6

图1.6中的上图是图1.5中信号的频谱图,下图为上图的放大,给出了我们关注部分的频率范围。注意四个频率10,25,50和100Hz的频谱分量。

与图1.5中的信号不同,下图所示的就是一个非平稳信号。图1.7中,信号的频率随着时间一直在变化,这种信号称为线性调频信号,是一种非平稳信号。

图 1.7

让我们再看一个例子,图1.8绘出的是一个包含四个频率分量的信号,它们分别在不同时刻出现,因此这是一个非平稳信号。0至300ms时是100Hz的正弦波,300-600ms时则是50Hz的正弦波,600-800ms时是25Hz的正弦波,最后的200ms内是10Hz正弦波。

图 1.8

下图是它的傅立叶变换:

图 1.9

不要介怀图中的那些小波纹,这是由信号中频率突变引起的,在这里并不重要。注意,高频分量的幅值比低频分量大,这是因为高频信号(300ms)比低频信号(200ms)持续时间更长。(频率分量幅值的精确值并不重要)。

除了那些波纹,图中的一切看起来都是正确的。频谱图有四个尖峰,对应原始信号中的四个频率分量,幅值也差不多是合理的…没错

错!

当然了,也不全错,但也不全对。对图1.5中的信号,考虑如下问题:各个频率分量都是在什么时刻出现的?

答案是

在所有时刻!还记得平稳信号吗,所有频率分量在信号的整个持续时间内一直存在。10Hz的频率分量一直存在,50Hz的分量也是,100Hz的分量依然是。

现在,让我们来考虑一下图1.7或1.8中的非平稳信号。

各个频率分量都是在什么时刻出现的?

对于图1.8中的信号,我们知道,第一个时间区间内出现的是频率最高的分量,最后一个时间区间内出现的是频率最低的分量。图1.7中,信号的频率成分随时间连续变化,因此,对这些信号来说,各个频率分量并未在所有时刻一直存在。

现在,对比图1.6和1.9,两幅频谱图的相似之处是显而易见的。两幅图中都包含了四个相同的频率分量,即10,25,50和100Hz。除了一些小波纹和两幅图中各频率分量的幅值(这些幅值可以做归一化处理)有所区别,两幅频谱图几乎是相同的,尽管两个信号在时域内差别很大。两个信号都包含了相同的频率分量,但是前者中,各频率分量存在于信号的整个周期内,而后者的频率分量则分别存在于不同的时间区间内出现。那么,为什么两个完全不同的信号,频谱图形这么相像呢?回想一下,傅立叶变换仅仅给出了信号的频谱分量,但却没有给出任何关于这些分量出现时间的信息。因此,傅立叶变换并不适用于分析非平稳信号,但有一个例外:

如果我们仅关心信号中包含哪些频率分量而不关心它们出现的时间,傅立叶变换仍可用于处理非平稳信号。但是,如果我们想知道频率分量出现的确切时间(区间),傅立叶变换就不再适用了。

实际应用中,由于平稳的和非平稳的信号都很多,很那将二者区分开来。例如,几乎所有的生物信号都是非平稳的,包括广为人知的心电图(ECG)、脑电图(EEG)和肌电图(EMG)。

再次注意,傅立叶变换仅能给出信号中包含哪些频率分量,仅此而已。

当需要对频谱分量进行时间定位时,我们就需要一个可以给出信号时-频表示的变换。

终极解决方案:小波变换

小波变换是这种类型的变换,它提供了信号的时频表示(还有一些变换也可给出这些信息,如短时傅立叶变化,魏格纳分布等等)。

特定的频谱分量在特定的时刻出现往往具有特殊的意义。这些情况下,了解这些特定的频谱分量出现的时间区间会非常有用。例如,在脑电图中,事件相关电位的延迟时间需要特别注意(事件相关电位是指大脑对某一特定刺激的反应,类似闪光灯,延迟时间是从接受刺激到作出反应之间耗费的时间)。

小波变换能够同时提供时间和频率信息,因此给出了信号的一种时频表示。

小波变换到底是如何奏效的完全是另外一个故事,需要在理解了短时傅立叶变换(STFT)之后再做解释。小波变换的出现是为了改进短时傅立叶变换(STFT)。STFT将在教程的第II部分详细阐述。现在暂时可以认为小波变换是为了解决STFT中遇到的有关分辨率的问题而发展起来的。

为了长话短说,我们略过时域信号处理中有关于各种高通和低通滤波器的相关内容。这些滤波器用来过滤信号中的低频和高频部分分量。这类方法被重复实施,每次都会从信号中滤除一些频率分量。

这里解释一下滤波是如何奏效的:设想我们有一个信号,其中频率最高的分量为1000Hz。第一步,我们通过高通和低通滤波器把信号分成两个信号(滤波器必须满足某些特定的条件,即容许条件),结果得到了同一信号的两个部分,0-500Hz的部分(低通部分)和500-1000Hz的部分(高通部分)。

然后,我们可以拿其中一部分(通常是低通部分)或者二部分,然后对每一部分继续进行相同的操作。这个过程叫做分解。

假设我们拿低频部分做了处理,现在我们就有了3组数据,分别为信号在0-250Hz,250-500Hz和500-1000Hz的部分。

然后再对低通部分的信号继续做高通和低通滤波处理;现在我们就有了4组数据,分别为0-125Hz,125-250Hz,250-500Hz和500-1000Hz。我们持续进行这个过程,直到将信号分解到一个预先定义的水平。这样我们就有了一系列信号,这些信号实际上都来自相同的信号,但是每一个都对应不同的频带。我们知道每个信号对应的频段,如果我们将这些信号放在一起画出三维图,一个轴表示时间,频率在另外一个轴上,幅度在第三个轴上。这幅图会告诉我们各个频率出现哪些时刻(这里有一个问题,叫做“不确定性原理”,即我们不能精确地知道哪个频率出现在哪些时间点,我们仅能知道某一频段出现在哪一时间区间内,后文中将有更多介绍)。

不过,我仍想简单地解释一下:

不确定性原理最早由海森堡发现并阐述,其表述为:移动粒子的动量和位置不可同时确定。在我们这个课题里则是这样:

时-频平面内的一个确定的点上,信号的频率和时间信息不能同时知道。换句话说:在任一时刻,我们无法确定存在哪个频谱分量。我们最多只能做到,在一个给定的时间区间内存在哪些频谱分量。这是一个分辨率的问题,也是研究者们从快速傅立叶变换(STFT)切换到小波变换(WT)的主要原因。快速傅立叶变换的分辨率随时间是固定不变的,而小波变换则能给出可变的分辨率:

高频信号在时域内很好分辨,低频信号则在频域内容易分辨。这意味着,相对于低频分量,高频分量更容易在时域内定位(有更小的相对误差)。反而言之,低频分量更容易在频域内定位。看下面的网格图:

对上图的解释是:最上面一行表明,高频信号有更多的采样点和较短的采样间隔。就是说,高频信号更容易在时域内分辨。最下面一行是对低频信号的采样,描述信号的特征点较少,因此,低频信号在时域内并不容易分辨。

在离散时间的情形中,信号的时间分辨率与先前相同,但是现在,频率信息的分辨率在每一个阶段都不同。注意到,低频信号更容易在频域内分辨,高频则不然。注意,相邻频率分量的间隔是如何随频率增高而增大的。

下面是连续小波变换的例子:

我们构造一个正弦信号,具有两个频率成分,分别处在两个不同的时间区间:

注意低频分量先出现,然后是高频分量。

图 1.10

图 1.11

注意,上图中代表频率的轴被标记为了“尺度”。“尺度”的概念将会在后续章节进行阐述,但这里需要注意的是,尺度是频率的倒数,即尺度越大频率越低,尺度越小频率越高。因此,图中的小的峰值对应的是信号中的高频分量,大的峰值对应的是信号中的低频分量(在时域内,低频分量先于高频分量出现)。

你可能被图中的频率分辨率搞晕了,因为高频信号似乎也有很好的频率分辨率。但请注意,高频(低尺度)信号处分辨率较好的是尺度分辨率,而非频率分辨率。尺度分辨率高意味着频率分辨率低,反之亦然。更多相关内容将在后续部分介绍。

未完,待续…

以上是此份教程的第一部分,我试着给出信号处理的简要概述——傅里叶变换和小波变换


转载自:http://blog.jobbole.com/101976/


其他参考资料:

http://blog.csdn.net/alihouzi/article/category/3132373

2014-10-26 21:53:38 EbowTang 阅读数 14952

本文代码的实现严重依赖前面的两篇文章:


1,一维信号的小波阈值去噪


2,小波变换一维Mallat算法的C++实现


注本文的大部分文字提取于参考论文
重要说明:本文代码是学习小波变换时写的,文中的代码有较严重的性能问题(但是运算结果是正确的),如你需要本代码,需要自行优化或者更改(一维阈值去噪那篇文章中的性能挺快的)

一,小波阈值去噪基本理论

      图像在获取或传输过程中会因各种噪声的干扰使质量下降,这将对后续图像的处理产生不利影响.所以必须对图像进行去噪处理,而去噪所要达到的目的就是在较好去除噪声的基础上,良好的保持图像的边缘等重要细节.在图像去噪领域得到广泛的应用.本博文根据小波的分解与重构原理,实现了基于硬阈值和软阈值函数的小波阈值去噪的C++版本,最终结果与matlab库函数运算结果完全一致。

1,小波阈值处理

小波阈值收缩法是Donoho和Johnstone提出的,一下便是养活不少学者的三篇基础论文:

【1】 Donoho D L. De-noising by soft-thresholding. IEEE Trans- actions on Information Theory, 1995, 41(3): 613−627
【2】 Donoho D L, Johnstone I M. Adapting to unknown smooth- ness via wavelet shrinkage. Journal of the American Statistic
Association, 1995, 90(432): 1200−1224
【3】 Donoho D L, Johnstone I M, Kerkyacharian G, Picard D. Wavelet shrinkage: asymptopia? Journal of Royal Statisti-
cal Society Series B, 1995, 57(2): 301−369

小波阈值去噪其主要理论依据是,小波变换具有很强的去数据相关性,它能够使信号的能量在小波域集中在一些大的小波系数中;而噪声的能量却分布于整个小波域内.因此,经小波分解后,信号的小波系数幅值要大于噪声的系数幅值.可以认为,幅值比较大的小波系数一般以信号为主,而幅值比较小的系数在很大程度上是噪声.于是,采用阈值的办法可以把信号系数保留,而使大部分噪声系数减小至零.小波阈值收缩法去噪的具体处理过程为:将含噪信号在各尺度上进行小波分解,设定一个阈值,幅值低于该阈值的小波系数置为0,高于该阈值的小波系数或者完全保留,或者做相应的“收缩(shrinkage)”处理.最后将处理后获得的小波系数用逆小波变换进行重构,得到去噪后的图像.


2,阈值函数的选取

  阈值去噪中,阈值函数体现了对超过和低于阈值的小波系数不同处理策略,是阈值去噪中关键的一步。设w表示小波系数,T为给定阈值,sign(*)为符号函数,常见的阈值函数有:


硬阈值函数:     (小波系数的绝对值低于阈值的置零,高于的保留不变)

     

软阈值函数:   (小波系数的绝对值低于阈值的置零,高于的系数shrinkage处理)

      


值得注意的是:

1) 硬阈值函数在阈值点是不连续的,在下图中已经用黑线标出。不连续会带来振铃,伪吉布斯效应等。

2) 软阈值函数,原系数和分解得到的小波系数总存在着恒定的偏差,这将影响重构的精度使得重构图像的边缘模糊等现象.

所以这里也添加一种简单的改进阈值函数,我们称之为半阈值

注意:图片中的公式有误,应该是sgn(w)*(|w|-P*T)(即将软函数中的阈值T缩小):

         

三种阈值处理策略见下图:



其实不少文章出现各种优秀的改进方案(于是有养活了不少学者,非本文重点):


3,阈值的确定

选取的阈值最好刚好大于噪声的最大水平,可以证明的是噪声的最大限度以非常高的概率低于,此阈值是Donoho和Johnstone提出的,其实我一直很想吐槽为什么阈值和信号的size相关呢?当然我的疑问也是大家的疑问,此问题有养活了一批学者,其中根号右边的这个参数(叫做sigma)就是估计出来的噪声标准偏差(第一级分解出的小波细节系数,即整个HH系数的绝对值的中值),本文将用此阈值去处理各尺度上的细节系数。


4,阈值策略

阈值策略有两种,一种是全局阈值策略,一种是分层阈值策略,从读论文了解到,全局阈值策略有他的缺陷:如果采用全局阈值进行处理,则会对所有尺度下的高频系数进行同一阈值处理,然而随着小波分解尺度的增加,有用信息分解后的小波系数将会越来越大,而噪声系数却会越来越小,所以为了防止高尺度下有用信息的分解系数被过度处理,分层阈值处理将会更合理。但是从我实际测试的结果来看------没什么区别!!!简便起见代码还是采用全局阈值。


为了不误导小伙伴们,请注意:这张图中阈值的估算不是来自CV1而是来自CD1,因为CD1是两次高通滤波的结果,里面含有最多的噪声,估计出来的噪声偏差也是最准确的。





5,小波阈值去噪过程


   有用信号经小波变换后, 其能量将集中在少数的小波系数上, 而噪声点的小波系数互不相关, 分布在各个尺度的所有时间轴上. 保留小波变换的各尺度下的模极大值点, 而将其他点置零或最大程度的减小, 然后将处理后的小波系数做小波逆变换, 即可达到抑制噪声的目的.阈值去噪是通过对变换域系数与阈值进行比较判断, 然后将处理后的系数进行逆变换重构去噪图像. 小波阈值去噪法的具体步骤如下:

步骤 1. 图像的小波分解: 确定小波函数和分解层次 N, 对图像进行 N 层的小波分解;
步骤 2. 阈值处理: 对分解得到的各层系数选择阈值, 并对细节系数进行阈值判断;
步骤 3. 图像重构: 对阈值处理后的系数通过小波逆变换重建图像.
   信号和噪声在小波域内具有不同的相关性. 信号在尺度间相应位置上的小波系数具有很强的相关性, 而噪声的小波系数则具有弱相关性或者不相关. 在阈值去噪中, 由于所选定的阈值通常固定, 不会随着小波系数的不同而变化, 这就不可避免地会对部分小波系数进行误判,于是又养活了一批学者。。。分析改进,分析改进,发文章!


1),小波的分解





2,小波的重构





二,Matlab库函数实现

1,核心库函数说明

1)wavedec2

图像的多级小波分解,将返回分解出来的小波系数以及小波系数的各级长度

2)waverec2

多级小波系数的重构,重构出原信号

3)wthresh函数

对系数进行指定类型(全局阈值或者分层阈值)的阈值去噪处理,软硬阈值处理函数。下图所示程序和运行结果可以比较清晰地看出该程序的执行过程。

clear;
y = linspace(-1,1,100);
thr = 0.4;
ythard = wthresh(y,'h',thr);
ytsoft = wthresh(y,'s',thr);
subplot(311);plot(y);
subplot(312);plot(ythard);
subplot(313);plot(ytsoft);



更具体的函数说明可以在,matlab里键入“doc 函数名”将得到很详细的说明,当然也可以百度

2,软,硬阈值处理效果:



哎哟,不错哦,半阈值去噪效果在视觉上的确有明显的改进效果




局部放大图像:

四幅图象均放大两倍,便于查看区别

这幅图像是源图像放大的效果:









3,完整的代码实现

说明:代码实现对图像一层分解后的细节系数进行全局阈值处理(即LHL,LH,HH均用同一阈值处理)。

% 获取输入参数  
w    = 'db3';%小波类型  
n    = 3;%分解层数   

sorh1 = 'hard';%硬阈值  
sorh2 = 'soft';%软阈值  
sorh3 = 'half';%半阈值   
% 对图像进行小波分解  
[c,l] = wavedec2(octimage,n,w);  
  

%求取阈值  
N = numel(octimage);  
[chd1,cvd1,cdd1] = detcoef2('all',c,l,1);   
cdd1=cdd1(:)';  
sigma = median(abs(cdd1))/0.6745;%提取细节系数求中值并除以0.6745  
thr = sigma*sqrt(2*log(N))/sqrt(1+sqrt(n)); %对阈值做了改进


% 对小波系数全局阈值处理  
cxchard = c;% 保留近似系数  
cxcsoft = c;% 保留近似系数  
cxchalf = c;% 保留近似系数  
justdet = prod(l(1,:))+1:length(c);%截取细节系数(不处理近似系数)  
  
% 阈值处理细节系数  
cxchard(justdet) = myWthresh(cxchard(justdet),sorh1,thr);%硬阈值去噪  
cxcsoft(justdet) = myWthresh(cxcsoft(justdet),sorh2,thr);%软阈值去噪  
cxchalf(justdet) = myWthresh(cxchalf(justdet),sorh3,thr);%软阈值去噪  

%小波重建  
xchard = waverec2(cxchard,l,w);  
xcsoft = waverec2(cxcsoft,l,w);  
xchalf = waverec2(cxchalf,l,w);   

figure(2);   
imshow(uint8(xchard(1:iCount, :)));title('硬阈值去噪图像');    
    
figure(3);   
imshow(uint8(xcsoft(1:iCount, :)));title('软阈值去噪图像');   

figure(4);   
imshow(uint8(xchalf(1:iCount, :)));title('半阈值去噪图像'); 




三,C加加实现

说明:如同一维的阈值去噪一样,在执行自己编写的wavedec2函数时必须先初始化,初始化的目的是为了获取信号的长度,选择的是什么小波,以及分解的等级等信息,然后计算出未来的各种信息,比如每个等级的系数的size,其中共有变量m_msgCL2D记录了这些信息。二维小波分解的初始化函数如下:

//初始化二维图像的分解信息,保存未来需要的信息
bool  CWavelet::InitDecInfo2D(
	const int height,//预分解的图像的高度
	const int width,//预分解的图像的宽度
	const int Scale,//分解尺度
	const int dbn//db滤波器编号,有默认值
	)
{
	if (dbn != 3)
		SetFilter(dbn);

	if (height < m_dbFilter.filterLen - 1 || width < m_dbFilter.filterLen - 1)
	{
		cerr << "错误信息:滤波器长度大于信号的高度或者宽度!" << endl;
		return false;
	}

	int srcHeight = height;
	int srcWidth = width;
	m_msgCL2D.dbn = dbn;
	m_msgCL2D.Scale = Scale;
	m_msgCL2D.msgHeight.resize(Scale + 2);
	m_msgCL2D.msgWidth.resize(Scale + 2);
	//源图像的尺寸
	m_msgCL2D.msgHeight[0] = height;
	m_msgCL2D.msgWidth[0] = width;

	//每一尺度上的尺寸
	for (int i = 1; i <= Scale; i++)//注意:每个尺度的四个分量的长宽是一样的
	{
		int exHeight = (srcHeight + m_dbFilter.filterLen - 1) / 2;//对称拓延后系数的长度的一半
		srcHeight = exHeight;
		m_msgCL2D.msgHeight[i] = srcHeight;

		int exWidth = (srcWidth + m_dbFilter.filterLen - 1) / 2;//对称拓延后系数的长度一半
		srcWidth = exWidth;
		m_msgCL2D.msgWidth[i] = srcWidth;
	}
	m_msgCL2D.msgHeight[Scale + 1] = srcHeight;
	m_msgCL2D.msgWidth[Scale + 1] = srcWidth;

	//计算总的数据个数
	int tmpAllSize = 0;
	int curPartSize = 0;
	int prePartSize = 0;
	for (int i = 1; i <= Scale; i++)
	{
		curPartSize = m_msgCL2D.msgHeight[i] * m_msgCL2D.msgWidth[i];
		tmpAllSize += curPartSize * 4 - prePartSize;
		prePartSize = curPartSize;
	}
	m_msgCL2D.allSize = tmpAllSize;
	m_bInitFlag2D = true;
	return true;
}


1,核心函数的实现

1)二维信号的单次分解

说明:本函数建立在一维的小波分解函数基础上(DWT)

// 二维数据的小波分解
void  CWavelet::DWT2(
	double *pSrcImage,//源图像数据(存储成一维数据,行优先存储)
	int height,//图像的高
	int width,//图像的宽
	double *pDstCeof//分解出来的图像系数
	)
{

	if (!m_bInitFlag2D)
	{
		cerr << "错误信息:未初始化,无法对信号进行分解!" << endl;
		return;
	}

	if (pSrcImage == NULL || pDstCeof == NULL)
	{
		cerr << "错误信息:dwt2数据无内存" << endl;
		Sleep(3000);
		exit(1);
	}
	
	int exwidth = (width + m_dbFilter.filterLen - 1) / 2 * 2;//pImagCeof的宽度
	int exheight = (height + m_dbFilter.filterLen - 1) / 2 * 2;//pImagCeof的高度

	double *tempImage = new double[exwidth*height];

	
	// 对每一行进行行变换
	double *tempAhang = new double[width]; 
	double *tempExhang = new double[exwidth]; // 临时存放每一行的处理数据
	for (int i = 0; i < height; i++)
	{
		for (int j = 0; j < width; j++)
			tempAhang[j] = pSrcImage[i*width + j];//提取每一行的数据

		DWT(tempAhang, width, tempExhang);

		for (int j = 0; j < exwidth; j++)
			tempImage[i*exwidth + j] = tempExhang[j];
	}

	// 对每一列进行列变换
	double *tempAlie = new double[height]; // 临时存放每一列的转置数据
	double *tempexlie = new double[exheight]; // 临时存放每一列的处理数据
	for (int i = 0; i < exwidth; i++)
	{
		// 列转置
		for (int j = 0; j < height; j++)
			tempAlie[j] = tempImage[j*exwidth + i];//提取每一列数据

		//执行变换
		DWT(tempAlie, height, tempexlie);

		// 反转置
		for (int j = 0; j < exheight; j++)
			pDstCeof[j*exwidth + i] = tempexlie[j];
	}

	AdjustData(pDstCeof, exheight, exwidth);//调整数据

	delete[] tempAlie;
	tempAlie = NULL;
	delete[] tempexlie;
	tempexlie = NULL;

	delete[] tempAhang;
	tempAhang = NULL;
	delete[] tempExhang;
	tempExhang = NULL;

	delete[] tempImage;
	tempImage = NULL;
	
}



2)二维信号的单次重构

说明:

//二维小波反变换
void  CWavelet::IDWT2(
	double *pSrcCeof, //二维源图像系数数据
	int dstHeight,//重构出来后数据的高度
	int dstWidth,//重构出来后数据的宽度
	double *pDstImage//重构出来的图像
	)
{
	int srcHeight = (dstHeight + m_dbFilter.filterLen - 1) / 2 * 2;
	int srcWidth = (dstWidth + m_dbFilter.filterLen - 1) / 2 * 2;//pSrcCeof的高度
	IAdjustData(pSrcCeof, srcHeight, srcWidth);//调整成LL,HL,LH,HH

	double *tempAline = new double[srcHeight]; // 临时存放每一列的数据
	double *tempdstline = new double[dstHeight]; // 临时存放每一列的重构结果

	double *pTmpImage = new double[srcWidth*dstHeight];
	// 列重构
	for (int i = 0; i < srcWidth; i++)//每一列
	{
		// 列转置
		for (int j = 0; j<srcHeight; j++)
			tempAline[j] = pSrcCeof[j*srcWidth + i];//提取每一列

		IDWT(tempAline, dstHeight, tempdstline);

		// 反转置
		for (int j = 0; j < dstHeight; j++)
			pTmpImage[j*srcWidth + i] = tempdstline[j];
	}


	// 对每一行进行行变换
	double *tempAhang = new double[srcWidth];
	double *tempdsthang = new double[dstWidth]; // 临时存放每一行的处理数据
	for (int i = 0; i < dstHeight; i++)
	{
		for (int j = 0; j < srcWidth; j++)
			tempAhang[j] = pTmpImage[i*srcWidth + j];//提取每一行的数据

		IDWT(tempAhang, dstWidth, tempdsthang);

		for (int j = 0; j < dstWidth; j++)
			pDstImage[i*dstWidth + j] = tempdsthang[j];
	}


	delete[] tempAline;
	tempAline = NULL;
	delete[] tempdstline;
	tempdstline = NULL;

	delete[] tempAhang;
	tempAhang = NULL;
	delete[] tempdsthang;
	tempdsthang = NULL;

	delete[] pTmpImage;
	pTmpImage = NULL;
}


3)二维信号的多级分解

说明:对于每一级分解都将调用单次二维分解函数来实现,所以本函数是建立在函数IDW2基础上

// 二维小波多级分解,需要先初始化获取未来数据信息
bool CWavelet::WaveDec2(
	double *pSrcData,//源图像数据,存储为一维信号
	double *pDstCeof//分解后的系数,它的大小必须是m_msgCL2D.allSize
	)
{
	if (!m_bInitFlag2D)
	{
		cerr << "错误信息:未初始化,无法对图像进行分解!" << endl;
		return false;
	}
	if (pSrcData == NULL || pDstCeof == NULL)//错误:无内存
		return false;

	int height = m_msgCL2D.msgHeight[0];
	int width = m_msgCL2D.msgWidth[0];
	int scale = m_msgCL2D.Scale;

	// 临时变量,图像数据
	double *tempImage = new double[height*width];
	int maxCoefSize =4 * m_msgCL2D.msgHeight[1] * m_msgCL2D.msgWidth[1];
	double *tempDst = new double[maxCoefSize];

	for (int i = 0; i < height*width; i++)
		tempImage[i] = pSrcData[i];

	int gap = m_msgCL2D.allSize - maxCoefSize;
	for (int i = 1; i <= scale; i++)
	{
		DWT2(tempImage, height, width, tempDst);

		// 低频子图像的高和宽
		height = m_msgCL2D.msgHeight[i];
		width = m_msgCL2D.msgWidth[i];
		
		for (int j = 0; j < height*width; j++)
			tempImage[j] = tempDst[j];//提取低频系数(近似系数)
		//
		for (int j = 0, k = gap; j < 4 * height*width; j++, k++)
			pDstCeof[k] = tempDst[j];//所有系数
		gap -= 4 * m_msgCL2D.msgWidth[i + 1] * m_msgCL2D.msgHeight[i + 1] - height*width;
	}
	delete[] tempDst;
	tempDst = NULL;
	delete[] tempImage;
	tempImage = NULL;
	return true;
}


4)多级分解系数的重构

// 根据多级分解系数重构出二维信号,必须先初始化获取分解信息
bool CWavelet::WaveRec2(
	double *pSrcCoef,//多级分解出的源系数
	double *pDstData//重构出来的信号
	)
{
	if (!m_bInitFlag2D)
	{
		cerr << "错误信息:未初始化,无法对信号进行分解!" << endl;
		return false;
	}
	if (pSrcCoef == NULL || pDstData == NULL)//错误:无内存
		return false;

	int height = m_msgCL2D.msgHeight[0];
	int width = m_msgCL2D.msgWidth[0];
	int decLevel = m_msgCL2D.Scale;

	int maxCeofSize = 4 * m_msgCL2D.msgHeight[1] * m_msgCL2D.msgWidth[1];

	double *pTmpImage = new double[maxCeofSize];

	int minCeofSize = 4 * m_msgCL2D.msgHeight[decLevel] * m_msgCL2D.msgWidth[decLevel];
	for (int i = 0; i < minCeofSize; i++)
		pTmpImage[i] = pSrcCoef[i];

	int gap = minCeofSize;
	for (int i = decLevel; i >= 1; i--)
	{
		int nextheight = m_msgCL2D.msgHeight[i - 1];//重构出来的高度
		int nextwidth = m_msgCL2D.msgWidth[i - 1];//重构出来的宽度
		IDWT2(pTmpImage, nextheight, nextwidth, pDstData);

		if (i > 1)//i==1已经重构出来了,不再需要提取系数
		{
			for (int j = 0; j < nextheight*nextwidth; j++)
				pTmpImage[j] = pDstData[j];
			for (int j = 0; j < 3 * nextheight*nextwidth; j++)
				pTmpImage[nextheight*nextwidth + j] = pSrcCoef[gap + j];
			gap += 3 * nextheight*nextwidth;
		}
	}
	
	delete[] pTmpImage;
	pTmpImage = NULL;

	return true;
}


2,函数正确性验证

1),二维单次分解与重构测试



2)二维多级分解与重构测试

说明:对二维数据进行了5层分解,选取的是小波族db3



3,阈值去噪结果:

说明:本测试只是模拟测试,对图像的处理也是一样的(完全一致)

硬阈值去噪结果




软阈值去噪结果





实际的VC++图像处理结果为:

源噪声图像为:



注意以下是采用:db6,3层分解,软阈值去噪,阈值是在前文提及的阈值基础上缩小2.5倍得到的效果:



注意以下是采用:db6,3层分解,硬阈值去噪,阈值是在前文提及的阈值基础上缩小2.5倍得到的效果:




附带上述matlab验证程序

clc;  
clear all;  
close all;  
% %%%%%%%%%%%%%%%%%%%%%%%%通过matlab的函数来实现阈值去噪%%%%%%%%%%%%%%%%%%%%%%%%%% %  
X=[ 10, 12, 30, 4, 5, 61, 2, 3;
    41, 5, 6, 27, 3, 4, 15, 6;
    72, 8, 41, 5, 6, 7, 8, 9;
    5, 64, 7, 8, 9, 14, 6, 27;
    8, 9, 40, 31,10, 12, 30, 4;
    50, 61, 2, 3, 41, 5, 6, 27];

X=double(X);  

% 获取输入参数  
wname    = 'db3';%小波类型  
n    = 3;%分解层数  
sorh1 = 'h';%硬阈值  
sorh2 = 's';%软阈值  

% 对图像进行小波分解  
[c,l] = wavedec2(X,n,wname);

%求取阈值
N = numel(X);
[chd1,cvd1,cdd1] = detcoef2('all',c,l,1); 
cvd1=cvd1(:)';
sigma = median(abs(cvd1))/0.6745;%提取细节系数求中值并除以0.6745
thr = sigma*sqrt(2*log(N)); 

% 对小波系数全局阈值处理  
cxch = c;% 保留近似系数  
cxcs = c;% 保留近似系数  
justdet = prod(l(1,:))+1:length(c);%截取细节系数(不处理近似系数)  

% 阈值处理细节系数  
cxch(justdet) = wthresh(cxch(justdet),sorh1,thr);%硬阈值去噪  
cxcs(justdet) = wthresh(cxcs(justdet),sorh2,thr);%软阈值去噪  

%小波重建  
xch = waverec2(cxch,l,wname);  
xcs = waverec2(cxcs,l,wname);  


鉴于不少人要代码,出于不误导人,给一些我自己的建议:
1,对程序要有自己的思考,要有自己的验证和学习过程,写时当时我还是个C++新手!点击下载完整程序
2,二维小波变换不建议用,因为速度太慢了(自己徒手写的),一维可以用。
3,推荐你另外的小波库网站,http://wavelet2d.sourceforge.net/

注:本博文为EbowTang原创,后续可能继续更新本文。如果转载,请务必复制本条信息!

原文地址:http://blog.csdn.net/ebowtang/article/details/40481539

原作者博客:http://blog.csdn.net/ebowtang


参考资源:

【1】《数字图像处理》(冈萨雷斯matlab第二版)
【2】 Donoho D L. De-noising by soft-thresholding. IEEE Trans- actions on Information Theory, 1995, 41(3): 613−627
【3】 Donoho D L, Johnstone I M. Adapting to unknown smooth- ness via wavelet shrinkage. Journal of the American Statistic
Association, 1995, 90(432): 1200−1224
【4】 Donoho D L, Johnstone I M, Kerkyacharian G, Picard D. Wavelet shrinkage: asymptopia? Journal of Royal Statisti-
cal Society Series B, 1995, 57(2): 301−369
【5】杨恢先,王绪四,改进阈值与尺度间相关的小波红外图像去噪
【6】《小波分析及其应用》,孙延奎著
【7】杨建国.小波分析及其工程应用[M].北京:机械工业出版社.2005
【8】毛艳辉.小波去噪在语音识别预处理中的应用.上海交通大学硕士学位论文.2010
【9】matlab各种函数说明,及其内部函数实现
【10】http://www.bearcave.com/misl/misl_tech/wavelets/haar.html

2019-10-13 21:00:10 weixin_44225182 阅读数 197

小波变换实现图像压缩

代码

X=imread('a5.jpg');
X=rgb2gray(X);
subplot(221); imshow(X);
title('原始图像');
%对图像用小波进行层小波分解
[c,s]=wavedec2(X,2,'haar');
%提取小波分解结构中的一层的低频系数和高频系数
cal=appcoef2(c,s,'haar',1);
ch1=detcoef2('h',c,s,1);      %水平方向
cv1=detcoef2('v',c,s,1);      %垂直方向
cd1=detcoef2('d',c,s,1);      %斜线方向
%各频率成份重构
a1=wrcoef2('a',c,s,'haar',1);
h1=wrcoef2('h',c,s,'haar',1);
v1=wrcoef2('v',c,s,'haar',1);
d1=wrcoef2('d',c,s,'haar',1);
c1=[a1,h1;v1,d1];
subplot(222),imshow(c1,[]);
title ('分解后低频和高频信息');
%进行图像压缩
%保留小波分解第一层低频信息
%首先对第一层信息进行量化编码
ca1=appcoef2(c,s,'haar',1);
ca1=wcodemat(ca1,440,'mat',0);
%改变图像高度并显示
ca1=0.5*ca1;
subplot(223);imshow(cal,[]);
title('第一次压缩图像');
%保留小波分解第二层低频信息进行压缩
ca2=appcoef2(c,s,'haar',2);
%首先对第二层信息进行量化编码
ca2=wcodemat(ca2,440,'mat',0);
%改变图像高度并显示
ca2=0.25*ca2;
subplot(224);imshow(ca2,[]);
title('第二次压缩图像');  

结果
在这里插入图片描述
在这里插入图片描述

2016-11-07 10:56:25 shenziheng1 阅读数 3861

1.前言:

图像去噪是信号处理的一个经典问题,传统的去噪方法多采用平均或线性法法进行,最常用到的就是维纳滤波,但是他的降噪效果并不是很明显。小波分析法开辟了非线性降噪的先河,小波能够降噪得益于小波变换的以下特点:低熵性(小波系数稀松分布,使图像变换后的熵降低)、多分辨率特性(极好的刻画了信号的非平稳性)、去相关性(噪声在变换后有白化趋势,小波域更有利于去噪)

目前,主流的小波去噪方法主要集中在三个方面:基于小波变换模极大值降噪、基于相邻尺度小波系数相关性去燥、基于小波变换域阈值去噪(硬阈值与软阈值、全局阈值与局部自适应阈值)。

2.小波图像去燥实现的步骤:

1.二维信号的小波分解。选择一个小波和小波分解的层次N,然后计算信号s到第N层的分解。

2.对高频系数进行阈值量化,对于从1~N的每一层,选择一个阈值,并对这一层的高频系数进行软阈值量化处理。

3.二维小波重构

3.小波系数阈值降噪

<span style="font-size:18px;">clear all;
load facets;
subplot(221);image(X);
colormap(map);
xlabel('(a)原始图像');
axis square
%产生含噪声图像
init=2055615866;randn('seed',init)
x=X+50*randn(size(X));
subplot(222);image(x);
colormap(map);
xlabel('(b)含噪声图像');
axis square
%下面进行图像的去噪处理
%用小波画数coif3对x进行2层小波分解
[c,s]=wavedec2(x,2,'coif3');
%提取小波分解中第一层的低频图像,即实现了低通滤波去噪
%设置尺度向量n
n=[1,2];
%设置阈值向量p
p=[10.12,23.28];
%对三个方向高频系数进行阈值处理
nc=wthcoef2('h',c,s,n,p,'s');
nc=wthcoef2('v',c,s,n,p,'s');
nc=wthcoef2('d',c,s,n,p,'s');
%对新的小波分解结构[nc,s]进行重构
x1=waverec2(nc,s,'coif3');
subplot(223);image(x1);
colormap(map);
xlabel('(c)第一次去噪后的图像');
axis square;
xx=wthcoef2('v',nc,s,n,p,'s');
x2=waverec2(xx,s,'coif2');%图像的二维小波重构
subplot(2,2,4);image(x2);
colormap(map);
xlabel('(d)第二次消噪后图解');
axis square;
</span>
降噪结果:


4.全局软阈值降噪

<span style="font-size:18px;"><span style="font-size:18px;">clear all;
load detfingr;
subplot(131);image(X);
colormap(map);
xlabel('(a)原始图像');
axis square;
init=255615866;
randn('state',init);  %添加随机值
x=X+20*randn(size(X));
subplot(132);image(x);
colormap(map);
xlabel('(b)含噪图像');
axis square;
[thr,sorh,kep]=ddencmp('den','wv',x);   %使用全局阈值选项进行图像消噪处理
xd=wdencmp('gbl',x,'sym5',2,thr,sorh,kep);
subplot(133);image(xd)
colormap(map);
xlabel('(c)消噪图像');
axis square;</span><span style="font-size:24px;">
</span></span>
降噪结果:


图像的小波变换

阅读数 11509

图像Haar小波变换

阅读数 3060

图像处理-小波变换

阅读数 4379

小波变换图像处理

阅读数 4721

没有更多推荐了,返回首页