小波算法c语言

1. 小波算法是什么

王卫国 郭宝龙

（西安电子科技大学机电工程学院，西安 710071）

摘 要 随着互联网的普及和图象应用范围的不断扩大,对图象的编码提出了新的要求,即不仅要求具有高的压缩比,还要求有许多新的功能,如渐进编解码、从有损压缩到无损压缩等。嵌入式零树小波编码较好地实现了这一思想,因此奠定了它在图象编码中的地位。近年来，在嵌入式零树小波编码(EZW)算法的基础上出现了许多新的改进算法，如多级树集合分裂算法(SPIHT),集合分裂嵌入块编码(SPECK),可逆的嵌入小波压缩法(CREW)等.本文对这些算法从原理到性能进行了比较和讨论,说明了嵌入式图象编码的研究方向。

关 键 词 图象编码 嵌入式 零树 小波变换

On Embedded Zerotree Wavelets Coding and other Improved Algorithms
WANG Wei-guo, GUO Bao-long

(School of Mechano-Electronic Engineering,Xidian Univ.,Xi’an 710071)

Abstract With the extensive application of internet and image,some new requirements on image coding,such as high compression rate ,pregressive codec,and compression from lossy to lossless ,are to be satisfied.These functions can be performed well by EZW(Embedded Zerotree Wavelets) coding.On the bases of EZW,many newly improved algorithms have been developed in recent years.They can illustrated by algorithms like SPIHT(Set Partitioning in Hierarchical Trees),SPECK(Set Partitioned Embedded block coder),In this paper,the writer discusses the principles and performances of these algorithms,thus explains the research tendency in the area of embedded image coding.

Keywords Image coding,Embedded,Zerotree,Wavelet transform

0. 引言

在基于小波变换的图象压缩方案中,嵌入式零树小波 EZW(Embedded Zerotree Wavelets)[1]编码很好地利用小波系数的特性使得输出的码流具有嵌入特性。它的重要性排序和分级量化的思想被许多编码算法所采用。近年来,在对EZW改进的基础上，提出了许多新的性能更好的算法,如多级树集合分裂算法(SPIHT :Set Partitioning In Hierarchical Trees)[2],集合分裂嵌入块编码(SPECK：Set Partitioned Embedded bloCK coder),可逆嵌入小波压缩算法(CREW:Compression with Reversible Embedded Wavelets)[3] 。本文对这些算法进行了原理分析、性能比较,说明了嵌入式小波图象编码的研究方向。

2. c语言怎么实现小波变换

#include<stdio.h>
#include<stdlib.h>
#defineLENGTH512//信号长度
/******************************************************************
*一维卷积函数
*
*说明:循环卷积,卷积结果的长度与输入信号的长度相同
*
*输入参数:data[],输入信号;core[],卷积核;cov[],卷积结果;
*			n,输入信号长度;m,卷积核长度.
*
*李承宇,[email protected]
*
*2010-08-18
******************************************************************/
voidCovlution(doubledata[],doublecore[],doublecov[],intn,intm)
{
	inti=0;
	intj=0;
	intk=0;

	//将cov[]清零
	for(i=0;i<n;i++)
	{
		cov[i]=0;
	}

	//前m/2+1行
	i=0;
	for(j=0;j<m/2;j++,i++)
	{
		for(k=m/2-j;k<m;k++)
		{
			cov[i]+=data[k-(m/2-j)]*core[k];//k针对core[k]
		}

		for(k=n-m/2+j;k<n;k++)
		{
			cov[i]+=data[k]*core[k-(n-m/2+j)];//k针对data[k]
		}
	}

	//中间的n-m行
	for(i=m/2;i<=(n-m)+m/2;i++)
	{
		for(j=0;j<m;j++)
		{
			cov[i]+=data[i-m/2+j]*core[j];
		}
	}

	//最后m/2-1行
	i=(n-m)+m/2+1;
	for(j=1;j<m/2;j++,i++)
	{
		for(k=0;k<j;k++)
		{
			cov[i]+=data[k]*core[m-j-k];//k针对data[k]
		}

		for(k=0;k<m-j;k++)
		{
			cov[i]+=core[k]*data[n-(m-j)+k];//k针对core[k]
		}
	}

}

/******************************************************************
*一维小波变换函数
*
*说明:一维小波变换,只变换一次
*
*输入参数:input[],输入信号;output[],小波变换结果，包括尺度系数和
*小波系数两部分;temp[],存放中间结果;h[],Daubechies小波基低通滤波器系数;
*g[],Daubechies小波基高通滤波器系数;n,输入信号长度;m,Daubechies小波基紧支集长度.
*
*李承宇,[email protected]
*
*2010-08-19
******************************************************************/
voidDWT1D(doubleinput[],doubleoutput[],doubletemp[],doubleh[],
		doubleg[],intn,intm)
{
//	doubletemp[LENGTH]={0};//?????????????
	
	inti=0;
/*
	//尺度系数和小波系数放在一起
	Covlution(input,h,temp,n,m);

	for(i=0;i<n;i+=2)
	{
		output[i]=temp[i];
	}

	Covlution(input,g,temp,n,m);

	for(i=1;i<n;i+=2)
	{
		output[i]=temp[i];
	}
*/

	//尺度系数和小波系数分开
	Covlution(input,h,temp,n,m);

	for(i=0;i<n;i+=2)
	{
		output[i/2]=temp[i];//尺度系数
	}

	Covlution(input,g,temp,n,m);

	for(i=1;i<n;i+=2)
	{
		output[n/2+i/2]=temp[i];//小波系数
	}

}

voidmain()
{

	doubledata[LENGTH];//输入信号
	doubletemp[LENGTH];//中间结果
	doubledata_output[LENGTH];//一维小波变换后的结果
	intn=0;//输入信号长度
	intm=6;//Daubechies正交小波基长度
	inti=0;
	chars[32];//从txt文件中读取一行数据

	staticdoubleh[]={.332670552950,.806891509311,.459877502118,-.135011020010,
					-.085441273882,.035226291882};
	staticdoubleg[]={.035226291882,.085441273882,-.135011020010,-.459877502118,
					.806891509311,-.332670552950};
	//读取输入信号
	FILE*fp;
	fp=fopen("data.txt","r");
	if(fp==NULL)//如果读取失败
	{
		printf("错误！找不到要读取的文件/"data.txt/"/n");
		exit(1);//中止程序
	}

	while(fgets(s,32,fp)!=NULL)//读取长度n要设置得长一点，要保证读到回车符，这样指针才会定位到下一行？回车符返回的是零值？是，非数字字符经过atoi变换都应该返回零值
	{
	//	fscanf(fp,"%d",&data[count]);//一定要有"&"啊！！！最后读了个回车符！适应能力不如atoi啊
		data[n]=atof(s);
		n++;
	}

	//一维小波变换
	DWT1D(data,data_output,temp,h,g,n,m);

	//一维小波变换后的结果写入txt文件
	fp=fopen("data_output.txt","w");

	//打印一维小波变换后的结果
	for(i=0;i<n;i++)
	{
		printf("%f/n",data_output[i]);
		fprintf(fp,"%f/n",data_output[i]);
	}

	//关闭文件
	fclose(fp);
}

3. 小波算法

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

4. 我想用多贝西小波对有奇异点的信号进行分析，然后提取奇异点发生的时间点信息，用C语言怎么写

(1)小波模极大值重构 MATLAB代码_天天向上_新浪博客 http://blog.sina.com.cn/s/blog_6c00b0e30100u1s5.html

function
[signal,swa,swd,ddw,wpeak]=wave_peak(points,level,Lo_D,Hi_D,Lo_R,Hi_R,offset)
%
该函数用于读取ecg信号，找到小波变换模极大序列

warning off;
ecgdata=load('ecg.txt');
%需要分析的信号，自己加
plot(ecgdata(1:points)),grid on,axis
tight,axis([1,points,-50,300]);
signal=ecgdata(1:points)'+offset;

% 信号的小波变换，按级给出概貌和细节的波形
[swa,swd] =
swt(signal,level,Lo_D,Hi_D);
figure;
subplot(level,1,1);
plot(real(signal)); grid on;axis tight;
for i=1:level
subplot(level+1,2,2*(i)+1);
plot(swa(i,:)); axis
tight;grid on;xlabel('time');
ylabel(strcat('a
',num2str(i)));
subplot(level+1,2,2*(i)+2);
plot(swd(i,:)); axis
tight;grid on;
ylabel(strcat('d ',num2str(i)));
end

%求小波变换的模极大值及其位置
ddw=zeros(size(swd));
pddw=ddw;
nddw=ddw;
posw=swd.*(swd>0);
pdw=((posw(:,1:points-1)-posw(:,2:points))<0);
pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))>0);
negw=swd.*(swd<0);
ndw=((negw(:,1:points-1)-negw(:,2:points))>0);
nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0);
ddw=pddw|nddw;
ddw(:,1)=1;
ddw(:,points)=1;
wpeak=ddw.*swd;
wpeak(:,1)=wpeak(:,1)+1e-10;
wpeak(:,points)=wpeak(:,points)+1e-10;

%按级给出小波变换模极大的波形
figure;
for i=1:level
subplot(level,1,i);
plot(wpeak(i,:)); axis tight;grid
on;
ylabel(strcat('j= ',num2str(i)));
end

注：运行此程序时一定要将待处理信号添加进去，程序中的红色部分。
追问：
ecgdata=load('ecg.txt');

5. 离散小波变换用c语言怎么写呀没有思路啊。求高人指点，一定感激不尽。

得有些理论基础

可加Q群：27896931

6. 关于小波db2和db4在c语言中的实现

小波db2和db4在c语言程序中具体差别体现在阶数不同，也就是存储单元占用量不同，同时信号的延时不同，延时单元的长度也不同

7. 用Mallat算法编写小波变换的C语言程序

www.google.cn

8. 有关C语言运行小波变换程序的时间问题。

我只听闻小波技术，从未试过哈，只做过傅里叶变换而已。。。但是，你可以分析下你输出的数据，正确与否，就用matlab的输出和C执行器的输出对比下，就知道靠不靠谱了啊。
如果嫌数据太多，可以找一些对比分析软件，如果实在无法全部分析，就弄几个特殊点，让程序就在那一些点输出对应值，对比matlab和C执行器的输出，多次换点，再对比，就知道靠不靠谱了。。。。你只是通过运算时间，怎么知道靠不靠谱呢，呵呵，要是你用多核DSP来算，只是一瞬间的事情，岂不是更不靠谱了。

9. 怎么把小波变换分成低频和高频部分，单独进行处理,高人指教，在线等，所有积分送上。。。

没具体做过，不过按照你的意图，应该就是使用mallat算法的离散小波变换（DWT），在matlab中有现成的函数很易实现。
如果你非要用C，那么需要C下解决以下几点：
1.信号处理中卷积的C程序实现；

2.当实现了卷积计算的C代码，你需要得到某个小波函数对应的高频（带通）和低通分解和重构滤波器组，这可以在matlab的wfilters函数和帮助中得到详细的滤波器系数；

3.以上两个问题都解决了，那么直接对信号用滤波器进行卷积再去掉结果越一半的数据即可得到DWT的高频细节和低频逼近系数，用这些系数补零插值再卷积一次就得到了高频细节和低频逼近的重构信号（这就是你要的高频和低频的两个单独部分，然后你想咋处理就咋处理）。对于多层的DWT你只要对小波系数逐层减半即可实现DWT每层尺度伸长一倍的效果。

4.这时你会碰到卷积滤波的边缘效应，你需要对原始信号进行扩展，你可以参看matlab中dwtmode函数的帮助文档；

5.对于二维DWT，你还要实现张量积C代码的实现，还有二维离散采样（行抽取，列抽取和斜抽取）的C代码（在dwt2函数的帮助文档中有详细的流程图）。

由于matlab函数通常是用C语言编制的，而且每个函数都可查询原始代码，如果你有一定的水平，可以试试。不过你这问题需要较为深入的小波分析，信号处理和数学理论的支持，能不能编出象样的C程序就不知道。这问题matlab很容易就解决了，可非要用C，这就是一个比较大的工程了，一个人要花费好长时间还不一定编的成，光数学理论就伤不起啊，技术问题就更别提了！到现在为止也没见有啥象样的小波C程序，你搜搜试试吧。能想到的就这么多了，水平有限，仅供参考，Good luck!

10. 如何用C语言实现小波多层变换wavedec2

C语言标准库当中没有这种函数，去OpenCV库当中找找有没有类似函数，如果自己实现想是相当麻烦