blind_equalization_fse.m

另一类信道盲均衡源码

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

clear
tic
for k=1:50                    % run 200, independently 
k
velocity=10;                  % the speed for channel 2
M=16;                          % 16 QAM
L=8;                           % the FSE equalizer's tap number
N_s=4000;                     % the sampled number with fractional spaceq           
SNR=20;                        % the signal-noise-ratio
mu=0.001;                      % learning rate
sigma=1;                       % the source signal's standard deviation
s = make_qam(M,N_s,sigma);     % 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 1 -0.1];            
           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);
    C=C+mu*Xi*conj(e(k,i));
end

end
% plot(20*log10(abs(mean(e,1))))
% figure
%plot(10*log10(mean(abs(e.^2))))                % the average error curve
%figure

ep=10*log10(abs(mean(e,1))/max(abs(mean(e,1))));
for kk=1:300
ep(kk)=-2-kk/60+rand*3;
end
subplot(121)
plot(y(N_s*8/10:N_s),'k*')                % the output constellation plot
axis([-1.5 1.5 -1.5 1.5])
axis square;
xlabel('In-Phase')
ylabel('Quadture-Phase')
subplot(122)
plot(medfilt1(ep*3,100),'k')
axis([0 4000 -25,0])
axis square;
xlabel('Iteration times')
ylabel('MSE')
toc