cma的仿真

clc
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 实现了信号的盲均衡算法
% CMA、MCMA算法
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%% 仿真条件 %%%%%%%%
i = sqrt(-1);
L = 6000;
SNR = 30;

h = [0.26-i*0.1, 0.93-i*0.2, 0.26];
h = h/norm(h);
Nh = length(h);

ISI_cma = 0;
ISI_mma = 0;

Nf = 17;
times = 20;

for loop = 1:times
    M = 4;
    iniphase = 0;
    sa = randi([0 M-1],L,1);
%     sn = (pskmod(sa,M,iniphase))';
%     sn = (qammod(sa,M,iniphase))';
    sn=pskmod(sa,M);
    rn = filter(h,1,sn);
    xn = awgn(rn,SNR,'measured');
    Ri = mean(imag(sn).^4)/mean(imag(sn).^2);
    Rr = mean(real(sn).^4)/mean(real(sn).^2);
    R = mean(abs(sn).^4)/mean(abs(sn).^2);
    
    %%%%%% CMA 算法 %%%%%%
    W_cma = zeros(Nf,1);
    W_cma((Nf 1)/2) = 1;
    muCMA = 0.00005;
    
    for k = 1:L-Nf
        X_cma(:,k) = fliplr(xn(k:k Nf-1));
        Y_cma(k) = W_cma'*X_cma(:,k);
        E_cma(k) = Y_cma(k)*(abs(Y_cma(k))^2-R);
        W_cma = W_cma - muCMA*conj(E_cma(k))*X_cma(:,k);
        
        tmpx = conv(W_cma',h);
        ISIx(k) = sum(abs(tmpx).^2)/max(abs(tmpx))^2-1;
        ISIlgx(k) = 10*log10(ISIx(k));
    end
    ISI_cma = ISI_cma ISIlgx;
    
    %%%%%% MCMA算法 %%%%%%%%%%
    W_mma = zeros(Nf,1);
    W_mma((Nf 1)/2) = 1;
    % muMMA = 0.005;%%8PSK时的步长选取
    muMMA  = 0.00005;%%16QAM时的步长选取
    
    for k = 1:L-Nf
        X_mma(:,k) = fliplr(xn(k:k Nf-1));
        Y_mma(k) = W_mma'*X_mma(:,k);
        Yr = real(Y_mma(k));
        Yi = imag(Y_mma(k));
        
        Er_mma = Yr*(Yr^2-Rr);
        Ei_mma = Yi*(Yi^2-Ri);
        E_mma(k) = Er_mma i*Ei_mma;
        W_mma = W_mma - muMMA*conj(E_mma(k))*X_mma(:,k);
        
        tmpy = conv(W_mma',h);
        ISIy(k) = sum(abs(tmpy).^2)/max(abs(tmpy))^2-1;
        ISIlgy(k) = 10*log10(ISIy(k));
    end
    ISI_mma = ISI_mma ISIlgy;
end
ISI_cma = ISI_cma/times;
ISI_mma = ISI_mma/times;

figure
k = 1:length(ISI_cma);
plot(k,ISI_cma,k,ISI_mma,'r');
xlabel('迭代次数');ylabel('码间干扰ISI(dB)');
figure
subplot(2,2,1)
plot(xn,'.');title('均衡前');
subplot(2,2,2)
plot(Y_cma(500:end),'.');title('CMA算法');
subplot(2,2,3)
plot(Y_mma(500:end),'.');title('RLS-CMA算法');