自适应奇异值分解(ASVD)局放信号去噪

%自适应奇异值分解（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;