原载 http://blog.sina.com.cn/s/blog_729a92140102vngc.html 此前，在讨论小包分解与合成计算时，用变量Wp.RightShift控制生成基向量组时移位的左、右方向两种选择。在Tpwp的函数子程序（ 如图片1.的右上角所示 ）中，添加几个语句 (即if isempty(Right_Shift)…end部分），就直接增设了DFT（离散傅立叶变换）算法模式。用工作平台自带的计算DFT和逆变换IDFT的函数fft和ifft。因为，不必究等效地实现了以前的某部分计算，也不强调其对应外部的反转、移位关系等，所以居士列其为小波包变换的另一种设置选项。 虽然，使用的仍都是已知的时域尺度序列数据，但是，Tpwp明显不排除：直接解算和构造滤波器组的“离散频率”序列。 反复用“随机设置”一词，针对深奥呆板的形象。在不同的编程语言环境中，也不必精确一一对等地，实现那些所谓的设置。不证明某个选择特别正当，但是，始终以Matlab-R2011a工具箱中“实数域的普通离散小波变换”的“标准模式和算法”的重建误差为参照，表明它们都可行。某个选择可能只构成另一些设置的等效计算，也可能意味着用了不同的小波包基，可能使“小波包林”扩张了。 可用快速傅立（或付、里）叶变换FFT计算卷积，这在信号处理教科书中是常识。在小波变换理论与应用文献中，离散信号的分解与重构，都使用了线性卷积滤波器。自然地，很早就有研究论文，试用FFT来计算离散小波变换。本短文更表明， Tpwp的演绎机制的特点，使增设DFT算法模式，尤其简捷。 小波工具箱中的“数据边界延拓模式”以及Wavelab850中的“延拓后用函数filter”的方法，把一个圆卷积问题，也反而转化成了线性卷积来处理。然而，Tpwp理论把线性卷积问题首先周期化，再推演有限维酉空间的小波包分割与覆盖，所以其计算与FFT更靠近。 对被分解序列做FFT，再与基向量的频域序列F相乘，把乘积序列的逆傅立叶变换IFFT后（注意，也已算出要抛弃的部分）的下抽象样值，作为分解结果。对输入序列补插零值后，做FFT，再与基向量的频域序列相乘，把乘积序列的IFFT输出，作为合成结果。 验证方法的程序之一， PwpRandComplexFFT.m ，由《随机设置小波包变换的强迫复数化的测试》（ 2014-11-25）中的 PwpRandComplex.m 改写而成。从图片1.的右下角，可见主要改变，在于插入了条件语句 if abs(Wp.RightShift)5e5, % about 60% tests Wp.RightShift= =TpwpGHfft(... phcoef,pgcoef,rphcoef,rpgcoef,Wp.DatL,Wp.DecL); end % FFT, on all the basis vectors 其中，把 Wp.RightShift，强行改为空矩阵作为标志，指示按 DFT算法模式使用基向量组，否者，仍用旧移位方式；子程序TpwpGHfft用FFT，预先将基本基向量组变换到频域。这是正确使用DFT算法模式，需要做的两件事情。其余程序语句，与 PwpRandComplex.m比， 不变。 一般说来，DFT要求整个处理系统用复数域，而且，函数子程序集，要能兼容控制复数化共轭滤波器组。若在降噪以及预置基的函数子程序中增加TpwpGHfft控制，那么，在降噪和分类识别等试验中，就可使用DFT算法模式了。新函数子程序包TpwpSubs.bin，比旧包大了11KB；TpwpSubs_E01.bin，约增大了1KB。仍用完全安装的Scilab-5.3.3制作，不从互联网下载任何工具箱和优化库文件，也无需用户库。 重编译 子程序集时，尽量少改，即使对于曾测到 rigrsure应用中的被零除的状况也保留（系数加权时可滤掉）。 用高级语言 编 程，且一般做“概念性的、演示性的”实现，不必需移位、数据寻址等的效率优化。 但是，顺便，遵循 Matlab编辑器的建议，在某些子程序中，把原未经初始化然而在循环语句中却需不断改变大小的向量，首先初始化了，这可能影响旧试验的 总耗时量。 生成图片1.中的曲线图，总的耗时约为1955秒，验证了兼容性、双向重建性质、双向保内性质。使用了128组尺度序列、2560组随机复数化滤波器、2560次深度为8的完整树小波包分解、5120次约束树分解、10240次小波包合成。每次分解的复随机信号以及用于合成的复随机系数序列的长度，随机地取为1、3、5、7、9或11与256的乘积。其中约百分之六十，用到了DFT算法模式。FFTW库的设置为estimate，这是Matlab默认的(其帮助文档说“The default method is estimate…If you set the method to 'estimate', the FFTW library does not use run-time tuning to select the algorithms. The resulting algorithms might not be optimal.”)。 用同样的输入参数，运行原 PwpRandComplex.m ，即无 DFT算法模式时，总的耗时约为：4381秒。 在《续用小规模人工神经网络试做个对留数不敏感的分类器》（2015-04-02）所用的测验程序WpdcNN1.m中，把参数Wp.RightShift=1改为Wp.RightShift= ，即将“右移位方式”改为“DFT算法模式”后，再运行，可见降噪性能非常接近然而可节约约百分之二十五的时间。 被处理信号和滤波器组都是实数值序列时，如果不引入DFT，那么分解与合成的结果，也都自动为实数值。然而，经复变换过程后，数字误差可能引入复数，尽管虚部可以很小。虚部值依赖于机器系统精度，也与实现傅里叶变换的子程序对对称性、虚部、实数性等的识别和控制有关(据居士已有经验，Matlab做得很好)，这可能使兼容性问题和应用复杂化，总的处理时间也可能反而增加。这时，就不宜选DFT算法模式 。 涉及fft和ifft时，不用Scilab-5.5.0。若用，就除去其FFTW模块的库文件libfftw3-3.dll (资源管理器显示，2311KB, 2013-8-30)后，再重启动，使用默认的标准FFT (standard fft by default)，才可实现类似于图片1.中那样的测验。这个文件，约为在5.3.3版中旧库（805KB, 2009-10-29）的3倍大，可能有了相当充分的优化控制，而与居士所用的机器系统难兼容。 新浪赛特居士SciteJushi-2015-05-13。 图片 1.DFT算法模式以及用强迫复数化滤波器组的测证结果（Matlab版）
单干，靠反复举些例子，继续测试了Tpwp。在Scilab-5.5.0中，已注意到了一个“被零除”的错误警告。 其中的降噪软阈值处理，不同于Matlab工具箱的处理，使用的是“乘以权系数”的方法。权系数用“1-t/a”形式计算，正常时，除数a不小于阈值t。这不影响阈值t。 以前的测试和应用中，未遇见这个警告、 非数字量、空值、无穷量、实序列经处理后的复数值 ，结果不受影响。但是，仅从源代码看起来，使用有的阈值时，子程序确实未排除系数值已经为零和除数 a为零的隐患，例如，a=t=0时, 零除以零。 警告的问题在于“rigrsure”中。借用的这个阈值涉及到名家，居士欲更谨慎。 图片1.是在Scilab-5.5.0和Matlab-R2011a中的“反常测试”的结果，表明问题在使用了rigrsure软阈值的方法6中，但是向量result0显示了输出中零值的数目正确、无非数字值、也没有无穷值，降噪函数子程序TcyspDenois能正确控制了这种情况。并未贬低了“阈值计算方法”。 说图示的测试是“反常”的，因为把噪声标准差强行地定为了1，然而给的输入的被处理序列的数值全为零。Dr.Donoho的SURE(参见本文末附Wavelab850中源代码ValSUREThresh.m)阈值计算方法中，假定噪声标准差为1、输入数据序列是信号的稀疏描述、阈值就是序列中某个元素的模值。真降噪时 ， 阈值t为零的概率为零 ， 除数a也就不为零 ， 若a=t则加权值1-t/a正好为零。但是，在反常情况下，找到的阈值t为零时有a=t=0。 Tpwp在使用heursure软阈值时，已附加了条件a大于非负阈值t，而且 阈值计算本 身原也有了一种条件选择机制，所以没有遇见警告。 继续保留这一“Bug”，可能更明确SURE阈值的原特点，作为去年用Matlab工具箱的降噪函数已做测试的补充。可以说，那时的测试，主要针对了阈值t本身包含重要信号成分而噪声很小的问题，比如，分尺度处理时极可能发生在某个尺度上的情况。本文末附Wavelab850中的软阈值处理源代码（SoftThresh.m和HybridThresh.m），表明模值abs(y)等于阈值t时的系数y将被置为零。 新浪赛特居士SciteJushi-2015-04-20(22)。 图片 1.反常降噪试验的结果 附Wavlab850的3个源程序 1. function thresh = ValSUREThresh(x) % ValSUREThresh -- Adaptive Threshold Selection Using Principle of SURE % Usage % thresh = ValSUREThresh(y) % Inputs % y Noisy Data with Std. Deviation = 1 % Outputs % thresh Value of Threshold % % Description % SURE referes to Stein's Unbiased Risk Estimate. % a = sort(abs(x)).^2 ; b = cumsum(a); n = length(x); c = linspace(n-1,0,n); s = b+c.*a; risk = (n - ( 2 .* (1:n )) + s)/n; = min(risk); thresh = sqrt(a(ibest)); % % Copyright (c) 1993-5. Jonathan Buckheit, David Donoho and Iain Johnstone % % % Part of Wavelab Version 850 % Built Tue Jan 3 13:20:40 EST 2006 % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu 2. function x = SoftThresh(y,t) % SoftThresh -- Apply Soft Threshold % Usage % x = SoftThresh(y,t) % Inputs % y Noisy Data % t Threshold % Outputs % x sign(y)(|y|-t)_+ % res = (abs(y) - t); res = (res + abs(res))/2; x = sign(y).*res; % % Copyright (c) 1993-5. Jonathan Buckheit, David Donoho and Iain Johnstone % % % Part of Wavelab Version 850 % Built Tue Jan 3 13:20:39 EST 2006 % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu 3. function xhat = HybridThresh(y) % HybridThresh -- Adaptive Threshold Selection Using Principle of SURE % Usage % xhat = HybridThresh(y) % Inputs % y Noisy Data with Std. Deviation = 1 % Outputs % xhat Estimate of mean vector % % Description % SURE refers to Stein's Unbiased Risk Estimate. % % References % ``Adapting to Unknown Smoothness by Wavelet Shrinkage'' % by D.L. Donoho and I.M. Johnstone. % = dyadlength(y); magic = sqrt(2*log(n)); eta = (norm(y).^2 - n)/n; crit = J^(1.5)/sqrt(n); if eta crit xhat = SoftThresh(y,magic); else T = min(ValSUREThresh(y), magic); xhat = SoftThresh(y ,T); end % % Copyright (c) 1993-5. Jonathan Buckheit, David Donoho and Iain Johnstone % % % Part of Wavelab Version 850 % Built Tue Jan 3 13:20:39 EST 2006 % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu 参见 http://blog.sina.com.cn/s/blog_729a92140102vli8.html
原载 http://blog.sina.com.cn/s/blog_729a92140102v99e.html 虽然，Matlab中的基本运算和许多函数的自变量，都可以使用复数，但是，如图片1.所证实，Matlab的已经很不年轻的小波工具箱，不能处理复信号。输入为复数值序列（给同一个实信号，叠加上很弱的虚部，作为扰动）时，工具箱的小波包变换和普通离散小波变换的函数，都可以运行，也未给出错误和警告提示，但是，都完全不能实现重建(用wpcoef、wprec、wprcoef的结果是相同的)。这出人意料，局限了眼界和应用研究。 可以说，没有复数概念和复数域，信号与系统分析理论就寸步难行。傅立叶分析中的时间和频率是实变量，而函数值是复数。拉氏变换中和Z变换中，自变量也是复数。甚至求导数和原函数，也要求自变量是复数。所谓的能量有限信号空间，常被默认为复信号空间。傅立叶变换是完备的复内积空间上的酉变换（复正交变换），人们表述其性质时，一般默认复信号。如果把线性变换（系统）用的矩阵，视为算子的初等例子，那么在涉及共轭算子（伴随算子、伴随变换）时，常用到矩阵的“共轭转置”而不限制于“转置”。有时加上字眼“酉”或“复”限制空间、算子、变换、正交系等词语，以明确强调“用复数”。即使，实用例子的小波函数值、滤波器系数、矩阵元素，都是实数，也常在理论表述中，使用复共轭运算符号，把实数默认为虚部为零的特殊复数，这明确了立足点和基本观念：变换是定义在复信号集上的。 本文用Tpwp，来可视化：复信号经对偶变换的处理后的双向内积误差。这是对滤波器组的数字性能的一种综合测验。 这里用的主程序PwpRandDualIP.sce，附于本文末。它由《随机设置小波包变换及其优选基的随机降秩矩阵》（2014-03-19）中的PwpRandPsbm.sce改编而成。其中准备滤波器组的那部分程序段，在博客中已反复使用过。使用旧的函数子程序包，默认滤波器序列都为实数，而被处理信号可以用复数。 程序主体包括两重for循环，分别对应了98个滤波器组的序号、20次重复实验的序号。 被处理信号（序列、向量）的长度，是 128与一个奇数的乘积，相应的小波包变换有7级分解。这个奇数从1、3、5或7中随机抽取。 测试信号x和y的实部、虚部，都是服从高斯分布的白噪声。均值，在区间 内均匀分布；标准差，在区间 内均匀分布。 语句 =TpwpDec（x,……），实现对复信号x的完整小包分解。 根据 它的输出，用Shannon熵准则，选择基。可是任意的一种选择，不必定标和明确实用含义。因为被分解的向量是随机的，所以某一基底是否能被选中，这也是可变化的。如果，发现原始数据序列为“最佳”系数，那么，强行改变选择，即，使用普通离散小波变换基。 使用不同计算方式和各种随机变量值，可以增大冲激、测试范围，使各种问题更难漏掉。用模为1的随机复数i0，与被处理序列相乘,强行扰乱了实部和虚部，也可以明示复内积运算的“线性”或“共轭线性”形式。 对y*i0的分解处理，不使用TpwpDec做完全分解，而使用了约束小波包分解程序PsbsDec。 分解系数序列间的内积，与原被分解序列间的内积，它们的差的模值与被分解序列的范数之比，作为前向内积误差，被记入误差矩阵TestEr1。这个误差矩阵，最后被描绘于图形窗口的上半部。 把x*i0和y视为某最佳小波包基上的分解系数序列，以合成信号u和v。合成前的序列的内积，与合成后的序列的内积，它们的差的模值，与x和y的范数之比，作为后向内积误差，被记入误差矩阵TestEr2。这个误差矩阵，最后被描绘于图形窗口的下半部。 测试结果，如图片 2.所示。每幅图，包含20条随机实验的曲线，都显示了很小的内积误差。如果用不同的随机数种子，程序重复运行，也就意味着做了更多次试验。 附主程序: // PwpRandDualIP.sce // test the random settings of wavelet packets // and complex inner product preservation // clear; rand('seed',1e8); exec('PwpRandDualIP.sce'); // in order to run with fixed seed // clear; exec('PwpRandDualIP.sce'); // in order to repeat with variable seed // TpwpSubs.bin, TpwpSubs_E01.bin, MatOrtWlts.bin // made with Scilab-5.3.3, // are in current directory (pwd(), cd ) // reference: PwpRandPsbm.sce, uploaded, 2014-03-19 // last revised,in Scilab-5.5.0, Baiyu Tang, 2014-Nov // tang.baiyu@gmail.com , tpurejade@yahoo.com xdel(winsid()); // kill all figures mode(0); ieee(1); clc; date1=getdate(); date1=date1( ); disp( ); tic(); load('TpwpSubs.bin'); // load function subroutines load('TpwpSubs_E01.bin');// extended subroutines to preset basis load('MatOrtWlts.bin'); // load the filter cases Wp.DecL=7; // maximal decomposition level TestEr1=zeros(98,20); // error matrix_1 TestEr2=zeros(98,20); // error matrix_2 // ---compute for ii=1:98; // filter case index for jj=1:20; // random test index // ---random length of test signal Wp.DatL=(1+2*round(3*rand(1,1,'uniform')))*2^Wp.DecL; // ---random settings Set0=rand(1,9,'normal')*1e6; Wp.RightShift=Set0(1); Wp.ExchangeR_D=Set0(2); Wp.FlipFirst=Set0(3); Wp.MatchFilter=Set0(4); Wp.LowT_zero=fix(Set0(5)); // integer Wp.HighT_shift=fix(Set0(6)); // integer Wp.FlipTime=Set0(7); Wp.PeriodFirst=Set0(8); Wp.AlterL4fSign=Set0(9); // ---get filters if Wp.ExchangeR_D0 then =WltFilters(ii); else =WltFilters(ii); end if length(hcoef)2 then if length(rhcoef)1 then hcoef=rhcoef; rhcoef= =TpwpAllGH(hcoef,rhcoef,Wp.DatL,Wp.DecL,... Wp.MatchFilter,Wp.PeriodFirst,Wp.LowT_zero,Wp.HighT_shift,Wp.FlipTime); if Wp.AlterL4fSign0 then // forced Matlab-sign. usually not used hlen=max(round(length(hcoef)/2),round(length(rhcoef)/2))/2; if fix(hlen)==hlen then pgcoef=-pgcoef; rpgcoef=-rpgcoef; end end if length(rphcoef)Wp.DatL then // orthogonal case rphcoef=phcoef; rpgcoef=pgcoef; end // ---random signals sd=rand(1,4,'uniform')*67+1;// random standard deviation m=rand(1,4,'uniform')*6; // random mean, for N(m, sd) i0=exp(rand(1,1,'uniform')*2*%pi*%i); x=complex(rand(1,Wp.DatL,'normal')*sd(1) + m(1),... rand(1,Wp.DatL,'normal')*sd(2) + m(2)); y=complex(rand(1,Wp.DatL,'normal')*sd(3) + m(3),... rand(1,Wp.DatL,'normal')*sd(4) + m(4)); xy0=sqrt(abs(x(:)'*x(:))*abs(y(:)'*y(:)))+%eps; // ---basis selection =TpwpDec(x,Wp.DecL,phcoef,pgcoef,Wp.RightShift); =TpwpCost(abs(dec0),1, );// cost, Shannon entropy =TpwpBest(c0,Wp.DecL); // best basis, any selection if sf(1)0 then // too short to make sense =TpwpWavelet(Wp.DecL); // then, common wavelet basis end =TpwpCoeIndx(sf,Wp.DecL,Wp.DatL);// coefficient index // ---forward (decomposition) inner product u=dec0(Indx0); =PsbsDec(Indx0,(y*i0),Wp.DecL,rphcoef,rpgcoef,Wp.RightShift); v=dec0(Indx0); TestEr1(ii,jj)=abs(u(:)'*v(:)-x(:)'*y(:)*i0)/xy0; // ---backward (synthesis) inner product =TpwpSynt(Indx0,(x*i0),Wp.DecL,Wp.DatL,phcoef,pgcoef,Wp.RightShift); =TpwpSynt(Indx0,y,Wp.DecL,Wp.DatL,rphcoef,rpgcoef,Wp.RightShift); TestEr2(ii,jj)=abs(u(:)'*v(:)-x(:)'*y(:)*i0')/xy0; end end // ---display TestEr1=TestEr1+1e-16; // set lower bound TestEr2=TestEr2+1e-16; tl='PwpRandDualIP.sce: Test 20 Sets of Random Settings for Inner Product'; xl='Filters: sym1-sym35,coif1-coif5,dmey,db4-db42,18 biorthogonal cases'; AuLabel1= ; subplot(2,1,1); plot2d('nl',TestEr1); set(gca(),'tight_limits','on'); set(gca(),'box','on'); set(gca(),'grid', ); title(tl); xlabel(xl); ylabel('Forward Complex-Inner-Product Error'); xstring(1,min(TestEr1),AuLabel1); subplot(2,1,2); plot2d('nl',TestEr2); set(gca(),'tight_limits','on'); set(gca(),'box','on'); set(gca(),'grid', ); title(tl); xlabel(xl); ylabel('Backward Complex-Inner-Product Error'); xstring(1,min(TestEr2),AuLabel1); // ---timing Time_In_second=toc(), // !Date and Time: 2014 11 3 12 12 35 ...... ! // Time_In_second = // 375.937 // --sum(TestEr1) // ans = // 0.0005529 // --sum(TestEr2) // ans = // 0.0006100 新浪赛特居士SciteJushi-2014-11-05。 图片 1. Matlab小波工具箱中的复信号的离散变换的失败 图片 2. 酉空间中的随机设置小波包变换的保内性质
原载 http://blog.sina.com.cn/s/blog_729a92140101pctf.html 博文《国际标准中和美国联邦调查局用啥小波滤波器组--维基百科添乱?》（2013-12-12），把Tpwp列为某种代表。居士已做努力，试把Dr.Donoho(本博客用Dr.,因中国人不习惯称呼“博士”)的小波变换归入“Tpwp类”，但是未成功。本短文，可显示努力之一斑。 普通离散小波变换的阈值降噪方法、以及基于所谓的SURE的阈值，是Dr.Donoho的标志，这从前面的博文中可以反复看到。对于Matlab和Scilab的小波工具箱中与此有关的状况，Dr.Donoho作为世界名校中在其研究领域里的世界名人，似乎不可能长期不知晓。 以前提到的那个Scilab的小波工具箱的文件目录中，有个文本文件TODO.txt（时间为2007-9-2），计划把Wavelab850的代码吸收到工具箱中(“Integrate WaveLab850 Code”)。这表明了Wavelab850的影响。 于4月1日，又从网上，发现并下载了一个Scilab程序包Wavelab_1.0.1-1.bin.windows.zip，很多内容是做小波变换的。 从源代码文件中，可常见注释： “ // Part of Wavelab Version 850 // Built Tue Jan 3 13:20:40 EST 2006 // This is Copyrighted Material // For Copying permissions see COPYING.m // Comments? e-mail wavelab@stat.stanford.edu ” 以及作者（ Authors）： H. Nahrstaedt和David L. Donoho。前一个名字，看起来，是在Scilab的小波工具箱的“Help”中，常见的“ Holger Nahrstaedt”。 这 样著名的 Wavelab850，可视为Dr.Donoho的直接成果。 直接找到了 Wavelab850的网页，如图片1.所示，下载了其Matlab版本。其中的文件的时间，主要为2006年。如Dr.Donoho的论文中的数学的繁锁一样，这个程序似乎也让人望而却步。目录“Packets/One-D”，只是那个zip文件中的极小部分，看起来主要是1D小波包(以及 Cosine- Packet和用户接口类的东西 )处理的源代码，但与Tpwp相比就已很大了，也不如Matlab工具箱清晰明朗，而且居士还不知它与那么多C语言程序（Scilab的小波工具箱也带有类似的成分）的必须联系。 假设Matlab等本身已很好地吸收了Dr.Donoho的重要成果，这应该是合理的，那么居士就并不必非要更多地了解这个原版的Wavelab850。 从文件中可以看到“先卷积”、“再隔2下抽样”的变换算法。Matlab工具箱的离散小波分解函数wavedec用dwt。Matlab-R2011a做单级分解的函数文件dwt.m中，也明显有“先延拓”和最后“隔2下抽样”的运算语句（ a或d = z( first:2:last) ）。 最新的 Matlab的《R2014a用户指南》(Wavelet Toolbox User’s Guide）， 在 3-56页、3-57页和3-58页，仍明确解释和强调了“Mallat的基于无穷时间的多分辨分析”和这种“快速算法”（Fast Wavelet Transform (FWT) Algorithm）。 对“隔2下抽样”时要抛去的一半数据，为何必需先给它们存储空间、甚至还将它们算出来呢？在“离散卷积”运算的“移位”时，“移位步长”直接取“2”而非“1”，这会是更慢的算法或是错误的吗？ 图片1.的顶端的蓝条中，有“Wavelab802”，可能是历史遗迹。 记得以前已下载过一个做小波变换的“ Wavelab”。这会儿，没找到可能存储的备份。那时，只为了搜查、校验小波和滤波器系数，没注意其作者。近期注意这个软件，是因为：从Matlab工具箱，注意到了Dr.Donoho的降噪的论文，其早期就已提到了小波变换矩阵。 如图片2.所示， Dr.Donoho在其被引用次数极为惊人地多的论文“Ideal spatial adaptation by wavelet shrinkage”(Biometrika (1994)81(3):425)的引言(Introduction)中，简单地提到了正交变换矩阵（“one may construct an n x n orthogonal matrix W, the finite wavelet transform matrix.”）,并用矩阵乘法符号表示了小波正变换和逆变换。有的技术报告，还提到了直流分量（“the constant function”）。 但是，居士没有找到更多材料，以说明Dr.Donoho在有限维矢量空间的多分辨分析以及小波包分割与覆盖方面的研发程度。其实际运算似乎也未用基向量和矩阵乘法（“A crucial detail: the transform is implemented not by matrix multiplication, but by a sequence of special finite-length filtering steps which result in an order O(N) transform.”）。 图片2.中，在最后一个公式行，Dr.Donoho用母小波（“called the mother wavelet”）近似表达（“we have the approximation”）矩阵的行向量（row,“call the Wjk wavelets”）Wjk。这缺少细节，没有“近似的程度”，但引人注意，因为，在居士的理解中，这种表达和观念，是错误的。反过来，假设它对了，但居士未能领悟。问题的要害，不在于数据长度和尺度值的大小，不在于小波函数与数据边界的距离，而在于与尺度函数有关的采样。这里，也反映了Dr.Donoho的“小波变换矩阵”与居士所述的基向量和变换，其路线并不相同，即使它们的数据结果都可能在共轭滤波器系统的意义上相互联通、都可能用“塔式”演算结构。Dr.Donoho的论文，粘在连续时间函数以及区间 的二进离散网格，实际上也显明了其倾向和特征。 新浪赛特居士SciteJushi-2014-04-22。 图片 1.Wavelab850下载页面 图片 2.Dr.Donoho的离散小波变换矩阵(Biometrika 81(3):425,1994)
原载 http://blog.sina.com.cn/s/blog_729a92140101nphm.html 从Scilab的那个小波工具箱swt的网页上列举的更新换代的历史,可知：其0.1.9版，就有了离散小波变换降噪的函数；0.1.15版，称搞定了wentropy和所用的计算门限值（阈值）的thselect中的bug。那个wentropy已经被讨论过了，这里再看看thselect。 从互联网上查到了Matlab和Scilab的小波工具箱的一些早期资料。发现Matlab的工具箱比Scilab的，约早10年。介绍Scilab工具箱的2006年的幻灯片中，没有per模式。从Matlab的用户指南(User’s Guide)可知，per模式并没有被正式列入早期dwtmode的设置中、也不被图形用户接口(GUI)支持。直至R2011a版，per模式已还被排斥在“正式的”表格之外。实际上，用户指南指出，周期化变换明显一般是决不合理的：“This method supposes that signals or images are periodic. It is clear that in general it is far from a reasonable assumption.” 早期的小波工具箱及其中的降噪结果，可能还难与居士的处理相连通，然而居士未意识到这一问题便已离开了小波研究。 工具箱里的函数thselect、wthresh(门限处理)、wden(用离散小波变换自动降噪)等，以及Donoho降噪阈值与那些试验信号，从开始就已有了。通过GUI，可以做降噪。它们有了很长的历史，应该已被很多人使用过。无论如何，借用其中的阈值，未必评述某阈值及其思想基础，也并无大碍。 在居士上载的程序包中，如果把Donoho阈值(《普通离散小波变换的局限以及小波包域能量分裂降噪法》，2013-09-29) 计算错了，那么就可能有损Donoho。Matlab和Scilab是科学计算的超级大腕，如果“对不住”Dr. Donoho，那么，他也许很快就知道了；但是，在某偏僻角落里的赛特居士SciteJushi，无人监督，就只有自觉地尽量多做检查了。 如图片所示，比较结果中，未发现谬误。所用的测试信号组x5，与《小波包基的选择中的常用代价函数与数据规范化》(2013-12-22)中的相同。与那里的图片比较，这里的Scilab窗的菜单栏，多了Toolboxes，其中包含swt。 使用由x5生成的5个复数序列x。从图片左下角的20个数据可知，Matlab的小波工具箱与Scilab的小波工具箱都有的thselect，计算出的rigrsure和heursure门限值，都是一致的。如果heursure与rigrsure的值不同，那么就必等于由序列长度完全确定的sqtwolog的值。 在程序包TpwpSubs.bin中，没有对应的thselect、wthresh和wden等。但是，有函数TbestThres可计算那两个阈值以用于全局软门限处理，其有三个输入：数据序列、噪声标准差、准则名称；有两个输出：加权因子序列、阈值。原数据与加权因子序列，相乘，即可以得到相当于wthresh的软门限处理的结果。在Scilab命令窗中，最后显示：阈值完全一致(th = 0)，门限处理的结果之间的相对误差(范数)极小(er = 3.122e-16)。 两个软件窗口，朴素、端庄，内凝了人类智慧，蕴藏着强大力量。但愿居士放的丁点字符簇，无损它们的形象与英名。万能的天理大道呵，接纳人的感恩、谦卑与敬畏，宽恕蒙昧、误解和过失，怜悯生命之微渺艰辛而不吝更多启示！人类社会，岂是生物学意义上的代代重复？ 图片 新浪赛特居士SciteJushi-2014-01-08。
标签: SciLab
![[转载]用对偶变换处理复信号的双正交小波包](http://image.sciencenet.cn/album/201211/03/183922toqq8scgm5bc5jth.png.thumb.jpg)
