adaptive_fractional_equalizer.m

信道自适应盲均衡算法源码

%%%%%%%%%%%%%% Blind Equalization using Constant Modulus Criterion (FSE)  %%%%%%%%%

clear
echo off
tic
for k=1:200    % run 200, independently 
k
velocity=10;                  % the speed for channel 2
L=5;                           % the FSE equalizer's tap number
N_s=5000;                     % the sampled number with fractional spaceq           
SNR=20;                        % the signal-noise-ratio
mu=0.02;                      % learning rate
sigma=1;                       % the source signal's standard deviation
%s = make_qam(M,N_s,sigma);     % generate source signal
s = 2*randint(1,N_s,[0,1])-1;     % generate source signal
R2=mean(abs(s).^4)/mean(abs(s).^2);              % the Godard constant
fs=10240;                      % the sampling frequency

  for i=1:2*N_s                % fill with zeros in order to make the source signal fractional spaced
      if (rem(i,2)==0)
          s_fs(i)=0;
      else
          s_fs(i)=s((i+1)/2);
      end
  end
 %%%%%%%%%%%%% channel 1(LTI) %%%%%%%%%%%%%%%%%%%%%%%%%%
           B=[0.2 0.5 0.1 -0.1];  
           %B=[0.05 -0.063 0.088 -0.126 -0.25 0.9047 0.25 0 0.126 0.038 0.088];
           A=1;
           r=filter(B,A,s_fs);               % the source passed through the channel, and this is the result

%%%%%%%%%%%%%% channel 2( Rayleigh LTV channel) %%%%%%%% 
%            h(1,:)=channel_fft_fun(fs,velocity,2*N_s);
%            h(2,:)=channel_fft_fun(fs,velocity,2*N_s);
%            h(3,:)=channel_fft_fun(fs,velocity,2*N_s);
%            h(4,:)=channel_fft_fun(fs,velocity,2*N_s);
%            Q=diag([1/sqrt(mean(abs(h(1,:)).^2)) 1/sqrt(mean(abs(h(2,:)).^2)/10^(-3/10)) 1/sqrt(mean(abs(h(3,:)).^2)/10^(-10/10)) 1/sqrt(mean(abs(h(4,:)).^2)/10^(-15/10))]);
%            h=Q*h;       
%            for i=1:2*N_s
%               if (i<4)
%                   r(i)=fliplr(s_fs(1:i))*h(1:i,i);
%               else
%                   r(i)=fliplr(s_fs(i-3:i))*h(:,i);
%               end
%           end    

w=10^(-SNR/20)*1/sqrt(2)*(randn(1,2*N_s)+j*randn(1,2*N_s));           % the noise corresponds to the given SNR
x=r+w;                         % the received signal
C=[1;zeros(L-1,1);0;zeros(L-1,1)];          % initialize the equalizer  

%%%%%%%%%% generate the fractionally spaced sampled signal %%%%%%%%%%%%
%%%%%%%%%% x_even: the even sampled signal                 %%%%%%%%%%%%
%%%%%%%%%% x_odd:  the odd sampled signal                  %%%%%%%%%%%%
 for i=1:2*N_s
     if (rem(i,2)==0)
         x_even(i/2)=x(i);
     else
         x_odd((i+1)/2)=x(i);
     end
 end
%%%%%%%%%%% calculate the recovered signal y(i)            %%%%%%%%%%%%
%%%%%%%%%%% the CMA error e(i)                             %%%%%%%%%%%%
for i=1:N_s
    if (i<L)
        Xi_odd=[zeros(1,L-i) x_odd(1:i)];
        Xi_even=[zeros(1,L-i) x_even(1:i)];
 
    else
        Xi_odd=[x_odd(i-L+1:i)];
        Xi_even=[x_even(i-L+1:i)];
    end
    Xi=conj([Xi_odd Xi_even]');
    y(i)=C'*Xi;
   % e(k,i)=y(i)*(R2-(abs(y(i)))^2);
      e(k,i)=s(i)-y(i); 
   C=C+mu*Xi*conj(e(k,i));
end
%  estimated_c=[0 0 0 0 0 1 0 0 0 0 0];	% initial estimate of ISI
%   for k=1:N-2*K,
%     y_k=y(k:k+2*K);
%     z_k(k)=estimated_c*y_k.';
%     e_k=info(k)-z_k(k);
%     estimated_c=estimated_c+delta*e_k*y_k;
%     mse(k)=e_k^2;
%     echo off ;
%   end;
end
% plot(20*log10(abs(mean(e,1))))
% figure
ep=10*log10(mean(abs(e.^2)))/(10*max(log10(mean(abs(e.^2)))));
 plot((1:length(ep))/5,medfilt1(ep,10),'k');  % the average error curve
 xlabel('Iteration times')
 ylabel('MSE(dB)')
 figure
plot(y(N_s*8/10:N_s),'*')                % the output constellation plot

toc