gradual_change.m

blind signal detection using the gradual change algorithm

%信道发生渐变的收敛性能比较图
clear all;
close all;
clc
num12=1;
MSE1(num12)=1;
num32=1;
MSE3(num32)=2;
while(MSE1(num12)>0.002|MSE3(num32)>MSE1(num12))

p=[1 0 0 1 0 1;1 1 0 1 1 1];
p=p(:,2:6);
init=[1 1 1 1 1];
state=kron(init,[1;1]);
c=[];
for i=1:31
    output=mod(sum(p.*state,2),2);
    state=[output state(:,1:4)];
    c=[c output];
end;
a=2*c(1,:).'-1;
b=2*c(2,:).'-1;
s=[];
for i=1:31
    b=[b(2:31);b(1)];
    s=[s a.*b];
end;
s1=s(:,3)./sqrt(31);
b1=rand(1,1000);


c1=s1;
for i=1:1000
    if(b1(i)>0.5)
        b1(i)=1;
    else
        b1(i)=-1;
    end
end
r1=c1*b1;


h=rand(1,10);
% h2=rand(1,10);
% r1=filter(h,1,r1); 
for i=1:10
    h2(i)=h(i)*(1+(rand()*2-1)/10);
end
r1(:,1:2)=filter(h,1,r1(:,1:2));
r1(:,2:1000)=filter(h2,1,r1(:,2:1000));
SNR=20;
r=awgn(r1,SNR);
lim=1*10^-6;


%---------------------LCDPFCMA-------------------------

w=s1;
u=10*10^-3;
B=eye(31)-s1*s1';
x=zeros(31,1);
H=eye(31);
numerr11=200;
MSE1=0;
MSEF1=0;
num11=1;
num12=201;

%---------------信道变化前
tic
for i=1:200
    X=r(:,1);
    y(i)=w'*X;
    e=2*(y(i)^2-1)*y(i);
    g=e*B*X;
    if norm(g)<0.00000000001
        numerr11=i
        break;
    end
    if i==1                                             
        p=-H*g;
    else
        ppast=p;
        j=x-x1;      
        k=g-g1;
        H=H-(H*j*j'*H)/(j'*H*j)+k*k'/(k'*k);
%         H=H-(H*k*k'*H)/(k'*H*k)+(j*j')/(k'*j);
%         H=H-(H*j*j'*H)/(j'*H*j)+k*k'/(k'*k)+(j'*H*j)*((k)/(k'*k)-(H*j)/(j'*H*j))'*((k)/(k'*k)-(H*j)/(j'*H*j));
%         H=H-(H*k*k'*H)/(k'*H*k)+(j*j')/(k'*j)+(k'*H*k)*((j)/(k'*j)-(H*k)/(k'*H*k))'*((j)/(k'*j)-(H*k)/(k'*H*k));
        p=-H*g;
    end
    g1=g;
    x1=x;
%     u=-((s1'*R*B*p+p'*B*R*B*x)/(p'*B*R*B*p));
    x=x+u*p;
    w=s1+B*x;
end

%---------------信道变化后

for i=200:1000
    X=r(:,2);
    y(i)=w'*X;
    e=2*(y(i)^2-1)*y(i);
    g=e*B*X;
    if norm(g)<0.00000000001
        numerr12=i
        break;
    end
    if i==1                                             
        p=-H*g;
    else
        ppast=p;
        j=x-x1;      
        k=g-g1;
        H=H-(H*j*j'*H)/(j'*H*j)+k*k'/(k'*k);
%         H=H-(H*k*k'*H)/(k'*H*k)+(j*j')/(k'*j);
%         H=H-(H*j*j'*H)/(j'*H*j)+k*k'/(k'*k)+(j'*H*j)*((k)/(k'*k)-(H*j)/(j'*H*j))'*((k)/(k'*k)-(H*j)/(j'*H*j));
%         H=H-(H*k*k'*H)/(k'*H*k)+(j*j')/(k'*j)+(k'*H*k)*((j)/(k'*j)-(H*k)/(k'*H*k))'*((j)/(k'*j)-(H*k)/(k'*H*k));
        p=-H*g;
    end
    g1=g;
    x1=x;
%     u=-((s1'*R*B*p+p'*B*R*B*x)/(p'*B*R*B*p));
    x=x+u*p;
    w=s1+B*x;
end
toc




for i=1:numerr11
    if i==1
        MSE1(i)=(y(i)-b1(1))^2;
    else
        MSE1(i)=MSE1(i-1)+(y(i)-b1(1))^2;
        MSEF1(i)=MSE1(i)-MSE1(i-1);
    end
end
for i=numerr11:200
    MSE1(i)=MSE1(i-1);
end
for i=201:numerr12
        MSE1(i)=MSE1(i-1)+(y(i)-b1(2))^2;
        MSEF1(i)=MSE1(i)-MSE1(i-1);
end
MSE1=MSE1./[1:numerr12];
y1=y;
for i=2:200
    if MSEF1(i)<lim
        num11=i;
        break;
    end
end
for i=201:numerr12
    if MSEF1(i)<lim
        num12=i;
        break;
    end
end
%---------------------LCCMA-------------------------

w=s1;
u=5*10^-3;
B=eye(31)-s1*s1';
x=zeros(31,1);
H=eye(31);
numerr21=200;
MSE2=0;
MSEF2=0;
num21=1;
num22=201;

%---------------信道变化前
tic
for i=1:200
    X=r(:,1);
    y(i)=w'*X;
    e=2*(y(i)^2-1)*y(i);
    g=e*B*X;
    if norm(g)<0.00000000001
        numerr21=i
        break;
    end
    if i==1                                             
        p=-H*g;
    else
        ppast=p;
        j=x-x1;      
        k=g-g1;
        H=eye(31);
        p=-H*g;
    end
    g1=g;
    x1=x;
%     u=-((s1'*R*B*p+p'*B*R*B*x)/(p'*B*R*B*p));
    x=x+u*p;
    w=s1+B*x;
end
%---------------信道变化后
for i=201:1000
    X=r(:,2);
    y(i)=w'*X;
    e=2*(y(i)^2-1)*y(i);
    g=e*B*X;
    if norm(g)<0.00000000001
        numerr22=i
        break;
    end
    if i==1                                             
        p=-H*g;
    else
        ppast=p;
        j=x-x1;      
        k=g-g1;
        H=eye(31);
        p=-H*g;
    end
    g1=g;
    x1=x;
%     u=-((s1'*R*B*p+p'*B*R*B*x)/(p'*B*R*B*p));
    x=x+u*p;
    w=s1+B*x;
end
toc




for i=1:numerr21
    if i==1
        MSE2(i)=(y(i)-b1(1))^2;
    else
        MSE2(i)=MSE2(i-1)+(y(i)-b1(1))^2;
        MSEF2(i)=MSE2(i)-MSE2(i-1);
    end
end
for i=numerr21:200
    MSE2(i)=MSE2(i-1);
end
for i=201:numerr22
        MSE2(i)=MSE2(i-1)+(y(i)-b1(2))^2;
        MSEF2(i)=MSE2(i)-MSE2(i-1);
end
MSE2=MSE2./[1:numerr22];
y2=y;
for i=2:200
    if MSEF2(i)<lim
        num21=i;
        break;
    end
end
for i=201:numerr22
    if MSEF2(i)<lim
        num22=i;
        break;
    end
end

%---------------------LSCMA-------------------------

w=s1;
R=eye(31);
MSE3=0;
MSEF3=0;
num31=1;
num32=201;

tic
for i=1:200
    X=r(:,1);
    y(i)=w'*X;
    d=y(i)./abs(y(i));
    e=y(i)-d;
    R=R+X*X';
    w=w-pinv(R)*X*e;
end
for i=201:1000
    X=r(:,2);
    y(i)=w'*X;
    d=y(i)./abs(y(i));
    e=y(i)-d;
    R=R+X*X';
    w=w-pinv(R)*X*e;
end
toc


for i=1:200
    if i==1
        MSE3(i)=(y(i)-b1(1))^2;
    else
        MSE3(i)=MSE3(i-1)+(y(i)-b1(1))^2;
        MSEF3(i)=MSE3(i)-MSE3(i-1);
    end
end
for i=201:1000
    MSE3(i)=MSE3(i-1)+(y(i)-b1(1))^2;
    MSEF3(i)=MSE3(i)-MSE3(i-1);
end



MSE3=MSE3./[1:1000];
for i=2:200
    if MSEF3(i)<lim
        num31=i;
        break;
    end
end
for i=201:1000
    if MSEF3(i)<lim
        num32=i;
        break;
    end
end
y3=y;



end
plot(MSE1,'k-');
hold on;
plot(MSE2,'k--')
plot(MSE3,'k:')
xlabel('迭代次数');
ylabel('MSE');
title('信道发生渐变的收敛性能比较图');
num12,num22,num32
plot(num12,MSE1(num12),'O')
plot(num22,MSE2(num22),'O')
plot(num32,MSE3(num32),'O')
legend('LC-DFPCMA','LC-CMA','LSCMA','收敛点',1)