小波演算法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庫當中找找有沒有類似函數，如果自己實現想是相當麻煩