求助韩敏《混沌。理论与方法》一书关于小波变换平滑噪声的代码？

cwhe_10

具体仿真代码如下：
%**************************************************************************
%    功能：               改进的小波变换噪声平滑方法
%           公式来源：韩敏-混沌时间序列预测理论与方法一P123-125
%**************************************************************************
clc;clear;close all;
odefun=@(t,x)[-10*(x(1)-x(2)); -x(1)*x(3)+28*x(1)-x(2); x(1)*x(2)-8*x(3)/3];
x0=[1 1 1];
tspan=[0:0.01: 50];
[t,y]=ode45(odefun,tspan,x0);

data=y(:,1)';data=data(1:5000);                                 %原始信号
subplot(1,2,1);
plot(data,'r');grid on;title('原始Lorenz信号数据');

V0=var(data);
V_noise=V0/25;                                                           %25%的白噪声方差，公式：P103
d_noise=randn(1,5000);                                              %生成5000个噪声序列
d_noise=d_noise/std(d_noise);                                    %对噪声序列的方差归一化，降为方差为1
d_noise=d_noise-mean(d_noise);                                %对噪声序列的均值去除，使之为0
d_noise=sqrt(V_noise)*d_noise;                                  %相对原始数据25%噪声水平的白噪声生成完毕！

data=data+d_noise;                                                    %含噪声信号
subplot(1,2,2);
plot(data,'r');grid on;title('25%噪声水平的Lorenz信号数据');
%--------------------------------------------------------------------------
%小波分解，提取小波分解系数
port=data;                                                                %被小波分解的信号
N=4;                                                                          %小波分解的层数
DB='db2';                                                                  %小波类型

a=[];b=[];                                                                 %小波分解系数存储矩阵，为N*5000矩阵
[c,l]=wavedec(port,N,DB);                                         %小波分解过程
for i=1:N
    a(i,  :   )=wrcoef('a',c,l,DB,i);                                       %近似系数（低频部分）提取过程
end
for i=1:N
    d(i,  :   )=wrcoef('d',c,l,DB,i);                                       %细节系数（高频部分）提取过程
end
%--------------------------------------------------------------------------
%小波细节系数阈值去噪处理--公式：P124
N1=length(data);                                                      %含噪声信号长度                                    
V=var(data);                                                             %含噪声信号方差

thr=zeros(1,N);
thr(N)=V*sqrt(2*log(N1))/sqrt(N);                            %最高频细节系数对应的阈值去噪值
for i=2:N-1                                                                %次高频细节系数对应的阈值去噪值
    thr(i)=V*sqrt(2*log(N1))/log(i+1);                       %问题出在这里？？？阈值与附件论文的不一样？？？求指导？
end
thr(1)=V*sqrt(2*log(N1));                                         %最低频细节系数对应的阈值去噪值

for j=1:length(data)
    for k=1:N
        if d(k,j)<thr(k)
           d(k,j)=0;
        end
    end
end
%--------------------------------------------------------------------------
%信号重构
c1=[a(N,  :    ) d(4,    :    ) d(3,   :   ) d(2,   :    ) d(1,   :   )];
s1=waverec(c1,l,DB);
figure
plot(s1,'r');grid on;title('25%含噪声信号的改进小波去噪优化方法曲线图');

上述代码的去噪阈值设置与论文的不一样，求大神指点一二！！！

若觉得信息不详细，可在研学论坛搜索：韩敏《混沌时间序列预测理论与方法》中基于小波改进

cwhe_10

研学论坛的那个帖子附件有，韩敏教授的混沌时间序列预测理论与方法一书的PDF下载链接，以及一篇非常不错的论文，均免费下载