matlab最优离散⼩波变换（DWT）信号去噪
%⾃适应⼩波去噪（能量最⼤原则-最优⼩波基）--效果还可以
%⼩波阈值选择（scale-dependent）,lambda=m/0.6745*sqrt(2*ln(n)) %参考⽂献：
% 李剑, 杨洋, 程昌奎，等. 变压器局部放电监测逐层最优⼩波去噪算法[J]. ⾼电压技术, 2007, 33(8):56-60. %主要函数：wavedec,wrcoef
function
y_denoised=DWTdenoising_optimumwavebase(noisydata,nlevel,threshol dtype,threshold)
% y_denoised--去噪后信号
% noisydata--含噪信号-待去噪信号
% ⼩波库为db⼩波
% nlevel--分解层数
% thresholdtype--去噪⽅法（hard-硬阈值，soft-软阈值）
% 读程序之前先了解[C,L]=wavedec得到的结果
if nargin==1
nlevel=1; %分解层数
thresholdtype='hard'; %阈值⽅法
threshold='Scale';
else
if nargin==2
thresholdtype='hard'; %阈值⽅法
threshold='Scale';
else
if nargin==3
threshold='Scale';
end
end
end
N=length(noisydata);
if size(noisydata,1)>size(noisydata,2)
noisydata=noisydata';
end
initialdata=noisydata;
bestindex=zeros(1,nlevel); %初始化每层最⼤能量对于的db⼩波detailenergy=0;
bestind=0;
for k=1:30 %母⼩波库为db1-db20
wavebase=strcat('db',num2str(k));
[C,L]=wavedec(initialdata,1,wavebase);
tmp1=sum(C(1+L(1):L(1)+L(2)).^2)+detailenergy;
tmp=sum(C(1:L(1)).^2)/(sum(C(1:L(1)).^2)+tmp1);
if tmp>maxenergy %最⼤能量计算
maxenergy=tmp;
bestind=k;
end
end
[C,L]=wavedec(initialdata,1,strcat('db',num2str(bestind)));
detailenergy=detailenergy+sum(C(1+L(1):L(1)+L(2)).^2);
initialdata=C(1:L(1));
bestindex(j)=bestind; %每层最⼤能量对于的db⼩波
end
% bestindex=ones(1,nlevel)*6;
initialdata=noisydata;
y_denoised=zeros(size(noisydata)); %初始化去噪结果
C=[];L=[];
C1=[];L1=[];
switch thresholdtype
case 'hard' %硬阈值去噪
for j=1:nlevel
[C{j},L{j}]=wavedec(initialdata,1,strcat('db',num2str(bestindex(j))));
CD=C{j}(L{j}(1)+1:end);
if strcmp(threshold,'Robust')
lambda=median(abs(CD))/0.6745*sqrt(2*log(N));
else
if strcmp(threshold,'Sqtwolog')
lambda=sqrt(2*log(N));
else
lambda=median(abs(CD))/0.6745*sqrt(2*log(L{j})); %scale-dependent threshold end
if abs(CD(k))<=lambda
CD(k)=0;
end
end
C{j}(L{j}(1)+1:end)=CD;
initialdata=C{j}(1:L{j}(1));
end
CA1=C{nlevel}(1:L{nlevel}(1));
for j=1:nlevel %逐层恢复⼩波系数-----逐层恢复
CA=wrcoef('d',C{j},L{j},strcat('db',num2str(bestindex(j))),1); %每⼀层都要提取⼀个细节系数
CA1=wrcoef('a',[CA1,C{nlevel+1-j}(L{nlevel+1-
j}(1)+1:end)],L{nlevel+1-j},strcat('db',num2str(bestindex(nlevel+1-j))),1);%先提取最后⼀层近似系数for k=j-1:-1:1 %对于每层的细节系数，都要在利⽤恢复近似系数的⽅式恢复
C1=[CA,C{k}(L{k}(1)+1:end)];
CA=wrcoef('a',C1,L{k},strcat('db',num2str(bestindex(k))),1);
end
y_denoised=y_denoised+CA;
end
y_denoised=y_denoised+CA1;
case 'soft' %软阈值去噪
for j=1:nlevel
[C{j},L{j}]=wavedec(initialdata,1,strcat('db',num2str(bestindex(j))));
CD=C{j}(L{j}(1)+1:end);
if strcmp(threshold,'Robust')
lambda=median(abs(CD))/0.6745*sqrt(2*log(N));
else
if strcmp(threshold,'Sqtwolog')
lambda=sqrt(2*log(N));
printf输出格式matlab
else
lambda=median(abs(CD))/0.6745*sqrt(2*log(L(j))); %scale-dependent threshold
end
if abs(CD(k))<=lambda
CD(k)=0;
else
if CD(k)>lambda
CD(k)=CD(k)-lambda;
else
if CD(k)<-lambda
CD(k)=CD(k)+lambda;
end
end
end
end
C{j}(L{j}(1)+1:end)=CD;
initialdata=C{j}(1:L{j}(1));
end
CA1=C{nlevel}(1:L{nlevel}(1));
for j=1:nlevel %逐层恢复⼩波系数-----逐层恢复
CA=wrcoef('d',C{j},L{j},strcat('db',num2str(bestindex(j))),1); %每⼀层都要提取⼀个细节系数
CA1=wrcoef('a',[CA1,C{nlevel+1-j}(L{nlevel+1-
j}(1)+1:end)],L{nlevel+1-j},strcat('db',num2str(bestindex(nlevel+1-j))),1);%先提取最后⼀层近似系数for k=j-1:-1:1 %对于每层的细节系数，都要在利⽤恢复近似系数的⽅式恢复
C1=[CA,C{k}(L{k}(1)+1:end)];
CA=wrcoef('a',C1,L{k},strcat('db',num2str(bestindex(k))),1);
end
y_denoised=y_denoised+CA;
end
y_denoised=y_denoised+CA1;
otherwise
printf('Wrong Input Parameters!\n'); end
测试：
noisydata=randn(1,1000);
hold on
plot(DWTdenoising_optimumwavebase(noisydata,6,'hard','Scale'),'Line Width',1.5)