cma_16qam_blind_equalization.m

此代码产生一个16QAM信号

%说明：CMA算法，性能由ISI表示
clc;
clear all;
close all;
K=11;%抽头个数
delta1=0.000015;
Num_of_realizations=1;%学习次数
SNR=25;%输入信噪比
sgma=1;%输入信号功率
N=5000;%数据长度
d=1;%最小欧式距离
Fd=1000;%波特率
Fs=1000;%采样速率
 h=[0.0410+0.0109j 0.0495+0.0123j 0.0672+0.0170j 0.0919+0.0235j 0.7920+0.1281j 0.3960+0.0871j 0.2715+0.0498j 0.2291+0.0414j 0.1287+0.0154j 0.1032+0.0119j];
ISI_av1=zeros(1,N);
M=16;%进制数调制
for jt=1:Num_of_realizations
    %产生信源
     xs=randint(1,N,M);
   a=qammod(xs,M);   
    a1=real(a);
    a2=imag(a);
    figure(1)
   plot(a1,a2,'*');
    title('(1)原始的星座图');
    if(jt==1) 
    end;
    R2=mean(abs(a).^4)/mean(abs(a).^2);%计算常模数R
    y1=conv(h,a); 
   x=awgn(y1,SNR,sgma);%加噪声
    x1=real(x);
    x2=imag(x);
    figure(2)
    plot(x1,x2,'*');
    title('(2)均衡器输入信号的星座图');
    W1(1:K)=0; 
    W1((K+1)/2)=1;%均衡器抽土系数的初始值
    %CMA
    for n=K:N
        x_k1=x(n:-1:n-K+1);
        z_k1=x_k1*(W1.');
        z1(n)=z_k1;
        s1(n)=4*z_k1*(abs(z_k1)^2-R2);
        isi1(n)=ISI(h,W1);
        W1=W1-delta1*s1(n)*conj(x_k1);
    end
    if(jt==1) 
    end;
    ISI_av1=ISI_av1+isi1;%均方误差的累加。
end
ISI_av1=ISI_av1/Num_of_realizations; %取平
I2=real(z1);
Q2=imag(z1);
figure(3)
plot(I2(N-2000:N),Q2(N-2000:N),'*');
axis('square');
figure(6)
plot(10*log10(ISI_av1));%画ISI，单位dB
xlabel('迭代次数');
ylabel('ISI(dB)');