xiangkjmethod1.m

强混沌干扰中的谐波信号提取程序.算法参考<物理学报>,50(6),2001.

%与XiangKJMethod.m相比，y0计算时包含了计算相点本身！
clc;clear all;close all;
t_final=100;%ts=0.01;%x0=[0;0;1e-10];
x0=[0;0;1e-10];
Q=109;%Q取奇数
threshold=0.97;
D=65;N=400;Amp=0.003;f=0.05;tao=1;%tao不能大，大了svd分解特征值都差不多
n0=500;%从混沌信号的第一个点就开始取值
[t,x]=ode45('lorenzeq1',[0,t_final],x0);plot(t,x),
figure(1);plot(x(:,2)./max(x(:,2)));
n=[1:1:N];
signal=Amp*sin(2*pi*f*n);
xMax=max(x(n0:tao:(N-1)*tao+n0,1));
r0=x(n0:tao:(N-1)*tao+n0,1)./xMax+signal';
%r0=x(n0:tao:(N-1)*tao+n0,2)+signal';
for i=1:N-D+1
    for j=1:D
        r(i,j)=r0(i-1+j);
    end
end
%以下计算y0
%前面行
for i=1:(Q-1)/2
    l=0;
    for s=i:i+Q-1
        l=l+1;
        for j=1:D
            a(l,j)=r(s,j);
        end
    end
    y0(i,:)=(sum(a,1)-r(i,:))/(Q-1);
end 
%中间行
for i=(Q+1)/2:1:N-D+1-(Q-1)/2
    l=0;
    for s=(i-(Q-1)/2):1:(i+(Q-1)/2)
        l=l+1;
        for j=1:1:D
            a(l,j)=r(s,j);
        end
    end
    y0(i,:)=(sum(a,1)-r(i,:))/(Q-1);
end   
%后部行   
for i=N-D+1-(Q-1)/2+1:1:N-D+1
    l=0;
    for s=i:-1:(i-Q+1)
        l=l+1;
        for j=1:1:D
            a(l,j)=r(s,j);
        end
    end
   y0(i,:)=(sum(a,1)-r(i,:))/(Q-1);
end      
ui=r-y0;
for i=1:1:N-D+1
    ui(i,:)=ui(i,:)./norm(ui(i,:));  
end    
%以下计算A(:,:,i),i表示与ui的第i行对应
%前面每行对应的A
for i=1:(Q-1)/2
    l=0;
    %for s=i+1:i+Q-1%不要包括i
    for s=i:i+Q-1    
        l=l+1;
        for j=1:D
            A(l,j,i)=ui(s,j);
        end
    end
end 
%中间每行对应的A
for i=(Q+1)/2:1:N-D+1-(Q-1)/2
    l=0;
    for s=(i-(Q-1)/2):1:(i+(Q-1)/2)
        %if s~=i  %不要包括i
           l=l+1;
            for j=1:1:D
                A(l,j,i)=ui(s,j);
            end
        %end
    end
end   
%后部每行对应的A  
for i=N-D+1-(Q-1)/2+1:1:N-D+1
    l=0;
    %for s=i-1:-1:(i-Q+1)%不要包括i
    for s=i:-1:(i-Q+1)
        l=l+1;
        for j=1:1:D
            A(l,j,i)=ui(s,j);
        end
    end
end   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r_transport=r';
for i=1:1:N-D+1
    [U,S,V]=svd((A(:,:,i))');
    Si(:,:,i)=S;%定义Si为了便于跟踪错对。
    sum0=sum(sum(S));
    summ=S(1,1);k=1;
    while summ/sum0<threshold
        k=k+1;
        summ=summ+S(k,k);
    end    
    TyMl=U(:,1:k);
    m = size(TyMl,1);
    TyMl_CZ=eye(m)-TyMl*TyMl';
    s_transport(:,i)=TyMl_CZ*r_transport(:,i);   
end 
ss=s_transport';%ss为方阵，每行代表一个向量
s=zeros(N-D+1,N);
%进行错位处理，好相加
for i=1:1:N-D+1
    for j=1:D
    s(i,j+(i-1))=ss(i,j);
    end
end 
NoZero=~~s;
NoZeroSum=sum(NoZero,1);
sSum=sum(s,1);
s_pie=sSum./NoZeroSum;
ff=(1:1:N)*f/2/N;
figure(2)
plot(ff,10*log10(abs(fft(x(n0:tao:(N-1)*tao+n0,1)))))
title('总信号的频谱')
figure(3)
plot(s_pie)
title('谐波信号的波形')
hold on;
plot(signal,':r');
figure(4)
plot(ff,10*log10(abs(fft(s_pie))))
title('谐波信号的频谱')