有趣的是,《维基百科》(en.wikipedia.org)中“Daubechieswavelet”的“Implementation”部分,也把某种Periodization与其它模式(被称为高级(尖端,sophisticated)的方法)“相对”而论,而且强调只影响“very ends”。但是,居士对整个滤波器系统的周期化的测试,并不管支集影响有多远;在256点数据的8级分解中,最后一个近似系数,实际上,涉及到了所有的256个数据点。
在Matlab中,做些测试,是应该的!
Matlab系统默认的模式是sym,它在表格中排在第一位。此外,系统有函数计算“合理的最大”分解深度。在命令窗中,运行
wmaxlev(256,'db33'), wmaxlev(64,'db20'),
得到的最大分解深度为1、0,即是说:只该分解一级、或干脆就不能做小波变换。分解深度越小,由滤波器系数本身的近似性和计算中的数字截断所引入的误差(不是理论上的重建问题!),一般应该更小。所以,测试时,用sym模式、系统计算的合理分解深度,作为比较标准。
本文后的图片所示(右边)的结果,用单位能量的正态分布的随机序列(左边),比较了Matlab的小波变换中的3种模式、它的65组滤波器(依次为Daubechies1至Daubechies 45、 Coiflets 1至Coiflets5、15组双正交)。模式sym、zpd、per的误差(未做任何变换域数据修正)结果,几乎完全重合,而且,也与居士在Scilab中周期化滤波器系统以分解单位矩阵的结果一致。用模式zpd、per时,分解深度都被强行设置为8;然而,用sym模式时的平均分解深度,仅约为3,不到8的一半,可以说sym模式很占便宜。
下面是居士写的计算图片中结果的Matlab程序temp8x.m:
function [Test8x]=temp8x
%test Matlab-Wavelet-Toolbox: dwtmode and itsFilters
%in Matlab-R2011a. 2013-August.
% tic;[Test8x]=temp8x;toc,
%------prepare
close all force;
biof1=['1.1'; '1.3'; '1.5'; '2.2'; '2.4'; '2.6';'2.8';...
'3.1'; '3.3'; '3.5'; '3.7'; '3.9'; '4.4'; '5.5'; '6.8'];
dwtmod1=['sym';'zpd';'per']; %3 modes
rng(1.0e9); x=randn(1,256); x=x/norm(x(); %test data
ModeSymLevel=0;
Test8x=zeros(65,3);
%------compute
fori1=1:65; % 65filters
if(i1>0)&&(i1<46)
FilterNam1=['db'int2str(i1)]; %Daubechies
elseif(i1>45)&&(i1<51)
FilterNam1=['coif' int2str(i1-45)]; %Coiflets
elseif(i1>50)&&(i1<66)
FilterNam1=['bior' biof1(i1-50,];%Biorthogonal
else
error('Undefined Filter.'); %unknown
end
forj1=1:3; %3mode error, from dwt and idwt
if j1==1
dl0=wmaxlev(length(x),FilterNam1);%level,only for Mode sym
ModeSymLevel=ModeSymLevel+dl0; %memorize the levels
else
dl0=8;
end
dwtmode(dwtmod1(j1,); %set mode
[xc1,xl1]=wavedec(x,dl0,FilterNam1);%transform
xr1=waverec(xc1,xl1,FilterNam1); %inverse
Test8x(i1,j1)=norm(xr1-x); %norm of error
end
end
%------display
ModeSymLevel=ModeSymLevel/length(Test8x(:,1)), %'sym' average level
subplot(1,2,1); plot(x); axis tight; %display raw data
title('The Normalized randn(1,256) forwavedec');
subplot(1,2,2); %display the result
semilogy(Test8x(:,1),'-g'); holdon;
semilogy(Test8x(:,2),'ob'); semilogy(Test8x(:,3),'+r'); hold off;
grid on; axis tight; legend(dwtmod1,'Location','NorthWest');
title('temp8x.m: Matlab dwt Modes and ItsFilters');
xlabel('db1-db45,coif1-coif5,15 BiorthogonalFilters');
ylabel('Reconstruction-Error Norm');
end
欢迎光临 电子技术论坛_中国专业的电子工程师学习交流社区-中电网技术论坛 (http://bbs.eccn.com/) | Powered by Discuz! 7.0.0 |