小波算法

Function wavelet(s,wname,n,options)

Begin

{

功能：

一维序列小波消噪。

参数：

s:一维序列

wname:小波函数名

现有小波函数名（小波函数的选取依靠经验）

Daubechies:

'db1' , 'db2', ... ,'db45' 'db1' 就是haar 小波函数

Coiflets :

'coif1', ... , 'coif5'

Symlets :

'sym2' , ... , 'sym8'

Biorthogonal:

'bior1.1', 'bior1.3' , 'bior1.5'

'bior2.2', 'bior2.4' , 'bior2.6', 'bior2.8'

'bior3.1', 'bior3.3' , 'bior3.5', 'bior3.7'

'bior3.9', 'bior4.4' , 'bior5.5', 'bior6.8'.

Reverse Biorthogonal:

'rbio1.1', 'rbio1.3' , 'rbio1.5'

'rbio2.2', 'rbio2.4' , 'rbio2.6', 'rbio2.8'

'rbio3.1', 'rbio3.3' , 'rbio3.5', 'rbio3.7'

'rbio3.9', 'rbio4.4' , 'rbio5.5', 'rbio6.8'.

n :分解层数

options ： 选项

选择字段说明

array('brief'颤孙:1, // 默认为1 采用简单剔除高频谐波 达到消噪的樱洞雹目的

// 如果为 0 采用估计序列噪音标准差剔除噪音，

'sigma':0, // 为0 默认采用 序列的高阶谐波估计标准差；也可自己输入值

'which':1, // 以 某一层谐波作为噪音估计的数据，默认第一层

'alpha':2, // 阈值惩罚系数，默认为2

"thr":0, // 阈值大小，默认0 采用谐波估计，也可以直接给出

'sorh':'s', // 阈值方式设置脊帆，'s' 软阈值，'h'硬阈值 默认为's'

)

返回结果：

一维数字数组，消噪后的序列。

范例：

s := array(2484.82690429688,2479.05493164063,2482.34301757813,2437.794921875,

2447.7548828125,2512.962890625,2443.05688476563,2433.15893554688,

2393.18310546875,2415.05395507813,2392.06201171875,2365.34301757813,

2359.21997070313,2344.787109375,2348.51611328125,2420.00,2438.7900390625,

2431.375,2440.40209960938,2383.48510742188,2377.51196289063,2331.36596679688,

2317.27490234375,2370.3330078125,2409.67211914063,2427.47998046875,

2435.61401367188,2473.40991210938,2468.25,2470.01904296875,2504.10791015625,

2508.09008789063,2528.2939453125,2509.79907226563,2503.8359375,2524.9189453125,

2479.53588867188,2481.083984375,2528.71411132813,2529.76098632813,2466.958984375,

2463.0458984375,2416.56201171875,2415.1298828125,2412.625,2395.06494140625,

2397.55395507813,2380.22412109375,2383.03393554688,2412.39306640625,

2333.4140625,2386.86010742188,2360.6640625,2333.22900390625,2325.90502929688,

2332.72998046875,2329.82006835938,2315.27001953125,2291.544921875,2248.59008789063,

2228.52490234375,2180.89501953125,2224.84008789063,2218.23510742188,2215.92993164063,

2191.14794921875,2186.29711914063,2204.78393554688,2190.11010742188,2166.205078125,

2170.01293945313,2173.56103515625,2199.4169921875,2169.38989257813,2148.45190429688,

2163.39501953125,2225.88989257813,2285.74389648438,2276.0458984375,2275.01000976563,

2244.580078125,2206.19311523438,2298.3759765625,2266.38403320313,2296.07495117188,

2319.11791992188,2285.0380859375,2292.61010742188,2268.080078125,2312.55590820313,

2330.40502929688,2331.13598632813,2291.90209960938,2347.53002929688,2349.58911132813,

2351.98095703125,2351.85498046875,2344.77099609375,2366.70190429688,2356.86010742188,

2357.18090820313,2363.59692382813,2381.42993164063,2403.5869140625,2409.55395507813,

2439.6279296875,2447.05688476563,2451.85693359375,2428.48706054688,2426.11499023438,

2460.69311523438)

n := 2

options := array('brief':1,'sigma':0,'which':1,'alpha':2,"thr":0,'sorh':'s')

return wavelet(s,wname,n,options)

天软数学组

20120627

}

if not ifarray(options) then options := array()

defaut := wavedefaut() union options

cout := 4

cl:=wavedec(s,n,wname) //小波分解

if defaut['brief']=1 then

ret :=wrcoef('a',cl[0],cl[1],wname,n)

else

begin

//***************小波消噪*************************************************

k := defaut['which'] //标准差估计选项 ，k 为 1 到 n的整数 默认为1；

if defaut['sigma']=0 then sigma := wnoisest(cl[0],cl[1],k)

else //通过小波第k层细节系数(谐波)估计 ，噪音标准差

sigma := defaut['segma']

if defaut['alpha']=0 then alpha :=2// alpha 惩罚因子 大于1 的数 一般为默认2；

else alpha := defaut['alpha']

if defaut['thr']=0 then

thr := wbmpen(cl[0],cl[1],sigma,alpha) //噪音信号全局阈值

else thr := defaut['thr']

sorh := defaut['sorh']

ret:=wdencmp('gbl',cl[0],cl[1],wname,n,thr,sorh)[0] //采用软阈值和近似信号进行消噪；

end //第一个参数为'gbl'为扩展接口备用，可以随意输入

return ret

end

function wavedefaut()

begin

return array('brief':1,'sigma':0,'which':1,'alpha':2,

"thr":0,'sorh':'s'

)

end

使用MATLAB设计小波变换程序中的若干问题3[转贴]使用MATLAB设计小波变换程序中的若干问题仍需探讨的问题：为什么使用PNG存储经小波变换后的重构图像变大？我曾在清华大学的多媒体课程的教师答疑中写了“老师：尊重事实:DB9阈值10的PNG文件就是比原文件大”和“续一：尊重事实:DB9阈值10的PNG文件就是比原文件大”，在林老师的鼓励和指导下，我进行了继续试验、分析，与刘赵璧（Anna）同学进行了探讨，并得到了Lily（姓名还不知道）同学的帮助，同时同学们也做了各自不同的实验，现在的实验结果可以说基本上比较明确，那就是有些图像就是会变大，这与图像的种类、纹理等密切相李银关。林老师曾经鼓励我去研究一下PNG的压缩方法，无奈我资质不够，至今在这方面的进展不大。由于临近期末考试，作业也要抓紧，所以我暂且将没有搞明白的内容搁置，待寒假期间再进行，希望对这些问题有各种看法也有兴趣研究的同学对此发表意见。以下是我最近试验、分析和阅读到的一些相关信息。试验结果我首先根据老师第三章的Haar矩阵算法推演出DB9的系数矩阵，并实现了分解重构及阈值处理程序，对几种照片进行了比较，然后使用3.1节的simplecmp进行了相同照片的实验，结果相当一致。细小差哪迅宴昌销别是因为我的程序对边界的扩展与MATLAB不一样，在设置阈值后引起了边界上小部分不一致造成的。表一：真彩色图像百合花的处理结果阈值PngHaar(Mat/Mine)0数Haat(Mat/Mine)PNGDb9(MAT/Mine)0数Db9(MAT/Mine)95973/95973095973/9597327524/2426895973/9597327/9574552/74292135838/136063101882/101992167412/1656621051976/51504163423/16374198411/98861199200/1957302032474/32346180167/18026792295/93660220629/217214从对比表中我们能够看到2个程序的结果相当一致，因此，我不再给出两种程序的对比，而是使用simplecmp直接处理的结果说明。将百合花图像使用[I,map]=rgb2ind(x,255)转换成为彩色图像处理，在将伪彩色图像转换为连续变换的灰度图像（如2.4常见问题中讨论的方法）进行处理：表二：百合花的伪彩色图像和处理后的灰度（gray）图像的处理结果阈值PngHaar(Index/Gray)0数Haar(Index/Gray)PNDb9(Index/Gray)0数Db9(Index/Gray)48535/43235048535/432356096/743048535/4323518/22553207/364509473/4362660362/499277009/528521058025/2360213362/5434464916/4781313202/658812060193/1434721948/6003966020/4601424468/73494其他伪彩色与进行加工的灰度图的结果与此完全一致，这也就说明了如果伪彩色文件的色板不是单调性递增就不适合小波分解。“Thecolorbartotherightoftheimageisnotsmoothanddoesnotmonotonicallyprogressfromdarktolight.Thistypeofindexedimageisnotsuitablefordirectwaveletdecompositionwiththetoolboxandneedstobepreprocessed.”。我对Facets进行同样的实验，结果与此一致。这种处理的结果可以从图像象素值的连续性来理解。这是处理与不处理的图像的中间一行的数据图。另外，不连续的图像质量在压缩后会被极大地破坏图2伪彩色文件变化前后的第128行数据的连续性情况对比分析多种试验图片基本能够反映类似的结果，虽然IndexedColorimage有时令Haar小波的分解重构图像出现增大现象，单经过处理之后，这种现象就会消失。然而对于DB9可以看到无论真彩色还是处理后的灰度图像都在阈值510处超过原始图像的大小，能不能因此得出DB9不适合进行图像压缩的结论呢？有一些同学确实这样认为，但我认为这种观点因为忽略了如何利用小波进行压缩和还原的过程，这也正是第四章老师为我们讲述的那些编码算法而造成的。在推荐材料[1]中也有类似的说明。图3、JPEG2000的基本结构看一下上图就可以明白为什么PNG不能衡量小波压缩的效率问题。上图的图像原始数据首先经过正变换（ForwardTransform）就是小波变换的得到小波系数，变换的小波系数经过阈值处理后进行量化，编码后得到压缩的图像文件JPEG2000，如果你没有JPEG2000的显示程序，那么你就不能看到它。它的显示程序就是由解码器从压缩数据中解出编码，进行反量化，得到小波系数，再实施逆变换（InverseTransform）就是小波系数重构。最终得到图像的原始数据。因此衡量小波变换的效率是应该看你选择的小波能不能分解出适合“编码器”压缩的小波系数，这种编码器不是PNG的LZ77，因为LZ77压缩小波分解系数的效率不是最好的。这种高效编码器在第四章可以找到。那么我们存储PNG文件的目的是什么呢？我认为压缩与去噪(de-noising)是同一种方法的两种提法。他们都使用了设置阈值的方法。我们可以仔细分析经过重构的PNG图片的质量来体会这种消除噪音的效果，也可以评定小波压缩后的图片的视觉质量，同时PNG的文件大小也可以让我们从LZ77算法的本质来理解小波变换压缩后的重构图像的内容变化情况。比如，我们可以从表2中的灰度图像在haar变换取阈值20时出现块状象素，文件大小变为14347，而db9却为46014，超过原始的PNG大小，但并不出现块状而是具有波状的特征。这本身说明了采用Harr小波压缩或去噪后重构的图像中相同的‘串’增多，便于PNG方式压缩，而db9则在相同阈值的情况下不会象Haar那样制造‘马赛克’，说明了它的平滑性，这也能帮助我们理解小波的特性。当然，当阈值继续增加后，超过某一界限，即使DB9也仍然会使PNG文件大小减小。这本身也就是由双尺度（Dyadic）小波变换的两种滤波器决定的。低通滤波结果相当于平均值，高通滤波结果相当于差值，差值能够保证重构图像的细节部分丢失最小，如果差值部分被阈值略去的过多，细节就会越来越少，平均意义的值约来越多，直到多到某一个临界值时（该图像的阈值取到40），重构的图像也可能出现较多的相同数字串，这就会提高PNG的压缩结果。下图是我对Haar（蓝色）小波取阈值为20，db9（红色）小波阈值取40时第128行1：32列的数据曲线与原始数据（黑色）曲线的对比。可以看出也db9在阈值=40时出现了较多的平均值，但比haar在阈值=20时的曲线要少的多。图4、haar（蓝色）和db9（红色）压缩后重构图像的第128行，1：32列的数据曲线不过，MATLAB给我们提供了量化的方法来决定如何选取阈值。在HELPWaveletToolbox:AdvanceConcepts:ChoosingtheOptimalDecomposition中提到了几种利用“熵”的概念来衡量如何选取合适的分解级。感兴趣的同学还可以参看wentropy,wdcbm2,wpdec的帮助。文献[1]中也提到了衡量压缩质量的客观化方法MSE,PNSR并指出小波的重构滤波器的长度越长，形状越规则越能够提供良好的压缩性能。上面对PNG的讨论因为没有足够的算法分析和程序解读，同时也没有准确的试验数据，因此只能作为猜测。但衡量小波压缩效率的方法我坚持认为不能以PNG文件大小来解说，如果采用图像文件大小来衡量，应该以JPEG2000来衡量。

您好，定义二维5/3小波提升方案的分解函数如下：

int multilevel_decomposition_2D( int * image, int width, int height, int * wavelet_data, int

levels )

定义二维5/3小波提升方案的重构函数如下：

int multilevel_reconstruction_2D ( int * image, int width, int height, int * wavelet_data, int levels )

类似地，对于9/7小波，我们用C++语言分别并蚂实现了二维的多层9/7小波提升方案的分解和重构函数如下：

void WaveletForward2D(int *data, int height, int width, int level)void WaveletInvert2D(int *data, int height, int width, int level)

值得注意的是，DSP平台上实现的5/3、9/7小波提升算法的输入和输出都是整型，这十分有利于直接困族对图像进行变换。另外，9/7小波提升算法中将对每次乘积的结果通过汪蔽弊floor函数进行取整，从而实现了整型到整型的无失真变换。DSP上只实现了CDF97小波。

在CCS 2 ('C6000)编译环境下编程，程序编译优化后，加载到DSP工作平台上。由于本算法使用的存储空间较多，程序中大量使用了堆空间，默认的存储分配方式无法满足要求，所以重新编写cmd文件分配堆栈空间大小如下：

-stack 0xf00 -heap 0x20000

---