clear;
SNR = 15;
M = 16;  %
%---Формирование опорного созвездия----------------

k = log2(M);
alphaRe = [-(2*sqrt(M)/2-1):2:-1 1:2:2*sqrt(M)/2-1];
alphaIm = [-(2*sqrt(M)/2-1):2:-1 1:2:2*sqrt(M)/2-1];
%
ip = [0:M-1];
ipBin = dec2bin(ip.');
%
ipDecRe = bin2dec(ipBin(:,[1:k/2]));
ipGrayDecRe = bitxor(ipDecRe,floor(ipDecRe/2));
%
ipDecIm = bin2dec(ipBin(:,[k/2+1:k]));
ipGrayDecIm = bitxor(ipDecIm,floor(ipDecIm/2));
%
modRe = alphaRe(ipGrayDecRe+1);
modIm = alphaIm(ipGrayDecIm+1);
% опорное созвездие
QAMmod = modRe + j*modIm;
%--------------------------------------------------
L = 9; % Длина эквалайзера
u_cma = 1e-4;
u_lms = 2e-3;

W = 3.2;% Свойства канала

N = 10000; %размер сигнала

% Радиус для СМА
qaml=(abs(QAMmod)).^4;
qam2=(abs(QAMmod)).^2;
R2=qaml/qam2;


a = zeros(N,1); % Исходна последовательность
x = zeros(N,1); % Последовательность после канала
y = zeros(N-(L-1),1);
y1 = zeros(N-(L-1),1);


% ------Формирование потока данных-------

Ind_d=randsrc(1,N,1:M);
a=QAMmod(Ind_d);

%-------------------------------------------------------
h = zeros(1,3); % ИХ канала
for n =1:3,
    h(n) = 0.5*(1 + cos(2*pi*(n-2)/W));
end

x=filter(h,1,a); %прохождение через канал
x = awgn(x,SNR);
%-------------------------------------------------------


%------Обработка данных--------------------
for k=1:N-L+1,
    X(k,:) = x(k:k+L-1); % Преобразование входного потока в матричное представление относительно тапов эквалайзера
end

w(1,1:L) = 0; % Подготовка эквалайзера
w(1,(L+1)/2) = 1; % Начальная установка эквалайзера

for k = 1:N-(L-1),
    y(k) = X(k,:) * w(k,:)';
    e(k) = y(k)*(R2 - abs(y(k))^2);
    w(k+1,:) = w(k,:) + u_cma * e(k)' * X(k,:); % Формирование новой оценки эквалайзера
    
% Нахождение ближайшей опорной точки для дальнейшего применения DD-LMS
raz=QAMmod-y(k);
min_raz=find(raz==min(raz));
y1(k)=QAMmod(min_raz(1));

%------------------------------------------------------------
    %DD-LMS
    if (abs(y1(k) - y(k)) < 0.75) % Запуск DD-LMS при определенном пороге
        e(k) = y1(k) - y(k);
        if (abs(e(k)) < 0.5)
            w(k+1,:) = w(k,:) + 2*u_lms * abs(e(k))^2 * e(k)' * X(k,:);
        else
            w(k+1,:) = w(k,:) + 2*u_lms * e(k)' * X(k,:);
        end
    end
end

%--Вспомогательные операции-------------------------------------


figure(1);
subplot(2,2,1),plot(x(fix(N/2):N-(L-1)),'*');
subplot(2,2,2),plot(y(fix(N/2):N-(L-1)),'*');

elog = 10*log(abs(e).*abs(e))/log(10);
subplot(2,2,3:4),plot(elog);

qaml=sum((abs(QAMmod)).^4);
qam2=sum((abs(QAMmod)).^2);
R2=qaml/qam2;