blind_equalization_bse.m

信道盲均衡源码

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

clear
echo off
M=16;                          % 16 QAM
for k=1:50                    % run 200, independently 
    k
N_s=4000;                      % the sampled number with fractional space
L=8;
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
velocity=10;                   % the speed for channel 2
fs=10240;                      % the sampling frequency

%%%%%%%%%%% channel 1 %%%%%%%%%%%%%%
  B=[0.2 0.5 1 -0.1];            % the channel's transfer function
  A=1;
  r=filter(B,A,s);               % the source passed through the channel, and this is the result

%%%%%%%%%%%%%% channel 2( Rayleigh LTV channel) %%%%%%%% 
%             h(1,:)=channel_fft_fun(fs,velocity,N_s);
%             h(2,:)=channel_fft_fun(fs,velocity,N_s);
%             h(3,:)=channel_fft_fun(fs,velocity,N_s);
%             h(4,:)=channel_fft_fun(fs,velocity,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:N_s
%                if (i<4)
%                    r(i)=fliplr(s(1:i))*h(1:i,i);
%                else
%                    r(i)=fliplr(s(i-3:i))*h(:,i);
%                end
%            end    

w=10^(-SNR/20)*1/sqrt(2)*(randn(1,N_s)+j*randn(1,N_s));           % the noise corresponds to the given SNR
x=r+w;                         % the received signal

C=[2;zeros(2*L,1)];          % initialize the equalizer  

%%%%%%%%%%% calculate the recovered signal y(i)            %%%%%%%%%%%%
%%%%%%%%%%% the CMA error e(i)                             %%%%%%%%%%%%

 for i=1:N_s
     if (i<L+1)
         Xi=[zeros(1,L+1-i) x(1:i+L)];
     elseif (i+L>N_s)
         Xi=[x(i-L:N_s) zeros(1,L+i-N_s)];
     else
         Xi=[x(i-L:i+L)];
     end
     Xi=conj(Xi');
     y(i)=C'*Xi;
     e(k,i)=y(i)*(R2-abs(y(i))^2);
     C=C+mu*Xi*conj(e(k,i));
 
 end
end
 subplot(121)
 plot(y(N_s*8/10:N_s),'k*')                % the output constellation plot
 axis square;
 xlabel('In-Phase')
 ylabel('Quadture-Phase')
 subplot(122)
 ep=10*log10(abs(mean(e,1))/max(abs(mean(e,1))));
 plot(0.8*medfilt1(ep,10),'k')
 axis square;
 xlabel('Iteration times')
 ylabel('MSE')