小波变换Mallat算法c++计算实例
小波分析实验：二维离散小波变换(Mallat快速算法)
小波分析实验：实验2二维离散小波变换（Mallat快速算法）实验目的：在理解离散小波变换原理和Mallat快速算法的基础上，通过编程对图像进行二维离散小波变换，从而加深对二维小波分解和重构的理性和感性认识，并能提高编程能力，为今后的学习和工作奠定基础。
实验工具：计算机，matlab6.5分解算法:重构算法: “"二工必（刃- 2上*［十三g （刃- 2k ）d ［ *分解算法写成矩阵的形式！ （lb g 的长度为4）4[0]如]力⑵ h[3] 0 0 0 '［勺【0】• 记"h[0] h[\]h[2]山⑶ …• ••••・ • •C J=勺【1] • •申[2] h[3] 00 0-.^[0] ^[1]_.勺［乃-1】_>[0] g[l] g ⑵ g[3] 0 • • •e=• 0 •g[0] g[l]g ⑵ • • g[3]■ • •・■ 0• D J =<[i]■•目2] ■g[3]0 0…茎0] 畀］|g[0] g[l] g[2] g[3] 0 0 0 I0 0 g[0] g[l]g[2] S [3] - 0• ••••• • ••・•・・■ • • g[2] g[3] 0 00 ...g[0] g[l]J |_勺4-1[叨］I二・(2»于是Mallat分解公式为矩阵变换？丄Cj- = PC^................. ⑶卩D j = Q D J-L..... .......... ⑷重构算法写成矩阵变换：-C J_I =C$ + Dj------------------------------------ (5) 4M NPPq. 一片『峰值信噪比计算公式：P沁沁逻竺皿E卢H耿V 屈E M ｛皿，00分别表示原始图像和重建图像，且本实验采取的一些小技乐P (I)分SW法…编程时用如下思想：（h, g 的长度为4）“今［1］勺［刀-1］■ V■■丐⑼£［1］ 4刀-1】将数据。
小波学习之一（单层一维离散小波变换DWT的Mallat算法C++和MATLAB实现）
⼩波学习之⼀（单层⼀维离散⼩波变换DWT的Mallat算法C++和MATLAB实现）1 Mallat算法离散序列的Mallat算法分解公式如下：其中，H(n)、G(n)分别表⽰所选取的⼩波函数对应的低通和⾼通滤波器的抽头系数序列。
从Mallat算法的分解原理可知，分解后的序列就是原序列与滤波器序列的卷积再进⾏隔点抽取⽽来。
离散序列的Mallat算法重构公式如下：其中，h(n)、g(n)分别表⽰所选取的⼩波函数对应的低通和⾼通滤波器的抽头系数序列。
2 ⼩波变换实现过程(C/C++)2.1 ⼩波变换结果序列长度⼩波的Mallat算法分解后的序列长度由原序列长SoureLen和滤波器长FilterLen决定。
从Mallat算法的分解原理可知，分解后的序列就是原序列与滤波器序列的卷积再进⾏隔点抽取⽽来。
即分解抽取的结果长度为(SoureLen+FilterLen-1)/2。
2.2 获取滤波器组对于⼀些通⽤的⼩波函数，简单起见，可以通过Matlab的wfilters(‘wavename’)获取4个滤波器；特殊的⼩波函数需要⾃⾏构造获得。
下⾯以db1⼩波函数(Haar⼩波)为例，其变换与重构滤波器组的结果如下：//matlab输⼊获取命令>> [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db1')//获取的结果Lo_D =0.7071 0.7071Hi_D =-0.7071 0.7071Lo_R =0.7071 0.7071Hi_R =0.7071 -0.70712.3 信号边界延拓在Mallat算法中，假定输⼊序列是⽆限长的，⽽实际应⽤中输⼊的信号是有限的采样序列，这就会出现信号边界处理问题。
对于边界信号的延拓⼀般有3种⽅法，即零延拓、对称延拓和周期延拓。
3种延拓⽅法⽐较情况如下：对于正交⼩波变换来说，前两种延拓⽅法实现起来⽐较简单，但重建时会产⽣边界效应，⽽且分解的层数越多，产⽣的边界效应越显著。
一维小波变换的C++实现
⼀维⼩波变换的C++实现 将⼩波展开系数当成离散信号，尺度函数和⼩波函数的MRA⽅程系数看成数字滤波器组，根据Mallat快速算法的原理，⼩波变换对数据的处理⽅法可简化成对信号逐级采样和滤波的过程。
图1 ⼩波变换的滤波器实现(a)分解算法 (b)重构算法 ⼀层⼩波分解算法流程如图2所⽰，信号将先经过⼩波分解低通滤波器和⾼通滤波器，随后被降采样，实现数据重构。
⽽滤波算法可简化为待处理信号与滤波器数组卷积的过程，为了保证卷积前和卷积后数组的长度相同，结合⼩波变换中数组延拓的思想，在实际编程过程中，可以将超过信号长度的那段数据以前端对齐的⽅式与前⾯⼀段数据相加。
将卷积后的数组每2个点采样⼀次，即可获得⼩波分解后的尺度系数和⼩波系数。
图2 ⼀层⼩波分解算法（X：待分解数组；H，G：⼩波分解滤波器；C，D：⼩波重构后数组） ⼩波重构算法是⼩波分解算法的逆运算，其流程为升采样和滤波，最后数据相加实现重构。
⼩波重构算法中滤波可视为系数与⼩波重构滤波器的卷积，与⼩波正变换类似，在重构算法中，需要将卷积后的数组末位对齐相加，获得与原数组长度相同的卷积结果。
将⼩波系数和尺度系数以2为步长进⾏升采样，将获得的新数组分别经过⼩波重构低通滤波器和⾼通滤波器，再将滤波后的两组数据相加，即实现了⼀层⼩波重构。
1#define LENGTH 5122#define LEVEL 43#define L_core 645static void Covlution(double data[], double core[], double cov[], int LEN)6 {7double temp[LENGTH + L_core - 1] = {0};8int i = 0;9int j = 0;1011for(i = 0; i < LEN; i++)12 {13for(j = 0; j < L_core; j++)14 {15 temp[i + j] += data[i] * core[j];19for(i = 0; i < LEN; i++)20 {21if(i < L_core - 1)22 cov[i] = temp[i] + temp[LEN + i];23else24 cov[i] = temp[i];25 }2627 }2829static void Covlution2(double data[], double core[], double cov[], int LEN)30 {31double temp[LENGTH + L_core - 1] = {0};32int i = 0;33int j = 0;3435for(i = 0; i < LEN; i++)36 {37for(j = 0; j < L_core; j++)38 {39 temp[i + j] += data[i] * core[j];40 }41 }4243for(i = 0; i < LEN; i++)44 {45if(i < L_core - 1)46 cov[i + LEN - L_core + 1] = temp[i] + temp[LEN + i];47else48 cov[i - L_core + 1] = temp[i];49 }5051 }5253static void DWT1D(double input[], double output[], double LF[], double HF[], int l)54 {55int i = 0;56double temp[LENGTH] = {0};57int LEN = LENGTH / pow(2, l - 1);5859 Covlution(input, LF, temp, LEN);60for(i = 1; i < LEN; i += 2)61 {62 output[i/2] = temp[i];63 }6465 Covlution(input, HF, temp, LEN);66for(i = 1; i < LEN; i += 2)67 {68 output[LEN/2 + i/2] = temp[i];69 }70 }7172static void DWT(double input[], double output[], double LF[], double HF[], int len[])73 {74int i;75int j;7677 len[0] = len[1] = LENGTH / pow(2, LEVEL);78for(i = 2; i <= LEVEL; i++) len[i] = len[i - 1] * 2;7980 DWT1D(input, output, LF, HF, 1);81for(i = 2; i <= LEVEL; i++)82 {83for(j = 0; j < len[LEVEL + 2 - i]; j++) input[j] = output[j];84 DWT1D(input, output, LF, HF, i);85 }86 }8788static void IDWT1D(double input[], double output[], double LF[], double HF[], int l, int flag) 89 {90int i = 0;91double temp[LENGTH] = {0};92int LEN = l * 2;9394if(flag) Covlution2(input, HF, temp, LEN);95else Covlution2(input, LF, temp, LEN);9697for(i = 0; i < LEN; i++)98 {99 output[i] = temp[i];103static void IDWT(double input[], double output[], double LF[], double HF[], int len[], int level)104 {105int i;106int j;107for(j = 0; j < len[LEVEL + 1 - level]; j++)108 {109 output[2 * j] = 0;110 output[2 * j + 1] = input[j];111 }112for(j = 0; j < 2 * len[LEVEL + 1 - level]; j++)113 {114 input[j] = output[j];115 }116 IDWT1D(input, output, LF, HF, len[LEVEL + 1 - level], 1);117118for(i = level - 1; i > 0; i--)119 {120for(j = 0; j < len[LEVEL + 1 - i]; j++)121 {122 input[2 * j] = 0;123 input[2 * j + 1] = output[j];124 }125 IDWT1D(input, output, LF, HF, len[LEVEL + 1 - i], 0);126 }127 }⽤C++算法实现的⼩波变换结果与MATLAB实现的⼩波变换结果对⽐（⼼电信号，db5⼩波，5层分解）。
小波变换(内附奇异值分析matlab程序)
2、算法及其应用实例
小波在信号的奇异性检测中的应用举例 信号的突变点和奇异点等不规则部分通常包含重要信息，一般信号 的奇异性分为两种情况： （1）信号在某一时刻其幅值发生突变，引起信号的非连续，这种类 型的突变称为第一类型的间断点； （2）信号在外观上很光滑，幅值没有发生突变，但是信号的一阶微 分有突变发生且一阶微分不连续，这种类型的突变称为第二类型的间 断点。 应用小波分析可以检测出信号中的突变点的位置、类型以及变 化的幅度。
程序代码
load nearbrk; x=nearbrk; %使用db4对信号进行2层分解 [c,l]=wavedec(x,2,‘db4’); subplot(411); subplot(4,1,i+2); plot(x); plot(d); ylabel('x'); ylabel(['d',num2str(3-i)]); %对分解的第六层低频系数进行重构 end a=wrcoef('a',c,l,'db4',2); subplot(412); plot(a); ylabel('a2'); for i=1:2 %对分解的第2层到第1层的高频系数 进行重构 a=wrcoef('a',c,l,'db4',3-i);
3、小波分析的优缺点
小波变换与Fourier变换相比，是一个时间和频域的局域变换因而能 有效地从信号中提取信息，通过伸缩和平移等运算功能对函数或信号进 行多尺度细化分析（Multiscale Analysis），解决了Fourier变换不能 解决的许多困难问题。 小波变换存在以下几个优点： 小波变换存在以下几个优点： (1)小波分解可以覆盖整个频域(提供了一个数学上完备的描述) (2)小波变换通过选取合适的滤波器，可以极大的减小或去除所提取得不 同特征之间的相关性 (3)小波变换具有“变焦”特性，在低频段可用高频率分辨率和低时间分 辨率(宽分析窗口)，在高频段，可用低频率分辨率和高时间分辨率(窄分 析窗口) 。 (4)小波变换实现上有快速算法(Mallat小波分解算法)。
(整理)小波变换课件第4章小波变换的实现技术.
第4章 小波变换的实现技术4.1 Mallat 算法双正交小波变换的Mallat 算法：设{}n h h =、{}n g g =、{}n h h =、{}n g g =为实系数双正交小波滤波器。
h ，g 是小波分析滤波器，h ，g 是小波综合滤波器。
h 表示h 的逆序，即n n h h -=。
若输入信号为n a ，它的低频部分和高频部分以此为1n a -和1n d -，小波分解与重构的卷积算法：11()()n n n na D a h d D a g --⎧⎪=*⎨=*⎪⎩ n11()()n n a Uah Ud g --=*+*先进行输入信号和分析滤波器的巻积，再隔点采样，以形成低频和高频信号。
对于有限的数据量，经过多次小波变化后数据量大减，因此需对输入数据进行处理。
4.1.1 边界延拓方法下面给出几种经验方法。
1. 补零延拓是假定边界以外的信号全部为零，这种延拓方式的缺点是，如果输入信号在边界点的值与零相差很大，则零延拓意味着在边界处加入了高频成分，造成很大误差。
实际应用中很少采用。
0121,0,,,,...,,0,0,......n s s s s -2.简单周期延拓将信号看作一个周期信号,即k n k s s +=。
简单周期延拓后的信号变为这种延拓方式的不足之处在于，当信号两端边界值相差很大时，延拓后的信号将存在周期性的突变，也就是说简单周期延拓可在边界引入大量高频成分，从而产生较大误差。
3. 周期对称延拓这种方法是将原信号在边界上作对称折叠，一般分二1）当与之做卷积的滤波器为奇数时，周期延拓信号为2）当与之做卷积的滤波器为偶数时，周期延拓信号为4. 光滑常数延拓在原信号两端添加与端点数据相同的常数。
0121,,,...,,n s s s s -0121,,,...,,n s s s s -0121,,,...,,n s s s s -0,...s 1,...,n s -01221,,,...,,,n n s s s s s --0121,,,...,,n s s s s -21012,...,,,,,...n s s s s s -321212,,,...,,,,...n n n s s s s s s ---10012,,...,,,,...n n s s s s s --10112,,,...,,,n n n s s s s s ---5. 平滑延拓在原信号两端用线性外插法补充采样值，即沿着信号两端包络线的一阶导数方向增加采样值。
小波变换C语言实现代码
#include <stdio.h> #include<math.h>#define height 256#define width 256void twoDDWT(double org[][width], double dwt[][width],int NumofBand); void oneDDWT(int flag, double org[][width], double dwt[][width]);void twoDIDWT(double dwt[][width], double rec[][width],int NumofBand); void oneDIDWT(int flag, double dwt[][width], double rec[][width]);int main()(int i, j;double org[height][width], dwt[height][width],rec[height][width];FILE *lena;FILE *lena3;FILE *lena4;lena = fopen("c:\\lena.raw", "rb");lena3 = fopen("c:\\lena3.raw", "wb");lena4 = fopen("c:\\lena4.raw", "wb");for (i = 0; i < height; i++) (for (j = 0; j < width; j++)( org[i][j] = fgetc(lena);))twoDDWT(org, dwt,4);for (i = 0; i < height/2; i++) (for (j = 0; j < width/2; j++)(org[i][j]=dwt[i][j];))twoDDWT(org, dwt,7);for (i = 0; i < height; i++) (for (j = 0; j < width; j++)( fputc(dwt[i][j], lena3);))twoDIDWT(dwt, rec,7);for (i = 0; i < height/2; i++)(for (j = 0; j < width/2; j++)( dwt[i][j]=rec[i][j];))twoDIDWT(dwt, rec,4);for (i = 0; i < height; i++)(for (j = 0; j < width; j++)( fputc(rec[i][j], lena4);))return 0;)void twoDDWT(double org[][width], double dwt[][width],int NumofBand)(if (NumofBand==4)(oneDDWT(0, org, dwt);oneDDWT(1, org, dwt);)if (NumofBand==7)(oneDDWT(2, org, dwt);oneDDWT(3, org, dwt);))void oneDDWT(int flag, double org[][width], double dwt[][width])(int i, j;double temp[height][width], temp1[height][width];double f_LPF[] = { -0.125, 0.25, 0.75, 0.25, -0.125 };double f_HPF[] = { -0.5, 1, -0.5 };if (flag == 0)//vertical direction{for (i = 0; i < height; i++){for (j = 0; j < width; j++){if (j == 0)temp[i][j] = f_LPF[0] * org[i][j + 2] + f_LPF[1] * org[i][j + 1] + f_LPF[2] * org[i][j] + f_LPF[3] * org[i][j + 1] + f_LPF[4] * org[i][j + 2];else if (j == 1)temp[i][j] = f_LPF[0] * org[i][j] + f_LPF[1] * org[i][j - 1] + f_LPF[2] * org[i][j] + f_LPF[3] * org[i][j + 1] + f_LPF[4] * org[i][j + 2];e se if (j > 1 && j < height-2)temp[i][j] = f_LPF[0] * org[i][j - 2] + f_LPF[1] * org[i][j - 1] + f_LPF[2] * org[i][j] + f_LPF[3] * org[i][j + 1] + f_LPF[4] * org[i][j + 2];else if (j == height-2)temp[i][j] = f_LPF[0] * org[i][j - 2] + f_LPF[1] * org[i][j - 1] + f_LPF[2] * org[i][j] + f_LPF[3] * org[i][j + 1] + f_LPF[4] * org[i][j];else if (j == height-1)temp[i][j] = f_LPF[0] * org[i][j - 2] + f_LPF[1] * org[i][j - 1] + f_LPF[2] * org[i][j] + f_LPF[3]* org[i][j - 1] + f_LPF[4] * org[i][j - 2];))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(dwt[i][j] = temp[i][j * 2 + 1];))for (i = 0; i < height; i++)(for (j = 0; j < width; j++)(if (j == 0)temp[i][j] = f_HPF[0] * org[i][j + 1] + f_HPF[1] * org[i][j] + f_HPF[2] * org[i][j + 1];else if (j > 0 && j < height-1)temp[i][j] = f_HPF[0] * org[i][j - 1] + f_HPF[1] * org[i][j] + f_HPF[2] * org[i][j + 1];else if (j == height-1)temp[i][j] = f_HPF[0] * org[i][j - 1] + f_HPF[1] * org[i][j] + f_HPF[2] * org[i][j - 1];))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(dwt[i][j + width/2] = temp[i][j * 2];)))if (flag == 1)//horizontal direction(for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(if (i == 0)temp[i][j] = f_LPF[0] * dwt[i + 2][j] + f_LPF[1] * dwt[i + 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i == 1)temp[i][j] = f_LPF[0] * dwt[i][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i > 1 && i < width-2)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i == width-2)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i][j];else if (i == width-1)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i - 1][j] + f_LPF[4] * dwt[i - 2][j];)for (i = 0; i < height; i++)(for (j = 0; j < width; j++)(temp1[i][j] = temp[i * 2 + 1][j];))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(if (i == 0)temp[i][j] = f_HPF[0] * dwt[i + 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i + 1][j];else if (i > 0 && i < width-1)temp[i][j] = f_HPF[0] * dwt[i - 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i + 1][j];else if (i == width-1)temp[i][j] = f_HPF[0] *dwt[i - 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i - 1][j];))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(temp1[i][j+width/2] = temp[i * 2][j];))for (i = 0; i < height; i++)(for (j = width/2; j < width; j++)(if (i == 0)temp[i][j] = f_LPF[0] * dwt[i + 2][j] + f_LPF[1] * dwt[i + 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i == 1)temp[i][j] = f_LPF[0] * dwt[i][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i > 1 && i < height-2)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i == height-2)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i][j];else if (i == height-1)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3]* dwt[i - 1][j] + f_LPF[4] * dwt[i - 2][j];))for (i = 0; i < height; i++)(for (j = width/2; j < width; j++)(temp1[i+height/2][j-width/2] = temp[i * 2 + 1][j];for (i = 0; i < height; i++)(for (j = width/2; j < width; j++)(if (i == 0)temp[i][j] = f_HPF[0] * dwt[i + 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i + 1][j];else if (i > 0 && i < height-1)temp[i][j] = f_HPF[0] * dwt[i - 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i + 1][j];else if (i == height-1)temp[i][j] = f_HPF[0] * dwt[i - 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i - 1][j];))for (i = 0; i < height; i++)(for (j = width/2; j < width; j++)( temp1[i +height/2][j] = temp[i * 2][j];))for (i = 0; i < height; i++)(for (j = 0; j < width; j++)(dwt[i][j] = temp1[i][j];)))if (flag == 2)//vertical direction(for (i = 0; i < height/2; i++)(for (j = 0; j < width; j++)(if (j == 0)temp[i][j] = f_LPF[0] * org[i][j + 2] + f_LPF[1] * org[i][j + 1] + f_LPF[2] * org[i][j] + f_LPF[3] * org[i][j + 1] + f_LPF[4] * org[i][j + 2];else if (j == 1)temp[i][j] = f_LPF[0] * org[i][j] + f_LPF[1] * org[i][j - 1] + f_LPF[2] * org[i][j] + f_LPF[3] * org[i][j + 1] + f_LPF[4] * org[i][j + 2];e se if (j > 1 && j < height/2-2)temp[i][j] = f_LPF[0] * org[i][j - 2] + f_LPF[1] * org[i][j - 1] + f_LPF[2] * org[i][j] + f_LPF[3]*org[i][j + 1] + f_LPF[4] * org[i][j + 2];else if (j == height/2-2)temp[i][j] = f_LPF[0] * org[i][j - 2] + f_LPF[1] * org[i][j - 1] + f_LPF[2] * org[i][j] + f_LPF[3]*org[i][j + 1] + f_LPF[4] * org[i][j];else if (j == height/2-1)temp[i][j] = f_LPF[0] * org[i][j - 2] + f_LPF[1] * org[i][j - 1] + f_LPF[2] * org[i][j] + f_LPF[3]*org[i][j - 1] + f_LPF[4] * org[i][j - 2];))for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(dwt[i][j] = temp[i][j * 2 + 1];))for (i = 0; i < height/2; i++)(for (j = 0; j < width/2; j++)(if (j == 0)temp[i][j] = f_HPF[0] * org[i][j + 1] + f_HPF[1] * org[i][j] + f_HPF[2] * org[i][j + 1];else if (j > 0 && j < height-1)temp[i][j] = f_HPF[0] * org[i][j - 1] + f_HPF[1] * org[i][j] + f_HPF[2] * org[i][j + 1];else if (j == height-1)temp[i][j] = f_HPF[0] * org[i][j - 1] + f_HPF[1] * org[i][j] + f_HPF[2] * org[i][j - 1];))for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(dwt[i][j + width/4] = temp[i][j * 2];)))if (flag == 3)//horizontal direction(for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(if (i == 0)temp[i][j] = f_LPF[0] * dwt[i + 2][j] + f_LPF[1] * dwt[i + 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i == 1)temp[i][j] = f_LPF[0] * dwt[i][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i > 1 && i < height/2-2)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i == height/2-2)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i][j];else if (i == height/2-1)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3]* dwt[i - 1][j] + f_LPF[4] * dwt[i - 2][j];))for (i = 0; i < height/4; i++)(for (j = 0; j < width/4; j++)(temp1[i][j] = temp[i * 2 + 1][j];))for (j = 0; j < width/4; j++)(if (i == 0)temp[i][j] = f_HPF[0] * dwt[i + 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i + 1][j];else if (i > 0 && i < height/2-1)temp[i][j] = f_HPF[0] * dwt[i - 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i + 1][j];else if (i == height/2-1)temp[i][j] = f_HPF[0] * dwt[i - 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i - 1][j];))for (i = 0; i < height/4; i++)(for (j = 0; j < width/4; j++)( temp1[i + height/4][j] = temp[i * 2][j];))for (i = 0; i < height/2; i++)(for (j = width/4; j < width/2; j++)(if (i == 0)temp[i][j] = f_LPF[0] * dwt[i + 2][j] + f_LPF[1] * dwt[i + 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i == 1)temp[i][j] = f_LPF[0] * dwt[i][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i > 1 && i < height/2-2)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i + 2][j];else if (i == height/2-2)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i + 1][j] + f_LPF[4] * dwt[i][j];else if (i == height/2-1)temp[i][j] = f_LPF[0] * dwt[i - 2][j] + f_LPF[1] * dwt[i - 1][j] + f_LPF[2] * dwt[i][j] + f_LPF[3] * dwt[i - 1][j] + f_LPF[4] * dwt[i - 2][j];))for (i = 0; i < height/4; i++)(for (j = width/4; j < width/2; j++)( temp1[i][j] = temp[i * 2 + 1][j];))for (i = 0; i < height/2; i++)(for (j = width/4; j < width/2; j++)(if (i == 0)temp[i][j] = f_HPF[0] * dwt[i + 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i + 1][j];else if (i > 0 && i < width/2-1)temp[i][j] = f_HPF[0] * dwt[i - 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i + 1][j];else if (i == width/2)temp[i][j] = f_HPF[0] * dwt[i - 1][j] + f_HPF[1] * dwt[i][j] + f_HPF[2] * dwt[i - 1][j];)for (i = 0; i < height/4; i++)(for (j = width/4; j < width/2; j++)(temp1[i + height/4][j] = temp[i * 2][j];))for (i = 0; i < height/2; i++)(for (j = 0; j < width/2; j++)(dwt[i][j] = temp1[i][j];))))void twoDIDWT(double dwt[][width], double rec[][width],int NumofBand) ( if (NumofBand==7)(oneDIDWT(2, dwt, rec);oneDIDWT(3, dwt, rec);)if (NumofBand==4)(oneDIDWT(1, dwt, rec);oneDIDWT(0, dwt, rec);))void oneDIDWT(int flag, double dwt[][width], double rec[][width])(int i, j;double temp[height][width], temp1[height][width];double i_LPF[] = { 0.5, 1,0.5 };double i_HPF[] = { -0.125, -0.25, 0.75, -0.25, -0.125 };if (flag == 2)//horizontal direction{for (i = 0; i < height/2; i++){for (j = 0; j < width/4; j++){if (i % 2 == 1)temp[i][j] = dwt[i / 2][j];elsetemp[i][j] = 0;}}for (i = 0; i < height/2; i++){for (j = 0; j < width/4; j++) {if (i == 0)rec[i][j] = i_LPF[0] * temp[i + 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i + 1][j];else if (i > 0 && i < height/2-1)rec[i][j] = i_LPF[0] * temp[i - 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i + 1][j];else if (i == height/2-1)rec[i][j] = i_LPF[0] * temp[i - 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i - 1][j];))for (i = height/4; i < height/2; i++)(for (j = 0; j < width/4; j++)(dwt[i - height/4][j] = dwt[i][j];))for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(if (i % 2 == 0)temp[i][j] = dwt[i / 2][j];elsetemp[i][j] = 0;))for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(if (i == 0)temp1[i][j] = i_HPF[0] * temp[i + 2][j] + i_HPF[1] * temp[i + 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i == 1)temp1[i][j] = i_HPF[0] * temp[i][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i > 1 && i < height/2-2)temp1[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i == height/2-2)temp1[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i][j];else if (i == height/2-1)temp1[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i - 1][j] + i_HPF[4] * temp[i - 2][j];))for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(rec[i][j] = rec[i][j] + temp1[i][j];))for (j = width/4; j < width/2; j++)(dwt[i][j - height/4] = dwt[i][j];))for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(if (i % 2 == 1)temp[i][j] = dwt[i / 2][j];elsetemp[i][j] = 0;))for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(if (i == 0)temp1[i][j] = i_LPF[0] * temp[i + 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i + 1][j];else if (i > 0 && i < height/2-1)temp1[i][j] = i_LPF[0] * temp[i - 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i + 1][j];else if (i == height/2-1)temp1[i][j] = i_LPF[0] * temp[i - 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i - 1][j];))for (i = height/4; i < height/2; i++)(for (j = width/4; j < width/2; j++)(dwt[i - height/4][j - width/4] = dwt[i][j];))for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(if (i % 2 == 0)temp[i][j] = dwt[i / 2][j];elsetemp[i][j] = 0;))for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(if (i == 0)dwt[i][j] = i_HPF[0] * temp[i + 2][j] + i_HPF[1] * temp[i + 1][j] + i_HPF[2] * temp[i][j] + HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i == 1)dwt[i][j] = i_HPF[0] * temp[i][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i > 1 && i < height/2-2)dwt[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i == height/2-2)dwt[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] *temp[i][j];else if (i == height/2-1)dwt[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i - 1][j] + i_HPF[4] * temp[i - 2][j];))for (i = 0; i < height/2; i++)(for (j = 0; j < width/4; j++)(rec[i][j + width/4] = temp1[i][j] + dwt[i][j];)))if (flag == 3)//vertical direction(for (i = 0; i < height/2; i++)(for (j = 0; j < width/2; j++)(if (j % 2 == 1)temp[i][j] = rec[i][j / 2];elsetemp[i][j] = 0;))for (i = 0; i < height/2; i++)(for (j = 0; j < width/2; j++)(if (j == 0)temp1[i][j] = i_LPF[0] * temp[i][j + 1] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i][j + 1];else if (j > 0 && j < height/2-1)temp1[i][j] = i_LPF[0] * temp[i][j - 1] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i][j + 1];else if (j == height/2-1)temp1[i][j] = i_LPF[0] * temp[i][j - 1] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i][j - 1];))for (i = 0; i < height/2; i++)(for (j = width/4; j < width/2; j++)(rec[i][j - width/4] = rec[i][j];))for (j = 0; j < width/2; j++)(if (j % 2 == 0)temp[i][j] = re可叩/2];elsetemp[i][j] = 0;))for (i = 0; i < height/2; i++)(for (j = 0; j < width/2; j++)(if (j == 0)rec[i][j] = i_HPF[0] * temp[i][j + 2] + i_HPF[1] * temp[i][j + 1] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i][j + 1] + i_HPF[4] * temp[i][j + 2];else if (j == 1)rec[i][j] = i_HPF[0] * temp[i][j] + i_HPF[1] * temp[i][j - 1] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i][j + 1] + i_HPF[4] * temp[而+ 2];else if (j > 1 && j < height/2-2)rec[i][j] = i_HPF[0] * temp[i][j - 2] + i_HPF[1] * temp[i][j - 1] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i][j + 1] + i_HPF[4] * temp[i][j + 2];else if (j == height/2-2)rec[i][j] = i_HPF[0] * temp[i][j - 2] + i_HPF[1] * temp[i][j - 1] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i][j + 1] + i_HPF[4] * temp[i][j];else if (j == height/2-1)rec[i][j] = i_HPF[0] * temp[i][j - 2] + i_HPF[1] * temp[i][j - 1] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i][j - 1] + i_HPF[4] * temp[i][j - 2];))for (i = 0; i <height/2; i++)(for (j = 0; j < width/2; j++)(rec[i][j] = rec[i][j] + temp1[i][j];)))if (flag == 1)//horizontal direction(for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(if (i % 2 == 1)temp[i][j] = dwt[i / 2][j];elsetemp[i][j] = 0;))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(if (i == 0)rec[i][j] = i_LPF[0] * temp[i + 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i + 1][j];else if (i > 0 && i < height-1)rec[i][j] = i_LPF[0] * temp[i - 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i + 1][j];else if (i == height-1)rec[i][j] = i_LPF[0] * temp[i - 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i - 1][j];))for (i = 0; i < height/2; i++)(for (j = width/2; j < width; j++)(dwt[i][j - height/2] = dwt[i][j];))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(if (i % 2 == 0)temp[i][j] = dwt[i / 2][j];elsetemp[i][j] = 0;))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(if (i == 0)temp1[i][j] = i_HPF[0] * temp[i + 2][j] + i_HPF[1] * temp[i + 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i == 1)temp1[i][j] = i_HPF[0] * temp[i][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i > 1 && i < height-2)temp1[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i == height-2)temp1[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i][j];else if (i == height-1)temp1[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i - 1][j] + i_HPF[4] * temp[i - 2][j];))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(rec[i][j] = rec[i][j] + temp1[i][j];for (i = height/2; i <height; i++)(for (j = 0; j < width/2; j++)(dwt[i - width/2][j] = dwt[i][j];))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(if (i % 2 == 1)temp[i][j] = dwt[i / 2][j];elsetemp[i][j] = 0;))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(if (i == 0)temp1[i][j] = i_LPF[0] * temp[i + 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i + 1][j]; else if (i > 0 && i < height-1)temp1[i][j] = i_LPF[0] * temp[i - 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i + 1][j];else if (i == height-1)temp1[i][j] = i_LPF[0] * temp[i - 1][j] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i - 1][j]; ))for (i = height/2; i < height; i++)(for (j = width/2; j < width; j++)(dwt[i - height/2][j - width/2] = dwt[i][j];))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(if (i % 2 == 0)temp[i][j] = dwt[i / 2][j];elsetemp[i][j] = 0;))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(if (i == 0)dwt[i][j] = i_HPF[0] * temp[i + 2][j] + i_HPF[1] * temp[i + 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i == 1)dwt[i][j] = i_HPF[0] * temp[i][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i > 1 && i < height-2)dwt[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] * temp[i + 2][j];else if (i == height-2)dwt[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i + 1][j] + i_HPF[4] *temp[i][j];else if (i == height-1)dwt[i][j] = i_HPF[0] * temp[i - 2][j] + i_HPF[1] * temp[i - 1][j] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i - 1][j] + i_HPF[4] * temp[i - 2][j];))for (i = 0; i < height; i++)(for (j = 0; j < width/2; j++)(rec[i][j + width/2] = temp1[i][j] + dwt[i][j];)))if (flag == 0)//vertical direction(for (i = 0; i < height; i++)(for (j = 0; j < width; j++)(if (j % 2 == 1)temp[i][j] = rec[i][j / 2];elsetemp[i][j] = 0;))for (i = 0; i < height; i++)(for (j = 0; j < width; j++)(if (j == 0)temp1[i][j] = i_LPF[0] * temp[i][j + 1] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i][j + 1];else if (j > 0 && j < height-1)temp1[i][j] = i_LPF[0] * temp[i][j - 1] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i][j + 1];else if (j == height-1)temp1[i][j] = i_LPF[0] * temp[i][j - 1] + i_LPF[1] * temp[i][j] + i_LPF[2] * temp[i][j - 1];))for (i = 0; i < height; i++)(for (j = width/2; j < width; j++)(rec[i][j - width/2] = rec[i][j];for (i = 0; i < height; i++)(for (j = 0; j < width; j++)(if (j % 2 == 0)temp[i][j] = rec[i][j / 2];elsetemp[i][j] = 0;))for (i = 0; i < height; i++)(for (j = 0; j < width; j++)(if (j == 0)rec[i][j] = i_HPF[0] * temp[i][j + 2] + i_HPF[1] * temp[i][j + 1] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i][j + 1] + i_HPF[4] * temp[i][j + 2];else if (j == 1)rec[i][j] = i_HPF[0] * temp[i][j] + i_HPF[1] * temp[i][j - 1] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i][j + 1] + i_HPF[4] * temp[而+ 2];else if (j > 1 && j < height-2)rec[i][j] = i_HPF[0] * temp[i][j - 2] + i_HPF[1] * temp[i][j - 1] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i][j + 1] + i_HPF[4] * temp[i][j + 2];else if (j == height-2)rec[i][j] = i_HPF[0] * temp[i][j - 2] + i_HPF[1] * temp[i][j - 1] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i][j + 1] + i_HPF[4] * temp[i][j];else if (j == height-1)rec[i][j] = i_HPF[0] * temp[i][j - 2] + i_HPF[1] * temp[i][j - 1] + i_HPF[2] * temp[i][j] + i_HPF[3] * temp[i][j - 1] + i_HPF[4] * temp[i][j - 2];))for (i = 0; i <height; i++)(for (j = 0; j < width; j++)(rec[i][j] = rec[i][j] + temp1[i][j];。
小波变换课件ch4 Mallat算法及二维小波47页PPT
Hale Waihona Puke 46、我们若已接受最坏的，就再没有什么损失。——卡耐基 47、书到用时方恨少、事非经过不知难。——陆游 48、书籍把我们引入最美好的社会，使我们认识各个时代的伟大智者。——史美尔斯 49、熟读唐诗三百首，不会作诗也会吟。——孙洙 50、谁和我一样用功，谁就会和我一样成功。——莫扎特
小波变换课件ch4 Mallat算法及二维 小波
1、纪律是管理关系的形式。——阿法 纳西耶 夫 2、改革如果不讲纪律，就难以成功。
3、道德行为训练，不是通过语言影响 ，而是 让儿童 练习良 好道德 行为， 克服懒 惰、轻 率、不 守纪律 、颓废 等不良 行为。 4、学校没有纪律便如磨房里没有水。 ——夸 美纽斯
小波变换快速算法及应用小结
离散小波变换的快速算法Mallat算法[经典算法]在小波理论中,多分辨率分析是一个重要的组成部分。
多分辨率分析是一种对信号的空间分解方法,分解的最终目的是力求构造一个在频率上高度逼近L2(R)空间的正交小波基,这些频率分辨率不同的正交小波基相当于带宽各异的带通滤波器。
因此,对于一个能量有限信号,可以通过多分辨率分析的方法把其中的逼近信号和细节信号分离开,然后再根据需要逐一研究。
多分辨率分析的概念是在构造正交小波基的时候提出的,并同时给出了著名的Mallat算法。
Mallat算法在小波分析中的地位相当于快速傅立叶变换在经典傅立叶变换中的地位,为小波分析的应用和发展起到了极大的推动作用。
MALLAT算法的原理在对信号进行分解时，该算法采用二分树结构对原始输入信号x(n)进行滤波和二抽取，得到第一级的离散平滑逼近和离散细节逼近，再采用同样的结构对进行滤波和二抽取得到第二级的离散平滑逼近和离散细节逼近，再依次进行下去从而得到各级的离散细节逼近对,…，即各级的小波系数。
重构信号时，只要将分解算法中的步骤反过来进行即可，但要注意，此时的滤波器与分解算法中的滤波器不一定是同一滤波器，并且要将二抽取装置换成二插入装置才行。
多孔算法[小波变换快速算法及其硬件实现的研究毛建华]多孔算法是由于1992年提出的一种利用Mallat算法结构计算小波变换的快速算法，因在低通滤波器和高通滤波器中插入适当数目的零点而得名。
它适用于的二分树结构，与Mallat算法的电路实现结构相似。
先将Mallat算法的电路实现的基本支路作一下变形。
令的z变换为与，下两条支路完全等价，只不过是将插值和二抽取的顺序调换一下罢了。
图中其它的上下两条支路也为等效支路，可仿照上面的方法证明。
这样，我们便可由Mallat算法的二分树电路结构得出多孔算法的电路级联图，原Mallat算法中的电路支路由相应的等效支路所取代，所以整个电路形式与Mallat算法非常相似。
小波mallat算法
⼩波mallat算法算法要求：输⼊序列是⼤于滤波器长度的偶数列确实可以通过编程的⼿段使算法适合所有的情况，但本⽂章的⽬的是展⽰mallat算法的过程，所以就⼀切从简了// Mallat.cpp : Defines the entry point for the console application.//#include "stdafx.h"#include "stdio.h"/*mallat算法分解* dSIn 输⼊的序列s，dH0尺度函数展开系数，dH1⼩波函数展开系数，dSOut输出低频部分，dDOut输出⾼频部分，* nSIn_Len 输⼊序列的长度，nH_Len 滤波器的长度。
*/int DwtFun(double *pdSIn,double *pdH0,double *pdH1,double *pdSOut,double *pdDOut,int nSIn_Len,int nH_Len) {int i,j,k;//延拓后的Len是⼀个本体长度加⼀个滤波器长度int nLen=nSIn_Len+2*nH_Len;//建⽴滤波前的序列pdSArray，滤波后的序列pdSAOut低频部分，pdDAOut⾼频部分double *pdSArray=new double[nLen];double *pdSAOut=new double[nLen];double *pdDAOut=new double[nLen];//对称延拓for(i=0;i<nLen;i++){if(i<nH_Len){pdSArray[i]=pdSIn[nH_Len-i-1];}else if(i>=nH_Len+nSIn_Len){pdSArray[i]=pdSIn[nH_Len+2*nSIn_Len-1-i];}else{pdSArray[i]=pdSIn[i-nH_Len];}}//求输出序列低频部分dSOut,⾼频部分dDOut.i->nLen,k->nH_Lendouble dSTemp,dDTemp;for(i=0;i<nLen;i++){dSTemp=0.0;dDTemp=0.0;for(k=0;k<nH_Len;k++){if((i-k)<0)continue;else{//低频部分dSTemp+=pdH0[nH_Len-k-1]*pdSArray[i-k];//⾼频部分dDTemp+=pdH1[nH_Len-k-1]*pdSArray[i-k];}}pdSAOut[i]=dSTemp;pdDAOut[i]=dDTemp;}//⼆抽取.先将pdSAOut前nH_Len长的⼀段舍弃,抽取偶数列for(i=nH_Len,j=0;i<nLen;i+=2,j++){pdSOut[j]=pdSAOut[i+1];pdDOut[j]=pdDAOut[i+1];}//返回输出序列的长度return j;delete pdSArray;pdSArray=NULL;delete pdSAOut;pdSAOut=NULL;delete pdDAOut;pdDAOut=NULL;}/*mallat 算法重构* psSIn 输⼊的低频序列，pdDIn输⼊的⾼频序列，g0，g1重构滤波器，pdOut输出序列，nSInLen输⼊序列的长度* nG_Len 滤波器长度*/int IDwtFun(double *pdSIn,double *pdDIn,double *pdG0,double *pdG1,double *pdOut,int nSInLen,int nG_Len) {int i,j,k;//建⽴⼀个数列存放插⼊后的数列int nTemp=2*nSInLen;double *pdInSertS=new double[nTemp];double *pdInSertD=new double[nTemp];//⼆插⼊j=0;for(i=0;i<nTemp;i++){if(i%2==0){pdInSertS[i]=0;pdInSertD[i]=0;}else{pdInSertS[i]=pdSIn[j];pdInSertD[i]=pdDIn[j];j++;}}//对称拓延//创建⼀个nTemp+nG_Len长的数列int nLen=nTemp+2*nG_Len;double *pdSAIn=new double[nLen];double *pdDAIn=new double[nLen];for(i=0;i<nLen;i++){if(i<nG_Len){pdSAIn[i]=pdInSertS[nG_Len-i-1];pdDAIn[i]=pdInSertD[nG_Len-i-1];}else if(i==nTemp+nG_Len){pdSAIn[i]=0.0;pdDAIn[i]=0.0;}else if(i>nTemp+nG_Len){pdSAIn[i]=pdInSertS[nG_Len+2*nTemp-i-1];pdDAIn[i]=pdInSertD[nG_Len+2*nTemp-i-1];}else{pdSAIn[i]=pdInSertS[i-nG_Len];pdDAIn[i]=pdInSertD[i-nG_Len];}}//⽤滤波器G0和G1对数列进⾏滤波double *pdSAOut=new double[nLen];double *pdDAOut=new double[nLen];double dSTemp,dDTemp;for(i=0;i<nLen;i++){dSTemp=0.0;dDTemp=0.0;for(k=0;k<nG_Len;k++){if((i-k)<0)continue;else{//低频部分dSTemp+=pdG0[nG_Len-k-1]*pdSAIn[i-k];//⾼频部分dDTemp+=pdG1[nG_Len-k-1]*pdDAIn[i-k];}}pdSAOut[i]=dSTemp;pdDAOut[i]=dDTemp;}//合并低频，⾼频for(i=2*nG_Len-1,j=0;i<nLen;i++,j++){pdOut[j]=pdSAOut[i]+pdDAOut[i];}return j;delete pdInSertS;pdInSertS=NULL;delete pdInSertD;pdInSertD=NULL;delete pdSAIn;pdSAIn=NULL;delete pdDAIn;pdDAIn=NULL;delete pdSAOut;pdSAOut=NULL;delete pdDAOut;pdDAOut=NULL;}int main(int argc, char* argv[]){int i;//db4⼩波，已经取反 h0,h1是分解滤波器，g0,g1是重构滤波器double dDb4h0[] = { 0.2303778133088964, 0.7148465705529154,0.6308807679398587, -0.0279837694168599,-0.1870348117190931, 0.0308413818355607,0.0328830116668852, -0.0105974017850690 };double dDb4h1[] = { -0.0105974017850690 , -0.0328830116668852, 0.0308413818355607 , 0.1870348117190931,-0.0279837694168599 , -0.6308807679398587,0.7148465705529154 , -0.2303778133088964};double dDb4g0[] = { -0.0105974017850690 , 0.0328830116668852,0.0308413818355607 , -0.1870348117190931,-0.0279837694168599 , 0.6308807679398587,0.7148465705529154 , 0.2303778133088964};double dDb4g1[] = { -0.2303778133088964 , 0.7148465705529154,-0.6308807679398587 , -0.0279837694168599,0.1870348117190931 , 0.0308413818355607,-0.0328830116668852 , -0.0105974017850690};//⽣成⼀个数列，本算法要求输⼊的数列是⽐滤波器长的偶数列double a[]={1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0}; //double a[]={1.0,4.0,5.5,8.2,2.7,5.2,2.0,2.0,2.0,3.0,3.0,4.0,4.0,14.0,17.0,11.0};//输出double *pdS=new double[100];double *pdD=new double[100];double *pdOut=new double[100];int l=DwtFun(a,dDb4h0,dDb4h1,pdS,pdD,16,8);for(i=0;i<l-1;i++){printf("%f\t",pdS[i]);printf("\n");}printf("*********************\n");for(i=0;i<l-1;i++){printf("%f\t",pdD[i]);printf("\n");}printf("*********************\n");int v=IDwtFun(pdS,pdD,dDb4g0,dDb4g1,pdOut,11,8);//i<v-nG_Len+1for(i=0;i<v-7;i++){printf("%f\t",pdOut[i]);printf("\n");}delete []pdS;pdS=NULL;delete []pdD;pdD=NULL;delete []pdOut;pdOut=NULL;return 0;}。
小波变换的实现技术
用Haar小波进行分解，
f 8 t P8 f t V 8
f8 t f7 t d 7 t f6 t d 6 t d 7 t f5 t d 5 t d 6 t d 7 t
c
J
用小波处理离散信号的基本步骤
b n n 0,1, , N
f t
1
其采样间距为
N
1
N 2 ，
L
n
a n L ,n t V L
L
使得 b n f N 1 n
N
1/ 2
a n bn
L
a
L n
对 a n 做小波分解、对小波系数处理以及对处理后的系数进行小波重构等
特点:
1) 能够实现重构. 2) 难以用于数据 压缩应用
• idwt()
X = idwt(cA,cD,Lo_R,Hi_R) X = idwt(cA,cD,Lo_R, Hi_R ，'mode',MODE)
对于周期延拓方法， lx 2 la 对于其他延拓方式， lx 2 la lf 2
1
)
小波分解与重构的多相位表示
• 滤波器的多相位矩阵
滤波器 h 和 g 的多相位矩阵为：
he ( z ) P(z) ho ( z )
he(z) P(z) ~ ho (z)
~ ~
ge(z) go ( z)
~ g o ( z)
~
滤波器 h 和 g 的对偶多相位矩阵为：
1) 为什么称为多孔算法（a’trous algorithm） ?
第7章小波变换2
n
(7.4)
j ,k hn 2 k j 1,n
c j ,k f , j ,k
j ,k g n 2 k c j 1,n
n
n
f , hn j 1,2 k n hn f , j 1,2 k n
证明如下：
g n j 1,2 k n n
n
x 2 g n 2 x n (7.3) n
n
x 2 hn 2 x n
（7.4）
j ,k 2 j / 2(2 j x k )
n
k k
VJ span{ Jk } Ck Jk
k
Mallat算法
一维多分辨分析Mallat算法
目标：
求出 （7.2）式中的 d j ,k ; c jk 寻找一种快速算法而不是内积
思路： 用上一级的c j1,k ,d j 1,k 求下一级的c j,k , d j,k
f x
J M j J 1
d
kZ
j ,k
j ,k cJ M ,k J M ,k
kZ
（ 7.2)
将双尺度方程 x an 2 x n 和小波方程
x 1 an1 2x n 写成如下形式：
n n
n
x 2 g n 2 x n n
n
x 2 hn 2 x n
(7.3)
an 与 hn、gn 的关系为 其中，
hn an 2
g n 1
n
a n 1 2
没必要记
第4章+小波变换的实现技术
或者查阅一些中文资料。
第17页 共22页
用小波处理函数/信号的基本步骤
用小波处理信号 f (t) 的基本过程包括：
1、初始化 2、小波变换 3、小波系数处理 4、小波逆变换
第9页 共22页
Mallat算法
l / 2
第10页 共22页
Mallat算法
设原始信号长为m，滤波器长为l，则延拓总长为l-1即可。 可以采用左、右、两端三种不同的方式。
（1）两端同时延拓：
第一种情况：当l为奇数时，每端延拓长度为 l / 2（下取整数）； 第二种情况：当l为偶数时，一端延拓长度为l/2，另一端延拓长度为l/2-1
在原信号两端补零
周期延拓
, sn1, s0 , s1, s2 ,, sn1, s0 , s1, s2 ,, sn1, s0 , s1, s2 ,, sn1, s0 ,
将信号看成一个周期信号
第7页 共22页
Mallat算法实现中的问题
1、边界延拓方法
周期对称延拓法
（1）与信号作卷积的滤波器长度为奇数
f t anLL,n t VL n
使得 bn f N 1n
L
N 1/ 2anL bn 2 2 anL bn
anL
对 anL 做小波分解、对小波系数处理以及对处理后的系数进行小波重构等 说明: 1) 对 bn 做小波分解，如何?
2) 若 bn 的采样间距为1，如何?
第23页 共22页
, sn2 , sn1,, s2 , s1, s0 , s1, s2 ,, sn1, sn2 ,, s1, s0 , s1,
小波变换VC6_0程序实现
(1) (2)
3 一维小波变换 VC++实现
由 2 可知 ， 离散小波变换实际是通过卷积运算完成的 ， 小波变换一维 DWT 或 图1 离散小波变换计算
IDWT 实现代码如下 ： 1) 函数输入 double * data ， 指向源数据的指针 。 int nCurLongth ， 当前处理数据长度 。 int IDWT ， 是否为 DWT ，1 表示为 IDWT ，0 表示 DWT 。 int nStep ， 当前分解层数 。 int nSupp ， 小波基的紧支集的长度 。 本文中 ， 小波基存储在 hCoef 这个二维数组中 ， 通过 nSupp 调用小波基 double s = sqrt(2); double* h = NULL;
Wavelet Transform Procedures for Implementation VC++6.0 LIU Chao1, XING Shu-guang1, YANG Xi-e2 (1.China University of Geosciences (Wuhan) The Faculty of Earth Resources, Wuhan 430074, China; 2.Northwest A&F University, Col lege of Humanities, Yangling 712100, China) Abstract: With the wavelet-depth study, the scope of application of wavelet transform more and more widely, but most of the research, the work is the use of Matlab wavelet toolkit for programming. This paper uses the Design and Implementation of VC++6.0 digital image processing software modules based on wavelet transform to prepare concrete realization of programs, design software interface to make it a convenient, for helping future use of wavelet transform to the practical work. Key words: wavelet transform; VC++6.0; Matlab
离散小波变换算法剖析及其通用程序实现
Abstract：For the application software development of wavelet transform (WT), some key steps and techniques for programming are introduced based on Mallat algorithm analysis and descriptions. Many data are employed to check the proposed algorithm programmed by the language of Delphi. The interval and final results shows no differences from the outputs produced by WT programs with the Matlab wavelet toolbox. The main algorithms are written as sub-functions that are helpful for the researchers to convert their achievements in Matlab to other computer languages and construct the WT software for specialized application independently. Key words：discrete wavelet transform; Mallat algorithm; algorithm description; programming; software development
第5章小波变换的matlab实现精品资料课件
· A3=wrcoef('a',C,L,'db1',3); · D1=wrcoef('d',C,L,'db1',1); · D2=wrcoef('d',C,L,'db1',2); · D3=wrcoef('d',C,L,'db1',3);
Approximation
A3
Detail D1
Detail D2
Wavelet Packet 2-D
Multiple 1-D Mutisignal Analysis 1 - D
Mutivariate Denoising Mutiscale Princ. Comp.Analysis
Wavelet Design New Wavelet for CWT
Specialized Tools 1-D SWT Denoising 1-D Density Estimation 1-D
低频系数
原始信号 高频系数
系数重构
·命令： upcoef
·格式：
1. Y=upcoef(O,X,'wname',N) 2.Y=upcoef(O,X,'wname',N,L) 3.Y=upcoef(O,X,'Lo_R,Hi_R',N)
4. Y=upcoef(O,X,'Lo R,Hi R ',N,L)
图形接口方式 (GUI)
let Toolbox Eile Window Help
Tain
Tenu
One-Dimensional Wavelet 1-D
Wavelet Packet 1-D Continuous Wavelet 1-D Complex Continuous WaveleWavelet 2 -D
二维离散小波分解的C语言实现 论文
高等教育自学考试毕业论文（设计）题目：二维离散小波分解的C语言实现摘要小波变换用于图像处理是小波变换应用效果比较突出的领域之一。
由于图像是二维信号，因此首先需要把小波变换由一维推广到二维。
本文在一维离散Mallat算法的基础上，用C语言实现了二维图像的离散小波变换。
这种二维变换是行列可分离的变换方式，即二维分解可以通过行和列依次作一维分解实现。
对图像作二维离散小波分解后得到一个低频子带和一系列高频子带，分别反映图像的基本信息和细节信息。
用这些子带也可以实现图像的重构。
目录第一章绪论 (1)1. 1小波理论与应用技术的发展概况 (1)1. 2图像技术的发展历程及面临的问题 (2)1. 3小波的特点及其在图像处理中的应用 (2)第二章Mallat算法由一维到二维的推广 (4)2. 1小波级数 (4)2. 2 Mallat算法 (5)2. 3二维离散小波变换 (7)2. 4二维离散小波变换后的系数分布 (8)第三章二维Mallat算法的C语言实现 (10)3. 1基本模块 (10)3．2 单层分解与重构 (10)3．3金字塔结构的多层分解和重构 (11)3．4小波系数的数据结构 (14)3．5 结果与分析 (14)参考文献 (19)致谢 (20)第一章绪论1. 1小波理论与应用技术的发展概况小波分析是当前数学中一个迅速发展的新领域，它同时具有理论深刻和应用十分广泛的双重意义。
小波分析的应用是与小波分析的理论研究紧密地结合在一起的。
现在，它已经在科技信息产业领域取得了令人瞩目的成就。
电子信息技术是六大高新技术中重要的一个领域，它的重要方面是图像和信号处理。
现今，信号处理已经成为当代科学技术工作的重要部分，信号处理的目的就是：准确的分析、诊断、编码压缩和量化、快速传递或存储、精确地重构(或恢复)。
从数学地角度来看，信号与图象处理可以统一看作是信号处理(图像可以看作是二维信号)，在小波分析地许多分析的许多应用中，都可以归结为信号处理问题。
小波分析与实例
小波分析的基本知识—二进小波变换
2、小波分析的基本知识—二进小波变换
定义：设yj，k(t)∈L2(R)，且满足 (1.64) 由此得到的小波yj，k(t)称为二进正交小波。
3、多尺度分析与Mallat算法
多分辨分析
为了改变信号的分辨率使得人们可以根据特定的目标处理相关的细节，1983年，与在计算机视觉的应用中引进了一个能够处理低分辨率图像，同时根据需要进一步提高图像分辨率的多分辨率Laplace塔式算法。1986年Mallat和Meyer构造了多分辨分析公式。随着多分辨分析的出现，构造小波的困难得到了较圆满的解决。为了对信号进行较高分辨率的处理，需要一种所谓的“增量信息”。为此，Mallat选用正交小波基作为对“增量信息”进行数学描述，并最终发展成为了多分辨分析。
load noissin c = cwt(noissin,1:48,'db4'); c = cwt(noissin,1:48,'db4','plot'); c = cwt(noissin,2:2:128,'db4','plot');
3、多尺度分析与Mallat算法
3、多尺度分析与Mallat算法
2、小波分析的基本知识—连续小波变换
2、小波分析的基本知识—连续小波变换
小波变换的系数如图所示的灰度值图表征，横坐标表示变换系数的系号，纵坐标表示尺度，灰度颜色越深，表示系数的值越大。
小波分析的基本知识—离散小波变换
离散小波变换：
在实际运用中，尤其是在计算机上实现，连续小波必须加以离散化。因此，有必要讨论一下连续小波ya，b(t)和连续小波变换Wf(a，b)的离散化。需要强调指出的是，这一离散化都是针对连续的尺度参数a和连续平移参数b的，而不是针对时间变量t的。 在连续小波中，考虑函数 这里，b∈R，a∈R＋，且a≠0，y是容许的，为方便起见，在离散化中，总限制a只取正值，这样相容性条件就变为