cma.m

恒模算法的简介

% [e-mail]

%%%%%%%%%%%%%%%%%%%% 
% control settings
%%%%%%%%%%%%%%%%%%%%

Le      = 2;       %equlizer length
N       = 4000;  % data length
var_s1  = 1;        % variance of signal
var_s2  = 1;

var_n   = 0;        % variance of noise

%%%%%%%%%%%%%%%%%%%% 
% channel
%%%%%%%%%%%%%%%%%%%%

h11 = [ 0       1];
h21 = [ 1       0];
h31 = [ 1       0];
h41 = [ 1       0];

h12 = [ 0       2];
h22 = [ 1       0];
h32 = [ 1       0];
h42 = [ 1       0];

%%%%%%%%%%%%%%%%%%%% 
% generate CM signals
%%%%%%%%%%%%%%%%%%%%

cs1 = [1 -1 j -j 1+j 1-j -1+j -1-j];
cs2 = [0.5 -0.5 j*0.5 -j*0.5 1.5 -1.5 1.5*j -1.5*j 1+j 1-j -1+j -1-j];

cs1 = cs1/sqrt(mean(abs(cs1).^2));
cs2 = cs2/sqrt(mean(abs(cs2).^2));

s1 = cs1([floor(rand(1, N)*length(cs1))+1]);
s2 = cs2([floor(rand(1, N)*length(cs2))+1]);

%%%%%%%%%%%%%%%%%%%% 
% noise
%%%%%%%%%%%%%%%%%%%%

n = randn(4, N);
n(1,:) = n(1,:)*sqrt(var_n/mean(n(1,:).^2));
n(2,:) = n(2,:)*sqrt(var_n/mean(n(2,:).^2));
n(3,:) = n(3,:)*sqrt(var_n/mean(n(3,:).^2));
n(4,:) = n(4,:)*sqrt(var_n/mean(n(4,:).^2));

%%%%%%%%%%%%%%%%%%%%
% channel output
%%%%%%%%%%%%%%%%%%%%

hs11 = chout(s1, h11)';
hs21 = chout(s1, h21)';
hs31 = chout(s1, h31)';
hs41 = chout(s1, h41)';

hs12 = chout(s2, h12)';
hs22 = chout(s2, h22)';
hs32 = chout(s2, h32)';
hs42 = chout(s2, h42)';

x1 = hs11+hs12+n(1,:);
x2 = hs21+hs22+n(2,:);
x3 = hs31+hs32+n(3,:);
x4 = hs41+hs42+n(4,:);

x = [x1; x2; x3; x4];

%%%%%%%%%%%%%%%%%%%%
% CM
%%%%%%%%%%%%%%%%%%%%

w_init = (rand(4, Le)-0.5)*5;

w=fminunc('J_CM', w_init,[], x, 1);

J_CM(w, x, 1)

%%%%%%%%%%%%%%%%%%%%
% recover signal
%%%%%%%%%%%%%%%%%%%%

y = zeros(1, size(x, 2));
for c=1:size(w, 1)
    y=y+chout(x(c, :), w(c, :))';
end

%%%%%%%%%%%%%%%%%%%%
% residual error for s1 and s2
%%%%%%%%%%%%%%%%%%%%

var_err_cvp1 = [];
var_err_cvp2 = [];
var_err_cvm1 = [];
var_err_cvm2 = [];

for (dly = 0:Le+4)

    s1_dly  = sig_dly(s1, dly);
    s2_dly  = sig_dly(s2, dly);
    
    err1    = y - s1_dly;
    err2    = y - s2_dly;
    
    var_err1 = mean(abs(err1).^2);
    var_err2 = mean(abs(err2).^2);

    var_err_cvp1 = [var_err_cvp1, var_err1];
    var_err_cvp2 = [var_err_cvp2, var_err2];

    err1     = y + s1_dly;
    err2     = y + s2_dly;
    
    var_err1 = mean(abs(err1).^2);
    var_err2 = mean(abs(err2).^2);
    
    var_err_cvm1 = [var_err_cvm1, var_err1];
    var_err_cvm2 = [var_err_cvm2, var_err2];
end

clf;
hold;
plot(0:Le+4, var_err_cvp1);
plot(0:Le+4, var_err_cvm1, 'r');

plot(0:Le+4, var_err_cvp2, 'g');
plot(0:Le+4, var_err_cvm2, 'k');

grid;

[min_valp1 dly_valp1] = min(var_err_cvp1);
[min_valm1 dly_valm1] = min(var_err_cvm1);

if (min_valp1 < min_valm1)
    dly_val1 = dly_valp1-1;
    min_val1 = min_valp1;
    sign1    = 1;
else
    dly_val1 = dly_valm1-1;
    min_val1 = min_valm1;
    sign1    = -1;
end

[min_valp2 dly_valp2] = min(var_err_cvp2);
[min_valm2 dly_valm2] = min(var_err_cvm2);

if (min_valp2 < min_valm2)
    dly_val2 = dly_valp2-1;
    min_val2 = min_valp2;
    sign2    = 1;
else
    dly_val2 = dly_valm2-1;
    min_val2 = min_valm2;
    sign2    = -1;
end

clf; hold;
if min_val1 < min_val2
    fprintf(1, '\n***** s1 recovered with delay %d, signal error value %f *****\n', dly_val1, min_val1);
    
    plot(sign1*real(sig_dly(s1, dly_val1)), 'r')
else
    fprintf(1, '\n***** s2 recovered with delay %d, signal error value %f *****\n', dly_val2, min_val2);
    
    plot(sign2*real(sig_dly(s2, dly_val2)), 'r')
end

plot(real(y));