ui_traj.m

有关信道估计和信道均衡的仿真程序

% this handles the "calc trajectory" action for the BERGulator

% Copyright 1997-1998 Phil Schniter

 % keep all variables local in scope
 function [] = ui_traj(f_init_manual);

 % declare global variables
 berg_global; ui_private_global;
 color = color;				% seems to fix matlab5.3 bug

 if nargin==1,				% manual initialization
   if length(f_init_manual)~=N_f, 
     error('specified initialization has incorrect length');
   end;
   f_init = f_init_manual(:);
 elseif ishandle(h_trace_c)&real_source,% initialize via ginput on surface
   set(h_msg,'BackgroundColor',[0 1 0]);
   set(h_msg,'FontAngle','italic', 'String','click on surface...'); drawnow;
   set(fig_menu,'HandleVisibility','off');
   figure(fig_surf); zoom off; f_init = ginput(1).'; zoom on;	% mouseclick
   set(fig_menu,'HandleVisibility','on');
   set(h_msg,'BackgroundColor',[.8 .8 .8]);
   set(h_msg,'FontAngle','normal', 'String',''); drawnow;
 else,					% initialize from "Init Spike"  
   f_init = zeros(N_f,1);
   spike = get(h_spike,'Value')-1;
   if is_FSE,
     f_init(2*spike+[1;2]) = [1,1]/sqrt(2);	% double-spike
   else
     f_init(spike+1) = 1;			% single-spike
   end
 end

 % get the relevant parameters from ui
 N = str2num(get(h_num,'String'));              % length of BS input sequence
   if N < 2, error('must be at least two iterations '); end;
   if N ~= ceil(N/D_u)*D_u,
     N = ceil(N/D_u)*D_u;   			% avoid non-integer ratios!
     set(h_num,'String',num2str(N));
   end;
 mu = str2num(get(h_step,'String'));            % algorithm step size

 % calculate godard radii for DSE-CMA
 if dse_alpha ~= dse_alpha_old,
   set(h_msg,'FontAngle','italic','String','calculating parameters...');drawnow;
   if real_source,                      	% PAM
     dse_gamma = calc_dsegamma_PAM(dse_alpha,alphabet,kappa);
   else,                                	% QAM
     dse_gamma = calc_dsegamma_QAM(dse_alpha,alphabet,kappa);
   end
   dse_alpha_old = dse_alpha;
   set(h_msg,'FontAngle','normal','String',''); drawnow;
 end

 % print ui message
 set(h_msg,'FontAngle','italic', 'String','calculating trajectory...'); 
 drawnow;

 % setup, allowing for arbitrarily long equalizers
 charge_dly = 2*ceil(N_c/2);
 s = zeros(1,(1+is_FSE)*N+2*charge_dly+N_f-2);	% BS/FS source
 N_s = floor(length(s)/(1+is_FSE));
 if real_source,
   s(1:1+is_FSE:N_s*(1+is_FSE)) = make_pam(M_s,N_s,1);
 else
   s(1:1+is_FSE:N_s*(1+is_FSE)) = make_qam(M_s,N_s,1);
 end
 r = conv(c,s);  r = r(charge_dly+1:length(r)-charge_dly+1);  % avoid charge
 clear s 					% save memory
 if real_noise,					% equalizer input
   r = r + sqrt(sigma2_n)*randn(1,length(r));	% ...with real noise
 else
   r = r + sqrt(sigma2_n/2)*(randn(1,length(r))...
	+ j*randn(1,length(r)));		% ...with cmplx noise
 end
 D_f = D_u*N_f*ceil(N/25000);			% tap record decimation factor
 D_e = D_u*ceil(N/25000);			% error data decimation factor

 % implement one of various algorithms
 alg = get(h_alg,'Value');
 if alg == 1,		% straight-up CMA
   [f,y,srcme] = berg_algs(N,N_f,D_f,D_e,is_FSE,kappa,mu,r,f_init,alg,D_u);
 elseif alg == 2,	% normalized CMA (N-CMA)
   [f,y,srcme] = berg_algs(N,N_f,D_f,D_e,is_FSE,kappa,mu,r,f_init,alg,0);
 elseif alg == 3,	% signed-error CMA (SE-CMA)
   [f,y,srcme] = berg_algs(N,N_f,D_f,D_e,is_FSE,se_gamma,mu,r,f_init,alg,0);
 elseif alg == 4,	% signed-regressor CMA (SR-CMA)
   [f,y,srcme] = berg_algs(N,N_f,D_f,D_e,is_FSE,kappa,mu,r,f_init,alg,0);
 elseif alg == 5,	% signed-error, signed-regressor CMA (SS-CMA)
   [f,y,srcme] = berg_algs(N,N_f,D_f,D_e,is_FSE,se_gamma,mu,r,f_init,alg,0);
 elseif alg == 6,	% dithered signed-error CMA (DSE-CMA)
   [f,y,srcme] = berg_algs(N,N_f,D_f,D_e,is_FSE,dse_gamma,mu,r,f_init,...
			   alg,dse_alpha);
 elseif alg == 7,	% weerackody's signed-godard algorithm (WSGA)
   [f,y,srcme] = berg_algs(N,N_f,D_f,D_e,is_FSE,wsga_gamma,mu,r,f_init,alg,0);
 elseif alg == 8,	% pure gradient descent (CME-GD)
   if real_source & real_noise, [f,y,srcme] =... 
	berg_cmegd(N,N_f,D_f,D_e,is_FSE,kappa,mu,r,f_init,0,sigma2_n,C);
   elseif real_source, [f,y,srcme] =...
	berg_cmegd(N,N_f,D_f,D_e,is_FSE,kappa,mu,r,f_init,1,sigma2_n,C);
   else, [f,y,srcme] =...
	berg_cmegd(N,N_f,D_f,D_e,is_FSE,kappa,mu,r,f_init,2,sigma2_n,C);
   end;
 else,
   error('unknown algorithm number in ui_traj.m');
 end
 f = [f_init,f];				% prepend initialization
 f_cur = f(:,length(f(1,:)));
 clear r 					% save memory

 % MS error history (back-calculated from equalizer history)
 h_best = C*f_cur;
 [dum,delta_des] = max(abs(h_best));		% find system delay/sign
 f_des = F(:,delta_des);			% local MMSE eq 
 h_des = zeros(N_h,1);  
 if isreal(c), 					% +/- 1 ambiguity
   h_des(delta_des) = sign(h_best(delta_des));
   f_des = f_des*sign(h_best(delta_des));
 else, 						% phase ambiguity
   h_des(delta_des) = exp(j*angle(h_best(delta_des)));
   f_des = f_des*exp(j*angle(h_best(delta_des)));
 end
 mse_des = abs(E(delta_des,delta_des));		% avoid negative mse's
 mserr = zeros(1,floor((N-1)/D_f)+1);
 for i=1:floor((N-1)/D_f)+1,
   mserr(i) = norm(h_des - C*f(:,i))^2 + sigma2_n*norm(f(:,i))^2;
 end;

 % clear user message and remember that a valid simulation has run
 set(h_msg,'FontAngle','normal', 'String','');
 sim_ran = 1;

 % plot trajectories on cost surface and set color
 if ishandle(h_trace_c)&real_source,
   figure(fig_surf);
   hold on; h_trace_new = plot(f(1,:),f(2,:)); hold off;
   set(h_trace_new,'Color',color);
   h_trace_t = [h_trace_t,h_trace_new];
   set(h_clr,'Enable','on'); drawnow; 		% enable clearing of plots
 end; 

 % exit routine if trajectory plots are not desired
 if no_traj_plots, 
   ui_clr_traj; 				% clear trajectory traces
   return; 					% exit now!!
 end;

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Optional Graphical Output... %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 % check for divergence
 if ~(min(isfinite(y))), 
   error(['ALGORITHM UNSTABLE!  Reduce step-size. (Try 0 < mu < ',...
	num2str(0.1/norm(c)^2/N_f,1),')']); 
 end;
 
 % establish max trace length for all histories
 if isempty(N_max), N_max = N; else, N_max = max(N_max,N); end;

 % smoothed error histories, for plotting
 [filt_num,filt_den] = butter(1,10^(-get(h_smooth,'Value')));
 cme_filt = 10*log10(abs(filtfilt(filt_num,filt_den,...		% CME 
	srcme(1:floor((N-N_h)/D_e)+1).^2)));
 mse_filt = 10*log10(abs(mserr));				% MSE

 % plot error history traces in a separate figure
 figure(fig_hist); 
 set(fig_hist, 'DeleteFcn','ui_del_hist');
 if ishandle(h_trace_h),
   subplot(211);
   axe1 = axis;				% retain axis properties 
   hold on;				% dont clobber traces
   if length(h_trace_mse)>2,
     if ishandle(h_trace_mse(3)), 
       delete(h_trace_mse(3)); 		% needed for matlab idiosynchracy
       h_trace_mse = h_trace_mse(1:2);	% ...where subplot obliterates legend
     end; 
   end; 
   subplot(212);
   axe2 = axis;				% retain axis properties 
   hold on;				% dont clobber traces
 else
   axe1 = [0,0,Inf,-Inf];		% force new axis, allow new traces
   axe2 = [0,0,Inf,-Inf];		% force new axis, allow new traces
 end
 subplot(211); h_trace_new=plot([0:length(cme_filt)-1]*D_e,cme_filt); hold off;
 set(h_trace_new,'Color',color); 
 h_trace_h = [h_trace_h,h_trace_new];
 subplot(212); h_trace_new = plot([0:floor((N-1)/D_f)]*D_f,mse_filt); hold off;
 set(h_trace_new,'Color',color); 
 h_trace_h = [h_trace_h,h_trace_new];

 % always replot MMSE boundary (since length may have changed)
 if ishandle(h_trace_mse), delete(h_trace_mse); end;
 subplot(212);
 hold on;
 h_trace_mse = [plot([0,N_max-1],[1,1]*10*log10(abs(mse_des)),'r:'),...
	plot([0,N_max-1],[1,1]*10*log10(abs(mse_opt)),'r-.')];
 hold off;
 if mse_opt < 10*eps, 
   disp(['Warning: calculated MSE bound of ',...
	num2str(10*log10(abs(mse_des)),3),...
	' dB unreliable due to numerical precision!']);
 end;

 % set history axis to avoid clipping old traces
 subplot(211);
 axis([0,N_max,min([axe1(3),cme_filt-3]),max([axe1(4),cme_filt+3])]);
 title('Smoothed CM-error history');
 xlabel('iteration number'); ylabel('dB');
 zoom on;
 subplot(212);
 axis([0,N_max,min([axe2(3),1.15*mse_filt,1.15*10*log10(abs(mse_opt))]),...
	max([axe2(4),mse_filt+5])]);
 title('Squared-error history'); 
 xlabel('iteration number'); ylabel('dB');
 h_trace_mse = [h_trace_mse, legend(h_trace_mse,'MSE_{mmse-loc}',...
	'MSE_{mmse-glob}')]; 
 zoom on;

 % plot linear tap trajectories in a separate figure
 figure(fig_taps);
 set(fig_taps, 'DeleteFcn','ui_del_taps');
 if ishandle(h_trace_tt),
   subplot(211); axe1 = axis; hold on;
   subplot(212); axe2 = axis; hold on;
 else
   axe1 = [0,0,Inf,-Inf];
   axe2 = [0,0,Inf,-Inf];
 end
 subplot(211); h_trace_new = plot([0:length(f)-1]*D_f,real(f)).'; hold off;
 axis([0,N_max,min([axe1(3),min(real(f))-0.1]),...
	max([axe1(4),max(real(f))+0.1])]);
 set(h_trace_new,'Color',color); h_trace_tt = [h_trace_tt, h_trace_new];
 title('Equalizer coefficient histories');
 xlabel('iteration number'); ylabel('real'); zoom on;
 subplot(212); h_trace_new = plot([0:length(f)-1]*D_f+1,imag(f)).'; hold off;
 axis([0,N_max,min([axe2(3),min(imag(f))-0.1]),...
	max([axe2(4),max(imag(f))+0.1])]);
 set(h_trace_new,'Color',color); h_trace_tt = [h_trace_tt, h_trace_new];
 xlabel('iteration number'); ylabel('imag'); zoom on;

 % excise scatter record
 scat_pct = 20;
 y_scat = y(floor([(1-scat_pct/100)*(floor((N-1)/D_e)+1):...
	floor((N-N_h)/D_e)+1]));

 % plot constellation diagram in a separate figure
 figure(fig_scat);
 set(fig_scat, 'DeleteFcn','ui_del_scat');
 if isreal(y),  % real data
   h_trace_s = plot(y_scat,zeros(1,length(y_scat)),'x');
 else		% cmplx data
   h_trace_s= plot(y_scat,'+');
 end
 y_scat_size = 1.1*max([max(abs(real(y_scat))),max(abs(imag(y_scat)))]);
 axis([-y_scat_size,y_scat_size,-y_scat_size,y_scat_size]);
 axis('square'); 
 ylabel('imaginary'); xlabel('real')
 title(['constellation diagram (last ',num2str(scat_pct),'% of output data)']);
 zoom on;
 clear y_scat	% save memory

 % enable clearing of plots
 set(h_clr,'Enable','on'); drawnow;