easi algorithm for blind source separation

clc
clear
SNR=25;
JSR=5;
M=2;
%% 产生BPSK信号
fd=2400;%  码元速率
fc=9600;%  载波速率 
fs=76800;% 采样速率
rolloff = 0.5;%滚降系数
nsamp = 4;
num=1200;
Sig= randint（1numM）;
source_1=M_MPSK（SigMfdfcfsrolloffnsamp）;
a1=norm（source_1‘fro‘）;
a1=a1/（sqrt（length（source_1）））;
source_1=source_1/a1/2; % 能量归一化
%% 产生单音干扰
f1=8400;f2=9600;f3=10600;
Fs = 1000;
T = 1/Fs;
L = length（source_1）;
t = （0:L-1）*T;
source_2=2*sin（2*pi*f1*t）+0.5*sin（2*pi*f2*t）+sin（2*pi*f3*t）;
a2=norm（source_2‘fro‘）;
a2=a2/（sqrt（length（source_2）））;
source_2=source_2/a2/2; % 能量归一化
%% 混合源信号并加入高斯白噪声
s=[source_1;source_2];
A1=[1 -0.6;0.6 1];
A2=[0.6 -0.8;0.8 0.6];
x=zeros（2L）;
x（:1:19200）=A1*s（:1:19200）;
x（:19201:38400）=A2*s（:19201:38400）;
x（1:）=awgn（x（1:）SNR‘measured‘）;%加入高斯白噪声
x（2:）=awgn（x（2:）SNR‘measured‘）;
%% 分离混合信号并计算PI
y=zeros（2L）;
z=zeros（2L）;
w=eye（2）;
w_old=eye（2）;
Error=zeros（1L）;
p=0.95;
PI_AI7=zeros（1L）;
% k=zeros（2L）;
% m=zeros（2L）;
% sigma=zeros（2L）;
% sum_sigma=zeros（1L）;
% yita=zeros（1L）;
yita_0=0.002;
% aerfa=1e-5;
for i=2:L
    z（:i）=x（:i）;
    y（:i）=w*z（:i）;  
%% 基于翘度的步长自适应
%     k（:i）=p*k（:i-1）+（1-p）*（y（:i）.^4-3*y（:i）.^2）;
%     m（:i）=p*m（:i-1）+（1-p）*k（:i）;
%     sigma（:i）=p*sigma（:i-1）+（1-p）*（ k（:i）- m（:i））.^2;
%     sum_sigma（i）=log10（sum（sigma（:i））/2）;
%     yita（i）=yita_0*exp（-aerfa*sum_sigma（i）^2）;
    
    fy=tanh（3*y（:i））;
    gy=sign（y（:i））;
    Error（i）=sum（abs（y（:i）-w_old*z（:i）））/L;
    w_old=w;
    if Error（i）<2e-3
        F=eye（2）-y（:i）*y（:i）‘-y（:i）*fy‘+fy*y（:i）‘;
        w=w+yita_0*F*w;
        w=w/max（max（w））;%max（max（w））为先找w每一列中最大的数，再找出这一行中最大的数；亦即找出w中最大的数。
    else
        F=eye（2）-y（:i）*y（:i）‘+gy*fy‘-fy*gy‘;
        w=w+yita_0*F*w;
        w=w/max（max（w））;
    end
  if i<19201
      A=A1;
  else A=A2;
  end
c=w*A;
  num_r=size（A1）; 
  for p=1:num_r;                             
        max1（p）=abs（c（p1））;                           
        for q=1:num_r                                
            if max1（p）<=abs（c（pq））
                max1（p）=abs（c（pq））;
            else max1（p）=max1（p）;
            end
        end
  end
    ss2=0;   
    for p=1:num_r
        ss1=0;
        for q=1:num_r
            ss1=ss1+abs（c（pq））/max1（p）;
        end
        ss2=ss2+abs（ss1-1）;
    end
    for q=1:num_r
        max2（q）=abs（c（1q））;
        for p=1:num_r
            if max2（q）<=abs（c（pq））
                max2（q）=abs（c（pq））;
            else max2（q）=max2（q）;
            end
        end
    end
    ss4=0;
    for q=1:num_r
        ss3=0;
        for p=1:num_r
            ss3=ss3+abs（c（pq））/max2（q）;
        end
        ss4=ss4+abs（ss3-1）;
    end
    PI_AI7（i）=ss2/num_r+ss4/num_r;   
end
figure（1）
plot（PI_AI7‘b-‘）;
hold on;
%% 固定步长easi
%% 混合源信号并加入高斯白噪声
s=[source_1;source_2];
A1=[1 -0.6;0.6 1];
A2=[0.6 -0.8;0.8 0.6];
x=zeros（2L）;
x（:1:19200）=A1*s（:1:19200）;
x（:19201:38400）=A2*s（:19201:38400）;
x（1