При SNR=15 и M=16 вроде бы работает.
При SNR=35 и М=64 (и при 256) уходит в разнос.
В чем причина может быть?
Код
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);