clc;
clear all;
delta1=0.0002;%迭代步長(zhǎng)因子
SNR=10;%輸入信噪比
sgma=1;%輸入信號(hào)功率 
N=1000;%數(shù)據(jù)長(zhǎng)度
d=1;%最小歐式距離
Fd=1000;%波特率
Fs=1000;%采樣速率
ceta=0;
a1=100;%調(diào)整randa
sum=0;Lep=2;L=3;%信道長(zhǎng)度
M=Lep+L-1;ISI_av=zeros(1,N-3);mse_av=0;
Num_of_realizations=2;%學(xué)習(xí)次數(shù)
%randa=0.9;%跌代因子
m=16;%進(jìn)制數(shù)調(diào)制 
for jt=1:Num_of_realizations
v=zeros(1,N);
Y=zeros(M,N);
w=rand(4,M);C=100000*eye(M);
gm=zeros(M,N);W=zeros(16,N);
b=zeros(1,N);bata=zeros(1,N);
S=qam_gen(m,N,d,Fd,Fs);  %產(chǎn)生信源

for n=Lep+L-1:N
s=conj(S(n:-1:n-Lep-L+2))'; 
h1=[0.6662-0.8427j,1.6323-0.2503j,-0.6617-0.4102j];%信道 
h2=[0.4607+0.5789j,0.5855-2.6912j,1.3273-0.4184j];
H=[h1,0;h2,0;0,h1;0,h2];
x(:,n)=H*s;%通過(guò)信道卷積后的信號(hào)
end
x1=conj([S(1),0,0,0])';x2=conj([S(2),S(1),0,0])';x3=conj([S(3),S(2),S(1),0])';
x(:,1)=H*x1;x(:,2)=H*x2;x(:,3)=H*x3;%求通過(guò)信道后的前3列
Y=awgn(x,SNR,sgma);%加噪聲后的信號(hào)

ym=[zeros(4,3),Y];%前貼零
for n=Lep+L-1:N
     for k=1:M
     z_k=w(:,k)'*ym(:,n+k-M);%信號(hào)一路估計(jì)值; %*exp(-j*ceta);
     z(k)=z_k;
    %z_d=decide_16QAM(z_k,1);
%ceta=ceta+delta1*imag(z_k*conj((z_d-z_k)));
    v(n-3)=v(n-3)+z(k);
     end
     c=norm(w'*H);%調(diào)整均衡后的輸出
ping(n-3)=v(n-3)/M;%信號(hào)估計(jì)值
mean1(n-3)=v(n-3)/M/c;%均衡后的輸出信號(hào)
 u(n-3)=Y(:,n-3)'*C*Y(:,n-3);
 for k=1:M
    e_k=ping(n-3)-bata(n-1)*z(k);%一路誤差
    b(n-3)=b(n-3)+e_k;
   end
 ping_e(n-3)=abs(b(n-3)/4);%平均誤差
 randa(n-3)=[exp(-a1*abs(ping_e(n-3)))+exp(-1/(n-3))]/2;%遺忘因子函數(shù)               
 gm(:,n-3)=C*Y(:,n-3)/(randa(n-3)+u(n-3));km=[zeros(4,3),gm];
   C=[C-gm(:,n-3)*Y(:,n-3)'*C]/randa(n-3); 
   for k=1:M
    e_k=ping(n-3)-bata(n-1)*z(k);
    r(:,k)=bata(n-1)*w(:,k)+km(:,n-M+k)*conj(e_k);%更新公式
   end 
 R=conj([conj(r(:,1))',conj(r(:,2))',conj(r(:,3))',conj(r(:,4))'])';
 bata(n)=norm(R);W=R/bata(n);
 t=norm(W);W=W/t;
 w(:,1)=W(1:4);w(:,2)=W(5:8);w(:,3)=W(9:12);w(:,4)=W(13:16);%更新權(quán)向量w
   mse(n-3)=mean(abs(b(n-3)/4)^2);%求MSE
   isi(n-3)=ISI_k(H,w);%求ISI
end
mse_av=mse_av+mse;
ISI_av=ISI_av+isi;
end
mse_av=mse_av/Num_of_realizations; 
ISI_av=ISI_av/Num_of_realizations; 

figure(1)
plot(10*log10(ISI_av),'-k')
xlabel('iterations');
ylabel('ISI(dB)');
i51=10*log10(ISI_av);
save x51 i51;