%自适应奇异值分解（Adaptive Singular Value Decomposition）去噪
%参考文献：Ashtiani M B, Shahrtash S M. Partial discharge de-noising
employing adaptive singular value decomposition[J]. IEEE Transactions
on Dielectrics & Electrical Insulation, 2014, 21(2):775-782.
function y_denoised=ASVDdenoising(noisydata,k)
if size(noisydata,1)>size(noisydata,2)
noisydata=noisydata';
end
N=length(noisydata);
y_denoised=zeros(N,1);
L=ceil(N/3);
% Y=[];
%
% for j=1:N-L
% Y=[Y;noisydata(j:j+L)]; %HANKEL矩阵 ,大矩阵时这种构
建方式（未事先确定矩阵尺寸）将会严重影响计算速度
% end
Y=HankelMatrixCons(noisydata,L);
[U,S,V]=svd(Y); %奇异值分解

g=diag(S);
g1=sort(g(1:100),'ascend');
g1=g1/g(1)*100;
n=length(g1);
std_g=zeros(n,1);

for j=1:n
if floor(std(g1(1:j)))>0
M=j;
break;
end
end
for j=1:n
std_g(j)=floor(std(g1(1:j)));
end

% subplot(2,1,1)
plot(std_g,'-ro')
if nargin==1
k=input('Input M\n');
% k=M;
end

Y1=U(:,1:(n+1-k))*S(1:(n+1-k),1:(n+1-k))*V(:,1:(n+1-k))';
order=n+1-k;
y_denoised=Subfun1STSVD(U,S,V,order);

% 构建(L*K,K=N-L+1)Hankelj矩阵
function X=HankelMatrixCons(data,L)
N=length(data);
if size(data,2)==N
data=data';
end
X=zeros(L,N-L+1); %创建Hankel矩阵
for j=1:L
X(j,:)=data(j:N-L+j)';
end

% 利用有效阶数提取信号（SVD处理） -对角线平均方式
function y=Subfun1STSVD(S,V,D,M)
Num_row=size(S,1); Num_column=size(D,1);
N=Num_row+Num_column-1;
y=zeros(1,N); %数据
Num=zeros(1,N); %每个数据对应的数据个数，最后用于求均值
Y3=zeros(Num_row,Num_column);
Y1=S(:,1:M)*V(1:M,1:M)*(D(:,1:M))';
for j=1:Num_row
Y3(j,:)=fliplr(Y1(j,:));
end
k=1;
for j=Num_column-1:-1:-(Num_row-1)
a=diag(Y3,j);
y(k)=sum(a);
Num(k)=length(a);
k=k+1;
end
y=y./Num;