bergulator.m

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

%
%                      **  The BERGulator **
%
% The BERGulator is a matlab5-based interactive simulation tool aimed 
% at studying the convergence behavior of the Constant Modulus Algorithm 
% (CMA) when used for direct blind linear channel equalization. 
%
% The BERGulator allows you to specify the channel impulse response 
% coefficients, the level of channel noise, and the source alphabet
% used for experimentation.  The channel may be complex-valued, and
% complex-valued noise is assumed if either the source or the channel 
% is complex-valued.  Furthermore, the channel may be baud-spaced or
% T/2 fractionally-spaced.
%
% Given such a channel/noise combination, The BERGulator is capable of 
% plotting the frequency/impulse responses of the channel, the MMSE 
% optimal equalizer, and the optimal system (chan+eq).  For length-2 
% real-valued channels, it can also plot contours of the CMA surface
% in equalizer space (with an overlay of the principle axes of the 
% quadratic-cost ellipsoid approximating the neighborhoods of CM minima,
% and with *'s and x's indicating the locations of MMSE and MSE minima).  
% When a valid surface is present, The BERGulator can superimpose the 
% CM-error-function sign-boundaries (for real-valued sources) and 
% indicate the magnitude and direction of the (constant) cost-gradient 
% within a given constant-sign region.  This feature is motivated by 
% the convergence behavior of signed algorithms.
%
% Finally, the BERGulator can plot the tap-trajectories of an equalizer 
% adaptation algorithm and plot the CM- and mean-squared error 
% histories.  When an error surface is present and the tap-trajectories
% will be real-valued, the user can select the equalizer intialization 
% from the surface and the resulting trajectory will be superimposed on 
% it. For higher-dimensional and/or complex-valued equalizers, the user
% selects the location of a single/double spike (BSE/FSE) initialization. 
% The equalizer-length, stepsize, and the number of iterations can be 
% varied.  Numerous variants of the CMA algorithm are available for 
% experimentation.
%
% After simulation run has completed, the frequency/impulse response
% of any equalizer along the adaptation history can be displayed.
% In addition the selected equalizer, the local MMSE equalizer, the
% corresponding system response and the local MMSE system response are
% also displayed.  If a surface trajectory is present, this equalizer 
% selection is obtained from the surface.  If only the coefficient 
% history is present, this equalizer is selected from a selected point 
% in the adaptation history.  A round/diamond marker is then placed on 
% graph to indicate a frequency/impulse response analysis.
% 
% When parameter changes are made that invalidate the currently 
% displayed frequency/impulse responses or error surface, the 
% offending object is automatically erased.  On the contrary, any
% analyses based on an simulation run is retained (for comparison to
% other experiments) until the "Clear Traces" button is pressed.
%
% Actions:
%   Cost Surf : Plots contours of the Godard 2-2 cost surface in eq space
%   Sign Bndr : Overlays CM-error-function sign-boundaries on surface
%   Facet Grad: Displays the 2D gradient of the signed-Godard cost surface
%   FREQ      : Frequency response analysis
%   IMP       : Impulse response analysis
%   Adapt Eq  : Runs a simulation using the selected parameters
%   Clr Trces : Erases all adaptation-based tap/error/freq/imp-rspn traces
%   Help      : Prints this help screen.
%   Quit      : Quit this godawful program!
%
% Channel Options:
%   Chan Type : Offers a choice between preselected channel coefficients,
%               user-selectable ("custom") ones, and ".mat" file input 
%               where the coefs must be stored in a vector named "c"
%   Src Type  : A choice between a number of QAM and PAM source alphabets
%   Spacing   : T/2 Fractionally-spaced (FS) or Baud-spaced (BS) 
%   dB SNR    : Signal-to-AWGN level in decibels
%   Normalize : Normalizes the channel coeficient vector (to unit norm)
%   Equal Axs : Forces the surface plot x/y-axes to be equally scaled
%		(unfortunately this interferes with zooming)
%
% Equalizer Algorithms:
%   CMA       : The Godard 2-2 constant modulus algorithm
%   N-CMA     : Normalized CMA
%   SE-CMA    : Signed-error CMA (i.e. Godard 1-1)
%   SR-CMA    : Signed-regressor CMA
%   SS-CMA    : Signed-error, signed-regressor CMA
%   DSE-CMA   : Dithered signed-error CMA
%   WSGA      : Weerackody's "Signed Godard Algorithm"
%   CMA-GD    : Godard 2-2 gradient descent (i.e. CMA avg-system behavior)
%
% Algorithm Options:
%   Eq Len    : The number of coefficients in the equalizer
%   Init Spike: The T-spaced location of single/double spike initialization
%               for BS/FS equalizers
%   Num Iter  : The number of (baud-spaced) symbols used in the adaptation
%   Trace Clr : The color of the trace used to display tap/error histories 
%               and frequency/impulse response analysis  
%   Stepsize  : The stepsize used by the coefficient update algorithm
%   Err Smooth: The amount of smoothing performed on the CM error history
%               (min = no smoothing)
%
%
% The BERGulator keeps all variables local in scope unless the user
% invokes the program "berg_global.m".  In the latter case, the important
% variables and user-interface handles are made available to the user
% so that he/she/it may access these quantites and/or construct an 
% m-file automating the BERGulator user-interface.
%
%
% Misc Details, for those who care...
%  - the T/2 implementation samples on "odd" outputs (assume {0,1,2,...})
%    when the equalizer length is even and "even" outputs when the equalizer
%    length is odd
%  - all "local-MSE" quantities are calculated based on the delay and phase
%    obtained from the maximum-magnitude entry of the latest simulation's 
%    final system response 
%  - the error-trace smoothing uses a 2nd-order butterworth filter
%  - all signed CMA algorithms (except WSGA) compute the sign of a 
%    complex number as the sign of the real part plus sqrt(-1) times
%    the sign of the imaginary part
%  - the dither used in DSE-CMA is distributed as U(-dse_alpha,dse_alpha]
%    where dse_alpha defaults to 1
%
% The most recent version of The BERGulator can be obtained from 
% "http://backhoe.ee.cornell.edu/BERG/sim_tools/BERGulator.html", which 
% is regularly updated as bug fixes and other changes are incorporated.  
% Send your comments to the author, whose address appears below.
%
%
% ----------------------------------------------------------------------
%
%    BERGulator written by:    Phil "Bert" Schniter 
%			       [e-mail]
%                              Blind Equalization Research Group
%                              (http://backhoe.ee.cornell.edu/BERG/)
%                              Cornell University
%                              June 1997
%
%			       Copyright 1997-1998 Phil Schniter
%
%    last updated 2/10/98
%
%
%
%
% See also BERG_GLOBAL

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

 % declare which variables are global
 berg_global; 		% (workspace needs to execute this too)
 ui_private_global;

 % declare output figure numbers and clear them
 fig_surf = 1;
 fig_freq = 2;
 fig_imp  = 3;
 fig_hist = 4;
 fig_taps = 5;
 fig_scat = 6;
 if ishandle(fig_surf), delete(fig_surf); end; 
 if ishandle(fig_freq), delete(fig_freq); end; 
 if ishandle(fig_imp),  delete(fig_imp);  end; 
 if ishandle(fig_hist), delete(fig_hist); end;
 if ishandle(fig_scat), delete(fig_scat); end;
 if ishandle(fig_taps), delete(fig_taps); end;

 % allocate control window
 if exist('fig_menu'), if ishandle(fig_menu), delete(fig_menu); end; end;
 fig_menu = figure('IntegerHandle','off', 'Name','The BERGulator',...
	'NumberTitle','off', 'MenuBar','none', 'Resize','off');
 pos = get(fig_menu,'Position'); 
 set(fig_menu,'Position',[pos(1)+400,pos(2)+100,400,330]);	
 clf; 

 % channel parameters
 uicontrol('Style','frame', 'Position',[140 175 245 135]);
 uicontrol('Style','text', 'Position',[150 275 75 25],'String',...
	'Channel Type');
 h_chan = uicontrol('Style','popupmenu', 'Position',[150 260 75 25],...
	'String','4-tap|6-tap|Custom|From file', 'Callback','ui_chan');
 uicontrol('Style','text', 'Position',[235 275 50 25], 'String','Spacing');
 h_fse = uicontrol('Style','popupmenu', 'Position',[235 260 50 25],...
	'String','FS|BS', 'Callback','ui_fse');
 uicontrol('Style','text', 'Position',[295 275 75 25], 'String','Source Type');
 h_source = uicontrol('Style','popupmenu', 'Position',[295 260 75 25],...
        'String',['BPSK|4-PAM|8-PAM|16-PAM|32-PAM|QPSK|16-QAM|64-QAM',...
	'|256-QAM|1024QAM'], 'Callback','ui_calc_params');
 uicontrol('Style','text', 'Position',[150 230 160 25], 'String',...
	'Channel Coefficients');
 h_coefs = uicontrol('Style','edit', 'Position',[150 215 159 25],...
	'BackgroundColor',0.9*[1 1 1], 'Callback','ui_chan');
 uicontrol('Style','text', 'Position',[320 230 50 25], 'String','dB SNR');
 h_SNR = uicontrol('Style','edit', 'Position',[320 215 50 25],...
	'String','50', 'BackgroundColor',0.9*[1 1 1],...
	'Callback','ui_calc_params');
 h_norm = uicontrol('Style','checkbox', 'Position', [150 185 120 20],...
        'String','Normalize Channel', 'Callback','ui_chan');
 h_axes = uicontrol('Style','checkbox', 'Position', [280 185 90 20],...
        'String','Equal Axes', 'Callback','ui_axes');

 % algorithm parameters
 uicontrol('Style','frame', 'Position',[140 20 245 135]);
 uicontrol('Style','text', 'Position',[150 120 50 25], 'String','Eq Len');
 h_N_f = uicontrol('Style','edit', 'Position',[150 105 50 25],...
        'String','2', 'BackgroundColor',0.9*[1 1 1], 'Callback','ui_N_f');
 uicontrol('Style','text', 'Position',[210 120 100 25], 'String',...
	'Equalizer Algorithm');
 h_alg = uicontrol('Style','popupmenu', 'Position',[210 105 100 25],...
        'String','CMA|N-CMA|SE-CMA|SR-CMA|SS-CMA|DSE-CMA|WSGA|CME-GD');
 uicontrol('Style','text', 'Position',[320 120 50 25], 'String','Init Spike');
 h_spike = uicontrol('Style','popupmenu', 'Position',[320 105 50 25],...
	'String','-1');
 uicontrol('Style','text', 'Position',[150 75 50 25], 'String','Num Iter');
 h_num = uicontrol('Style','edit', 'Position',[150 60 50 25],...
	'String','5000', 'BackgroundColor',0.9*[1 1 1]);
 uicontrol('Style','text', 'Position',[210 75 100 25], 'String','Trace Color');
 h_col = uicontrol('Style','popupmenu', 'Position',[210 60 100 25],...
        'String',['nova lox|sabotage red|orange peel|lemonade',...
	'|fudge brownie|key lime|john deere grn|sky blue|cop blue',...
	'|80s neon|kingpin|slush gray|johnny cash blk'], 'Callback','ui_color');
 set(h_col,'Value',6); 
 uicontrol('Style','text', 'Position',[320 75 50 25], 'String','Stepsize');
 h_step = uicontrol('Style','edit', 'Position',[320 60 50 25],...
	'String','5e-3', 'BackgroundColor',0.9*[1 1 1]);
 uicontrol('Style','text', 'Position',[149 30 97 16], 'String',...
	'Error Smooth:  min');
 h_smooth = uicontrol('Style','slider', 'Position',[250 30 100 16],...
	'Min',0, 'Max',2.6, 'Value',1.3, 'Callback','ui_smooth');
 uicontrol('Style','text', 'Position',[351 30 25 16], 'String',...
	'max');

 % action buttons 
 h_surf = uicontrol('Style','pushbutton', 'Position',[20 280 100 30],...
	'Callback','ui_surf', 'String','CM Cost Surface');
 h_sign = uicontrol('Style','pushbutton', 'Position',[20 240 100 30],...
	'Callback','ui_sign', 'String','Sign Boundaries');
 h_freq = uicontrol('Style','pushbutton', 'Position',[20 200 45 30],...
        'Callback','ui_freq', 'String','FREQ');
 h_impr = uicontrol('Style','pushbutton', 'Position',[75 200 45 30],...
        'Callback','ui_imp', 'String','IMP');
 h_traj = uicontrol('Style','pushbutton', 'Position',[20 125 100 30],...
	'Callback','ui_traj', 'String','Adapt Equalizer');
 h_clr  = uicontrol('Style','pushbutton', 'Position',[20 85 100 30],...
        'Callback','ui_clr', 'String','Clear Traces');
 h_help = uicontrol('Style','pushbutton', 'Position',[20 45 45 30],...
        'Callback','helpwin bergulator', 'String','Help');
 h_quit = uicontrol('Style','pushbutton', 'Position',[75 45 45 30],...
        'Callback','ui_quit', 'String','Quit');

 % text messages
 h_msg2 = uicontrol('Style','text', 'Position',[10 175 120 15],...
        'BackgroundColor',[.8 .8 .8]);
 h_msg =  uicontrol('Style','text', 'Position',[15 20 115 15],...
        'BackgroundColor',[.8 .8 .8]);

 % initialize 
 sim_ran = 0;			% no simulations run yet
 origin = [0,0];		% origin for surface plots
 alg_surf = 0;			% do not plot algorithm-specific surfaces
 N_v = 15;			% number of contour lines on surface
 D_u = 1;			% CMA update decimation factor
 is_FSE = 2-get(h_fse,'Value');	% remember current setting of FSE/BSE 
 SNR_ok = 80;			% an SNR value avoiding numerical difficulties
 N_f = 2;			% a temporary value... 
 N_max = [];			% max trace length
 h_trace_c = [];		% cost-contour trace record
 h_trace_b = [];		% sign boundaries trace record
 h_trace_f = []; 		% freq-resp trace record
 h_trace_ff = [];		% globally MMSE freq-resp trace record
 h_trace_fff = [];		% locally MMSE freq-resp trace record
 h_trace_i = [];		% impulse-resp trace record
 h_trace_ii = [];		% globally MMSE impulse-resp trace record
 h_trace_iii = [];		% locally MMSE impulse-resp trace record
 h_trace_mse = [];		% MSE bound trace record
 h_trace_h = []; 		% CMA/MSE-history trace record
 h_trace_t = []; 		% surface tap-trajectory trace record
 h_trace_tt = []; 		% tap trajectory trace record
 h_trace_s = [];		% scatter plot (const. diagram) trace record

 % initialize algorithm parameters with defaults
 no_traj_plots = 0;		% return trajectory-related plots by default
 dse_alpha = 1;			% DSE-CMA dither amplitude
 dse_alpha_old = 0;		% ...force calculation of dse_gamma
 ncma_alpha = 0.001;		% N-CMA regularization term
 rand('state',sum(100*clock));	% randomize seed

 % apply initializations
 ui_fse;			% handles fse vs. bse
				% ... which propagates through ui_chan.m
				% ... which propagates through ui_calc_params.m
 ui_N_f;			% handles 2-tap vs. longer equalizers
				% ... and propagates through ui_calc_params.m

 % dim invalid actions, and provide something for the cluemasters out there... 
 set(h_clr,'Enable','off');
 set(h_sign,'Enable','off');
 set(h_msg,'String','click a button...');