fig_bse.m

张贤达老师 通信信号处理4－9章的程序 这本书是张老师的经典之作。

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

clear
M=16;                          % 16 QAM
N_s=5000;                      % 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(i)=y(i)*(R2-abs(y(i))^2);
     C=C+mu*Xi*conj(e(i));
 
 end
% plot the command window 
FONTSIZE=0.315;         % 单位为厘米;  
LINEWIDTH_AXES=1.2;%140/42;  % corresponding to 0.140mm 
LINEWIDTH_WAVE=0.5;%180/42;  % corresponding to 0.180mm
figNumber=figure('Name','The Output of BSE equalizer', 'Visible','on');
set(gcf,'Color','white');
plot(y(N_s*8/10:N_s),'r*','LineWidth',LINEWIDTH_WAVE);  
set(gca,'XColor','blue','YColor','blue','FontUnits','centimeters','FontName','Times New Roman','FontSize',FONTSIZE,'LineWidth',LINEWIDTH_AXES);
set(figNumber,'PaperPosition',[0.634518 6.34518 20.3046*0.8 15.2284*0.8]);     %Shrink the figure size