单干，靠反复举些例子，继续测试了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_729a92140102uwlk.html 进入小波工具箱后，“SURE”,是令人不舒服的东西之一。某种阈值处理，而已，若属居士写的东西，不被喜欢就拉倒。但它与“名校名家”、“无偏”、“自适应”、“最佳”相关联，却易使人面临混乱。 《 小波工具箱中 sure不比shannon更危险吗 》 (2013-09-16)称： 可以肯定，现在 Matlab的主要问题不在于“Shannon”本身。后来《改善Matlab的小波包处理的一个简单例子》（2013-12-14）已表明：Shannon熵本身，是可以使用的，人们常用它是不奇怪的，虽然Wavelab的作者Dr.Donoho等质疑了其合理性、认为它最差，也许它确实有不足。 那两个熵准则的Matlab的定义式中，都有对“未取模值”的数据的直接取平方的运算，表明未准备复信号的处理。当然，把复信号的实部和虚部分开处理，是一种方法，但是，小波理论，实际上一般使用复数域、复信号空间。早在Matlab小波工具箱诞生之前，小波包基系数的Shannon熵的计算表达式，就常有数据取模和能量归一化运算。然而，工具箱的文档资料中，一直没有采用这些。 小波包变换用的那个SURE阈值的计算形式，不如普通离散小波变换降噪中的rigrsure似的繁琐，而类似sqtwolog的计算。虽然，普通的小波变换，只是小波包变换树中的一个局部，但是两个阈值计算相似而不相同。这里，反映了Dr.Donoho的小波包降噪的一个基本特点，它与Tpwp循着不同的思路（不论是非）。 按Matlab-R2011a的帮助文档所述，把它用于一个可与“Shannon”并列的最佳小波包基的选择准则，这很抬举了，所以也使其不得不面临信号定标和作为代价函数的那种一般问题的拷问。 现在宽容退步。依据用户手册，可认为它是只为降噪而设计的。即使如此，有同样问题：实际噪声标准差，和要求输入的阈值参数，都乘以同一个适当的常数后，类似《用Matlab的小波包变换试验能量阈值降噪方法》（2014-06-26）中那样，结果会怎样？ 可能试一试，但是注意到，最早和最新的手册都称：在白噪声的方差为“1”时，其工具箱的小波包降噪的性能良好（For now, suffice it to say that this method works well if your signal is normalized in such a way that the data fit the model x ( t ) = f ( t ) + e ( t ), where e ( t ) is a Gaussian white noise with zero mean and unit variance.）。 所以，直接再厚道些，就使试验中的实际噪声标准差只为1！要求用户对被处理数据用噪声定标，也不一定过分。当然，这已略使人不自在，因为，既如此，如果有了被处理信号时，阈值实际也就定了，那么，真正的SURE准则，实际无需其它“输入参数”。 干脆，由工具箱的函数ddencmp，自动确定降噪参数。 依据手册和帮助文档的长期的内容重复，可以假定：SURE、不寻常的代价准则、Donoho、自动设置参数的函数ddencmp、小波包降噪和压缩的函数wpdencmp、试验信号生成函数wnoise、常见数据长度、滤波器选择，在小波工具箱中应该已经磨合为一体了，可以构成小波包降噪的“金标准”。 从可选基的数目、对DFT的频率分辨率的逼近方面看，在累积误差、存储空间、时间开销都允许时，分解越深的小波包处理，对不同信号的自适应性就越强。但是，在与Dr.Donoho有关的小波变换和小波包变换降噪方法中，常见：选择（人为）适当的分解深度参数。 为保险起见，把好人做到底：做“金标准”时，用循环语句试验不同深度分解的处理结果，可看其中最好的。 设置一个“无所谓最不最佳的伪SURE”，当作灯泡：依据《 小波工具箱举例不当与 Donoho阈值》（2013-09-26），用“3倍标准差”作为硬阈值，强行替换自动的SURE的默认参数。 试验程序，如图片1.所示。其主体有4层嵌套for循环，从外到内，依次对应：由wnoise给出的6种基本测试信号，无噪信号的6种强度，10次(变量Te0)随机试验，不同的分解深度。 误差数据阵列Ea的三维寻址指标，依次为：无噪信号的编号，信噪比条件的编号，分解深度和准则参数的编号。 显示时，以整数1或0开始的行中，这第一个数字，“1”表示直接用ddencmp给的参数；“0”表示以“3倍标准差”做“硬阈值”的“伪SURE”。同一行中，第二个整数，表示分解深度；第三、四、五、六个数，分别表示，紧随其后的6个行显示的误差（范数比）矩阵中的元素的均值、标准差、最小值、最大值。 保存这一行数据。在程序结束时，以“IsDefault_Depth_Mean_Std_Min_Max”为名，再集中显示它们，如图片2所示。 此前的一行，显示一下默认参数，确认“SURE”。 从均值 、标准差、最大值看来，默认处理的结果，明显最差；“伪 SURE”的结果，最好；《用Matlab的小波包变换试验能量阈值降噪方法》的结果，处于二者之间。从可以实现的很低信噪比时的最小的误差值看，那个rigrsure最差，这里的默认处理中，存在比其它方法都更好的少数例子。 新浪赛特居士SciteJushi-2014-06-30。 图片 1. 设置Matlab小波包降噪的“金标准”的程序段 图片 2. 图片1.中的程序运行的结果（末尾部分）
原载 http://blog.sina.com.cn/s/blog_729a92140101p6ob.html 在信噪比很高的实验条件下，降噪后的波形图可能看起都很“完美”，但是，实际上，失真和畸变可能已经引入了相当大的绝对误差、和信噪比下降。例如，即使在变换域不做任何处理，经某“正变换、逆变换对”的处理后，也可能造成了误差。 在低信噪比条件下，对降噪后的波形的主观评价，更可能相当不一致。看起来，更好（如，光滑，舒服感）的结果，实际上均方根误差却更大、信噪比较低；均方根误差减小时，而视觉效果反而很差，例如残余的成簇噪声可给人以变态、面目全非的印象。在陈述降噪结果时，专门找“坏的”波形例子，一般不必要，只找“好的”，则可能使人感到“吹牛”——某名人在谈到写自传时，表达了类似的意思。 单个实验的波形图，可作为初级示例，是易于被普通观众感知的直观表达，但是专业参考价值太低，不能反映成功的几率，更不能取代批量处理的统计结果。幸运的是，计算机技术的发达，使研究者可以更快地浏览和统计大量的实验数据样本，尽可能少地引入主观选择、主观判断，发现问题，解决问题，明确各方法的适用条件，步步高升。 实际上，信噪比（SNR）这个概念术语的使用中，“信号”和“噪声”的内涵也并未绝对统一。例如，直流成分的问题。 在Matlab-R2011a的小波工具箱中，用函数wnoise生成小波变换降噪的测试信号（Donoho与Johnstone），其用输入参数SQRT_SNR的平方表示SNR（ % = WNOISE(FUN,N,SQRT_SNR) returns the previous vector X % rescaled such that std(x) = SQRT_SNR. The returned vector XN % contains the same test vector X corrupted by an % additive Gaussian white noise N(0,1). Then XN has a % signal-to-noise ratio of (SQRT_SNR^2)）。 SQRT_SNR即是居士在博客中使用的NFSA。 在最新的Matlab的《R2014a用户指南》(“Wavelet ToolboxTM User’s Guide”）中，峰值信噪比“ Peak Signal to Noise Ratio ( PSNR )”，在5-67页，用像素的最大值的平方，表示图像信号的功率。用均方误差表示“噪声”功率，与居士用RMSE评估降噪的结果相一致。 该手册，在5-17页，谈到Cycle Spinning和普通离散小波变换降噪的经典函数wden时，使用的信噪比SNR定义式，把本该在分母上的表示序列长度的N，误放到了分子上。居士认为，其程序不可能这样计算，这属于“行文Bug”，是小问题。它不致理论障碍、和方法致命缺陷。但是，那里说的也正是wnoise的bumps信号，该信号有较大的直流成分，其SNR定义式，用序列值的模的平方和，表示信号能量，对直流的处理未明确交代，那么已不必就与过去的SQRT_SNR一致。 在以分贝（ dB）为单位时，使用MSE（均方误差、误差的模平方的均值）与使用RMSE(均方根误差、MSE的平方根)是一样的，使用“幅度”和使用“幅度的平方”表示无噪信号的强度是相同的，也都有0的对数log(0)的问题。 在前面讨论小波包时，关于最佳基的选择的代价函数，提到了一个习惯：E(0)=0，即0的熵（log(0)无意义）被约定为0。这在文献中很常见。问题在于，log(1)=0，这使1的熵也为0。为避免这个麻烦，数据常需要适当地规范化。为叙述方便，不妨称这种“冲突”为，“零熵模糊”。 如果， 与上面零熵的情形类似，令SNR式中的log(0)=0，那么也不能区别于log(1)=0。不妨称这种现象为，“零分贝模糊”。使用NFSA为横轴（dB）、负RMSE为纵轴(dB)，可图示降噪结果。一般说来纵轴可表示SNR的增量，但是，如果降噪处理前的NFSA（横坐标值）小于零，横、纵坐标值之和也不“显著”大于零，那么，即使纵坐标值非常大，这种结果的意义，已变得模糊。例如，直接令降噪结果信号恒为0值，可能得RMSE约等于NFSA，也就得增量为NFSA的相反数。宜避免用小于零的NFSA值。 在居士的降噪程序中设置了“无信号”警告，可能强调这一问题。然而，在大批量重复实验中，太多的警告烦人时，不得不先运行warning(‘off’)予以屏蔽。 纵坐标值为负，不管怎样，都表明“处理性能差”。计算引入的数字误差、门限处理时引入的损伤等，已起主导作用，超过了原有的白噪声的“危害”。在极高信噪比的试验中，也可能看到在程序中设置的“无噪声”警告。相应方法已不能降噪。 使用博文《结合移位(Cycle-Spinning)平均的小波包域降噪的程序》（2013-10-15）中的程序，TDen3Signals.sci，仅修改NFSA的赋值（80：20：180）和极少的图形标注字符，然后运行，可得结果如图片1.所示。用软门限方法或阈值偏大时，例如“Soft Power”，在高信噪比的实验条件下，更易出现信噪比不增加反而下降的情形。 在用Donoho阈值的处理中，不做“无噪声”和“无信号”警告。但并不允许随意而为，而且情况更复杂。在图片1.中，基于Donoho的“SURE”的结果，从一个信号到另一个信号之间，表现出了更剧烈的性能“跳变”。过去，在利用“先验基”和整顿“额外极点”的试验中，也可以察觉性能“跳变”。这些不是居士将Donoho阈值用于小波包的特有问题，而可能是Donoho处理方法固有的隐患。 Donoho阈值的计算，假设了噪声标准差是已知固定的（ %THSELECT Threshold selection for de-noising. % THR = THSELECT(X,TPTR) returns threshold X-adapted value % Threshold selection rules are based on the underlying % model y = f(t) + e where e is a white noise N(0,1).）。 数据序列的长度，卷入其中，未必显示了对 SNR的考虑结果。早期的Matlab小波工具箱，就已把Donoho的离散小波变换降噪，置于重要地位，但其变换的主体，并不支持per模式（此外：它也还不能直接做居士已离开的“连续小波变换”——这绝不意味着谁更好或更差，但最新的手册，已显示了丰富的连续变换的内容）。函数wden中，用输入参数SCAL对噪声重定标，不同尺上用不同的阈值，非per模式中的变换结果的数据序列的长度随滤波器和分解深度的不同而不同，这些都更使阈值计算难被“直观地”把握。 在Matlab命令窗中，运行： clear; format shortG; format compact dwtmode('zpd'); % common, in the Help tR={'rigrsure';'sqtwolog';'heursure';'minimaxi'}; init = 2055615866; Er0 = zeros(100,4); % error lv = wmaxlev(2^8, 'sym8');% allowed maximal level for tN = 1:4, % soft-threshold-rule index for sp = 1:100, % wnoise,signal 6,snr=(sp/2)^8 = wnoise(6, 8, (sp/2)^8, init); xd = wden(xn, tR{tN},'s','one', lv, 'sym8'); Er0(sp,tN) = norm(xd-x)/norm(xn-x); end % error norm ratio, in Er0 end % almost invariable norm(xn-x) disp(date); disp(max(Er0)); % display subplot(2,1,1);loglog(((1:100)/2).^8',Er0);% plot axis tight; grid on; xlabel('snr in =wnoise(6,8, snr, init)'); ylabel('norm(xd-x)/norm(xn-x)'); % RMSE ratio title('Error Increase in xd=wden(xn, ...) ?'); legend(tR,'Location','NorthWest'); subplot(2,1,2); plot(x); axis tight; title('Last x in =wnoise(...)'); 可得图片 2.的曲线图。 这里使用了Matlab文档资料中常见的降噪函数及其主要参量设置（如随机数种子、zpd模式、计算得分解深度）。所用的信号，其时间带宽积大，是与Donoho有关的测试信号中难以处理的一个例子。大部分情况下，Er0的几条曲线，都“在 1 以上”，表明了，“噪声”不减小，反而“被放大”。尤其引人注意的是，与“rigrsure”有关的“值超过2000的短暂突起”，以及相应的阈值“heursure”在阈值“rigrsure”与阈值“sqtwolog”之间的复杂切换。 把上段程序中的dwtmode（‘zpd’）改为dwtmode（‘per’），即用最小周期化模式；分解深度改为lv=8，即做最彻底的分解。再运行，得图片3.的结果。“突起”部分，更宽；更高，达几百万。 用前面提到的信号bumps, 用R2014a最新手册中的长度参数（10）、小波名称（sym4），且强行减小一级分解（10-1=9）。再实验，得图片4.，比图片3.的结果，还更糟十万倍。 新浪赛特居士SciteJushi-2014-04-12。 图片 1. 高信噪比条件下的Tpwp降噪使信噪比下降 图片 2. 小波变换（zpd,wmaxlev）降噪的性能对信噪比条件的依赖 图片 3. 小波变换（per,max level）降噪的性能对信噪比条件的依赖 图片 4. 小波变换（per,max-1 level）降噪的性能对信噪比条件的依赖
原载 http://blog.sina.com.cn/s/blog_729a92140101obm9.html 附件 Den6Comp.zip 在居士看来，以前的博文中提及的小波包降噪、参数估计等的子程序，已可以处理复(值)信号，但是，要处理复信号，涉及加噪声、RMSE和NFSA等的主程序还需要修改。 如图片1.所示，Scilab-5.3.3中的计算方差(variance,其帮助文档表明其输入可以为复数序列)和标准差(stdev, 其帮助文档指明输入为实数序列或矩阵)的函数的输出结果(format(‘v’,8),mode(0))，是复数。它们与Matlab-R2011a的std和var的结果(都为实数)，基于不同的概念。据居士所知，Matlab给出的这些结果，更符合通常的信号处理文献的表述，看起来与零均值的复高斯噪声的平均功率、功率谱密度的概念也一致。 对于复信号，计算RMSE值时，必须插入abs以取“差”的模值；而对于实(值)信号，可以直接求差值的平方。 试验所加噪声的方差，已知，为1。实部、虚部独立同分布，均分能量。与实信号时一样，处理结果相对于无噪信号的RMSE值的倒数(分子1，可视为处理前的RMSE、原加噪声的平均功率的平方根)，表示了信噪比的“增量”。 无噪信号被放大(NFSA, dB)前，其交流功率被设定为1，那么，NFSA值就是忽略直流功率时的降噪前的信噪比。但是，在功率归1化时，信号不能除以stdev输出的结果，而除以相对于其均值的RMSE。 在命令窗中，运行 clear; TestNum0=50; Gr2=0; RealOnly=0; exec('TDen6Comp.sci'); 可以得到 (约需2500秒)图片2.的结果。复制、粘贴后，在回车前，把Gr2的赋值改为Gr2=3，可以得到图片3.的结果。在warning(‘on’)状态，处理强阻尼信号时，可能看到一些“警告”提示:未分解的数据序列已被认为是最佳的；可事先用warning(‘off’),将“警告”关闭。 在回车前，把 RealOnly的赋值改为RealOnly =1，则强行限制只使用信号的实部做试验。这可以重复博文 《试试小波包变换降噪与复指数和信号模型拟合的串联》 (2014-01-22)的处理，但是可看到比在哪里更好的结果。尤其是低信噪比时用Donoho阈值的处理，变化更明显。产生差别的原因是，哪里在做线性拟合时未限定“阻尼”因子。结果变化越明显，那么，处理越可能具有某种不稳定性。这里上传的拟合子程序，强制调整单位圆外的极点，“掩饰”了这一问题，但是实用的程序，必需这种控制。 降噪的不少部件集成在了一起。简单地修改TDen6Comp.sci，可以进行若干试验。没有证据表明，各个部件或方面选择，都达到了最佳。实际上，某种“最佳”的定义，也未必在任何情形都合理、适用；技术环节的强强组合，其结果未必更佳，反而可能更失败。 程序包Den6Comp.zip上传在科学网附件中。主程序如下： // TDen6Comp.sci // wavelet packet denoising with cycle-spins,+State_space_SVD,+linear fitting // copy/paste following line modify value press Enter: // clear; TestNum0=50; Gr2=0; RealOnly=0; exec('TDen6Comp.sci');// Gr2=0 or 3 // TestNum0: total averages of tests // Gr2 =0: use first 3 signals; = 3: use other 3 signals // RealOnly = positive number: use the real part of complex signal vector; // otherwise, use complex signal // !! TDen6Comp.sci, TpwpSubs.bin, MatOrtWlts.bin, SsSvdLfit.bin, // must be in Scilab current directory !! // reference: TDen3Signals.sci, uploaded in 2013-10-15 mode(0); xdel(winsid()); // kill all figures ieee(1); date1=getdate(); date1=date1( ); disp( ); tic(); // timing function =TwltDenoise(Ordr,RawX00,WavNum,NFSAs,TestNum,PlotF,PlotP); // Wp: used settings // AveRMSE: result, RMSE; row: NFSA; column: method // Ordr: pole-model order // RawX00: original signal, SD( like std, RMSE )=1 // WavNum: filter case No. // NFSAs: values of NFSA for tests // TestNum: TestNum noisy tests per NFSA // PlotF: NFSA value for plotting waveform // PlotP: plot position // in Scilab-5.3.3,Baiyu Tang,2013-Dec // last revised,in Scilab-5.3.3,Baiyu Tang,2014-Feb // Tang.Baiyu@gmail.com // ---signal and wavelet settings--- Wp.WltNum=WavNum(1); // filter case No. Wp.DatL=length(RawX00); if Wp.DatL32 then error('Too Short Data Length.'); end // length of signal vector, not less than 32 Wp.DecL=0; tem=Wp.DatL; while tem==(fix(tem/2)*2), tem=tem/2; Wp.DecL=Wp.DecL+1; end if Wp.DecL2 then error('Odd Data Length.'); end // tree height // ---basis selection settings--- Wp.Thres=3; Wp.NormP=1.5; Wp.BestRule=5; // 5: l_1 norm // 4: log-energy // 6: l_p (Wp.NormP) norm // 3: threshold (Wp.Thres) entropy // otherwise: Shannon entropy Wp.FixWavelet=-1; // 0 : directly use wavelet decomposition // otherwise: find the best basis // ---biorthogonal specific settings--- Wp.ExchangeR_D=-1; // 0 : reversed filters h and h1 // otherwise: usual use Wp.AdjBioNorm= : no adjustment // 0 : unit analysis // 0 : unit synthesis // = 0 : equal Wp.DualYes=1; // 0 : dual transforms for selection. // otherwise: usual transform // ---common settings--- Wp.PeriodFirst=1; // 0 : G from periodized H // otherwise: first make g from h, then periodize them Wp.FlipFirst=-1; // 0 : flip h1 and h, before nothing else is done // otherwise: do not flip Wp.MatchFilter=-1; // 0 : make filter match first, equal length // and make time-zero at the first // otherwise: normal Wp.LowT_zero=0; // low-pass filter time-zero shift // arbitrary integer Wp.HighT_shift=0; // 2M+1 shifts, from h to g // arbitrary integer Wp.FlipTime=-1; // 0 : -t, like convolution // otherwise: like correlation // ! after periodization Wp.AlterL4fSign=1; // 0 : Alter the sign of periodized high-pass basis // prepared for simple compatibility with 'wavedec' of Matlab-R2011a // otherwise: do not alter. Wp.RightShift=1; // 0 : right shift to generate all basis vectors // otherwise: left shift // ---noise reduction settings--- Wp.MaxSpin=max(min(Wp.DecL*2, 16)-1, 0); // integer, 16, =0, maxima shifts, for cycle-spinning Wp.NoiseSD=1; // noise, std deviation // ======================== DenMode1= ; // threshold mode index, integers in ModeNam1= ; // threshold mode names AveRMSE=zeros(length(NFSAs),length(DenMode1)); // RMSE // =========label string======== AuLabel1= ; // ==============test========== Ir0=isreal(RawX00); // real or complex vector ? for kkk=1:length(NFSAs), SNR=NFSAs(kkk); // SNR in dB RawX=RawX00(:).'*10^(SNR/20); // amplify rand('seed',1.0e9); for iii=1:TestNum, if Ir0 then // real xn0=RawX+rand(1,Wp.DatL,'normal'); // a row vector else // complex xn0=RawX+complex(rand(1,Wp.DatL,'normal'),rand(1,Wp.DatL,'normal'))/sqrt(2); end =TcyspDenois(xn0,Wp,DenMode1); // xpwpd0: denoised vectors, one row means one mode // xn0 : noisy data // Wp : settings // DenMode1: thresholding mode No., vector. for jjj=1:size(xpwpd0,'r'), =PoleStSp(xpwpd0(jjj,:),Ordr); // pole estimate =PoleFitRes(pole1,xpwpd0(jjj,:)); // linear fitting RMSE=sqrt(mean(abs(yn0(:)-RawX(:)).^2)); AveRMSE(kkk,jjj)=AveRMSE(kkk,jjj)+RMSE; if (iii==1)(SNR==PlotF(1))(jjj==1) then // save yn0 for plot xpwpd0(jjj,:)=yn0(:).'; end end if (iii==1)(SNR==PlotF(1)) then // plot waveforms scf(10000); subplot(3,2,PlotP*2-1); if Ir0 then // real signal plot(xn0,'g'); plot(real(xpwpd0(1,:)),'r'); plot(RawX,'b'); xtitle( ,' ',' '); else // display the imaginary part of complex signal xn0=imag(xn0); plot(xn0,'g'); plot(imag(xpwpd0(1,:)),'r'); plot(imag(RawX),'b'); xtitle( ,' ',' '); end set(gca(),'tight_limits','on'); set(gca(),'box','on'); ypos=0.99*min(xn0)+0.01*max(xn0); xstring(1,ypos,AuLabel1); end end AveRMSE(kkk,:)=-20*log10(AveRMSE(kkk,:)/TestNum+%eps); // 1/average, to dB end scf(10000); subplot(3,2,PlotP*2); // plot the RMSE LinStyl= ; NumPlot=size(AveRMSE,'c'); for jjj=1:NumPlot(1), plot(NFSAs,AveRMSE(:,jjj),LinStyl(DenMode1(jjj))); end if PlotP==3 then legend(ModeNam1(DenMode1),-3,%f); end set(gca(),'tight_limits','on'); set(gca(),'box','on'); set(gca(),'grid', ); xtitle( ,... 'NFSA (dB): Noise-Free Signal Amplification',... '1/RMSE (dB) of Final Result'); ypos=0.99*min(AveRMSE)+0.01*max(AveRMSE); xstring(min(NFSAs),ypos,AuLabel1); endfunction // ********************* load('TpwpSubs.bin'); // load all the functions for WP processing load('MatOrtWlts.bin'); // load filter cases // ***** orthogonal filter cases,in MatOrtWlts.bin, // mostly extracted from Matlab-R2011a: // 1-35, Symlets 1 ... Symlets 35 // 36-40, 35+ ..., Coiflets 1,2,3,4,5 // 41, Meyer // 42-80, 38+ ..., Daubechies 4 ... Daubechies 42 // 81-98, 80+ ..., biorthogonal cases. // =WltFilters(filter-case-No), // one empty: forward/inverse transforms use the same basis vectors. load('SsSvdLfit.bin'); // signals,functions for pole estimate,linear fitting =PoleSignal(); // get 6 test signals format('v',8); for j1=1:3, // 3 signals per group WltNum=96; // 96, equivalent: CDF9/7, Matlab_bior4.4 NFSAs0= ; PlotFlg=NFSAs0(1); j2=j1+Gr2(1); ordr0=length(s0(j2).P); if RealOnly0 then // use the real part only RawX0=real(s0(j2).S); ordr0=ordr0*2; else RawX0= ; end std0=sqrt(mean(abs(RawX0-sum(RawX0)/length(RawX0)).^2)); // std., RMSE RawX0=RawX0/std0; =TwltDenoise(ordr0,RawX0,WltNum(1),NFSAs0,TestNum0,PlotFlg,j1); TestRMSE, end NFSAs0, Tic_Toc_In_Second=toc(), 新浪赛特居士SciteJushi-2014-02-20。 图片 1.Scilab与Matlab的方差计算的差别 图片 2. 处理复信号1至3的结果 图片 3. 处理复信号4至6的结果
原载 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。
原载 http://blog.sina.com.cn/s/blog_729a92140101ly6e.html 前面已经提到“sure”出身不凡，在Matlab小波工具箱中有重要的地位，并以处理图例证实居士的机器系统工作正常。不妨，再试段程序： e=zeros(1,5); format short g; for ii=1:1000 Num0=randi(256,1,1); x=zeros(1,Num0); e1=wentropy(x,'shannon'); e(1)=abs(e1)+e(1); e2=wentropy(x,'log energy'); e(2)=abs(e2)+e(2); thr=rand(1,1); e3=wentropy(x,'norm',1+thr); e(3)=abs(e3)+e(3); e4=wentropy(x,'threshold',thr); e(4)=abs(e4)+e(4); e5=wentropy(x,'sure',thr); e(5)=abs(e5+Num0)+e(5); end; e, 可得 e= 。 这表明：随机(随机试验，1000次)长度的“零序列”，除去“sure熵”外，其余所有熵始终为零；然而，其“sure熵”e5始终等于，序列长度的负值。 在Matlab的函数wentropy的帮助文档资料中，在列出几种熵的定义式之前，如本文后图片所示，明确指出了：The entropy E must be an additive cost function such that E (0) = 0。 所以，其计算的“sure熵”e5并不满足熵的条件： E (0) = 0。但是Matlab的小波工具箱为什么收容了它呢？ 依居士看，根据图片中的定义，e5似乎应该为零，而不该是这个wentropy的结果。可以写新的计算程序，使定义的“sure熵”始终为零。那么，Matlab的函数wentropy的计算可能是错误的，但是，前一博文中的曲线图表明这一问题似乎不是新的孤立的，而且最早发现Matlab用“shannon熵”处理居士的信号失败时信号范数实际上也是归一化了的。 系统的很多例子用zpd模式。熵的定义式中，包含了序列长度的话，那么使用不同滤波器或不同分解深度时，小波包方法的性能是否变化，这还不清楚。 新浪赛特居士SciteJushi-2013-09-24。
标签: 降噪
