%自动小波阈值选择去噪（scale-dependent）,lambda=m/0.6745*sqrt(2*ln(n))
%参考文献：
% Ma X, Zhou C, Kemp I J. Automated wavelet selection and thresholding for PD detection[J]. Electrical Insulation Magazine IEEE, 2002, 18(2):37-45.
%主要函数：wavedec,detcoef,wrcoef
function y_denoised=DWTdenoising_automaticthreshold(noisydata,wavebase,nlevel,thresholdtype,threshold)
% y_denoised--去噪后信号
% noisydata--含噪信号-待去噪信号
% wavebase--母小波、小波基函数
% nlevel--分解层数
% thresholdtype--去噪方法（hard-硬阈值，soft-软阈值）
% threshold--去噪时使用的阈值策略Scale: median(C)/0.6745*sqrt(2*ln(nj)); Robust: median(C)/0.6745*sqrt(2*ln(N)); Sqtwolog: sqrt(2*ln(N))
% 读程序之前先了解[C,L]=wavedec得到的结果
N=length(noisydata); %数据长度
if size(noisydata,1)>size(noisydata,2) %处理行向量，列向量转置
noisydata=noisydata';
end
if nargin==1 %查看输入的参数个数(nargin)，比如只输入了noisydata，后面的wavebase等参数均为输入wavebase='db1'; %小波基
nlevel=1; %分解层数
thresholdtype='hard'; %阈值方法
threshold='Scale';
else
if nargin==2
nlevel=1; %分解层数
thresholdtype='hard'; %阈值方法
threshold='Scale';
else
if nargin==3
thresholdtype='hard';
threshold='Scale';
else
if nargin==4
threshold='Scale';
end
end
end
end
[C,L]=wavedec(noisydata,nlevel,wavebase);
L1=fliplr(L(1:nlevel+1));
switch thresholdtype
case 'hard' %硬阈值去噪
for j=1:nlevel
CD=detcoef(C,L,j);
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(L1(j))); %scale-dependent threshold end
end
for k=1:L1(j)
if abs(CD(k))<=lambda
CD(k)=0;
end
end
C(sum(L1(j+1:end))+1:sum(L1(j+1:end))+L1(j))=CD;
end
case 'soft' %软阈值去噪
for j=1:nlevel
CD=detcoef(C,L,j);
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(L1(j))); %scale-dependent threshold end
end
for k=1:L1(j)
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(sum(L1(j+1:end))+1:sum(L1(j+1:end))+L1(j))=CD;
end
otherwise
printf('Wrong Input Parameters!\n');
end
y_denoised=zeros(size(noisydata));
for j=1:nlevel
y_denoised=y_denoised+wrcoef('d',C,L,wavebase,j); %逐层恢复细节系数end
y_denoised=y_denoised+wrcoef('a',C,L,wavebase,nlevel); %恢复近似系数。