lms.txt

自适应LMS盲均衡算法

N = 5000;
M = N/10;
p = 0.5;
x =2*binornd(1, p, 1, N)-1;
ss = zeros(1,N);
sss = zeros(1,N);
p = 2;
for i = 1:N;
    ss(i) = x(i)^(2*p);
    sss(i) = x(i)^(4*p);
end

m = mean(ss);
l = mean(sss);
z = l/m;
Rz = z;

n = 1:3;
W = 3.5;
h = (1/2) * (1 + cos((2 * (pi/W))*(n-2)));
ww = 0:2 * pi/128:2 * pi;
ww = ww(1:128);
figure(1), plot(ww, abs(fft(h,128))); 
title('channel frequency response');

xlabel('frequncy between 0-2*pi'); ylabel('amplitude');

v = 0.0618 * rand(1,N+2);
u = conv(x,h);
r = u + v;

mu = 0.068;
w = rand(1,N);
y = zeros(1,N);
s = zeros(1,N);
k = 11;

for n = (k+1):M;
    
       
    e(n) = x(n-6) - w((n-1):(n+9)) * r(n:-1:(n-10))';
    
    w((n):(n+10)) = w((n-1):(n+9)) + mu*e(n)*r((n):-1:(n-10));
    
    y(n-6) = w((n):(n+10)) * r((n):-1:(n-10))';
    
    Rp = (y(n-6)/abs(y(n-6))) * (abs(y(n)) + Rz*(abs(y(n)))^(p-1) + (abs(y(n)))^(2*p-1));
    
    if y(n-6) > 0
        s(n-6) = Rp;
    elseif y(n-6) < 0
        s(n-6) = -Rp;
    end
  
       
    error1(n) = (x(n)-y(n))^2;
    error2(n) = (x(n)-s(n))^2;
end

figure(2);
subplot(211);
plot(x(1:M));
title('sent sequence Ak');
subplot(212);
plot(r(1:M));
title('the sequence passed through the channel and noise is added, Rk');

figure(3);
subplot(211);
plot(y(1:M));
title('estimated sequence before slicer, Qk');
subplot(212);
plot(s(1:M));
title('estimated sequence after slicer, A_hat_k');

figure(4);
subplot(211);
plot(error1(1:M));
title('error signal before slicer');
subplot(212);
plot(error2(1:M));
title('error signal after slicer');

www = 0:2 * pi / 512:2 * pi;
www = www(1:512);

figure(5);
subplot(211);
plot(www,abs(fft(w((M):(M+k-1)),512)));
title('frequence response of equalizer at 1000th sample');
subplot(212);
plot(w((M):(M+k-1)));
title('impulse response of equalizer at 1000th sample');


mu = 0.0068;
k = 11;

for n = M:N;
          
    e(n) = s(n-5) - w((n-1):(n+9)) * r(n:-1:(n-10))';
    
    w((n):(n+10)) = w((n-1):(n+9)) + mu*e(n)*r((n):-1:(n-10));
    
    y(n-5) = w((n):(n+10)) * r((n):-1:(n-10))';
    
    Rp = (y(n-6)/abs(y(n-6))) * (abs(y(n)) + Rz*(abs(y(n)))^(p-1) + (abs(y(n)))^(2*p-1));
    
    if y(n-5) > 0
        s(n-5) = Rp;
    elseif y(n-5) < 0
        s(n-5) = -Rp;
    end
  
       
    error1(n) = (x(n)-y(n))^2;
    error2(n) = (x(n)-s(n))^2;
end

figure(6);
subplot(211);
plot(x(M:N));
title('sent sequence Ak');
subplot(212);
plot(r(M:N));
title('the sequence passed through the channel and noise is added, Rk');

figure(7);
subplot(211);
plot(y(M:N));
title('estimated sequence before slicer, Qk');
subplot(212);
plot(s(M:N));
title('estimated sequence after slicer, A_hat_k');

figure(8);
subplot(211);
plot(error1(M:N));
title('error signal before slicer');
subplot(212);
plot(error2(M:N));
title('error signal after slicer');

www = 0:2 * pi / 512:2 * pi;
www = www(1:512);

figure(9);
subplot(211);
plot(www,abs(fft(w((N):(N+k-1)),512)));
title('frequence response of equalizer at 1000th sample');
subplot(212);
plot(w((N):(N+k-1)));
title('impulse response of equalizer at 1000th sample');