teqdemo.m

多载波调制的仿真程序

   % Dmin cannot be smaller than 1 and larger than Dmax
   if vv<1 | vv > Dmax, 
      vv = v; 
   end
   set(gco,'Userdata',vv,'String',num2str(vv))
   return
   
   % set Dmax   
elseif strcmp(action,'setDmax'),
   hndlList=get(gcf,'Userdata');
   Dmin = get(hndlList(7),'Userdata');
   v = get(gco,'Userdata');
   s = get(gco,'String');
   vv = eval(s,num2str(v(1)));
   % Dmax cannot be smaller than 1 and Dmin
   if vv < 1 | vv < Dmin, 
      vv = v;
   end
   set(gco,'Userdata',vv,'String',num2str(vv))
   return
   
   % set N   
elseif strcmp(action,'setN'),
   hndlList=get(gcf,'Userdata');
   v = get(gco,'Userdata');
   s = get(gco,'String');
   vv = eval(s,num2str(v(1)));
   % N has to be a power of 2
   if (floor(log(vv)/log(2))-(log(vv)/log(2))) ~= 0
      vv = v; 
   end
   set(gco,'Userdata',vv,'String',num2str(vv))
   return
   
   % design TEQ   
elseif strcmp(action,'design'),
   set(gcf,'Pointer','watch');
   % get handles
   axHndl=gca;
   fhndlList=get(gcf,'Userdata');
   freqzHndl = fhndlList(1);
   NbHndl = fhndlList(2);
   NwHndl = fhndlList(3);
   NHndl = fhndlList(4);
   CGHndl = fhndlList(5);
   MarginHndl = fhndlList(6);
   DminHndl = fhndlList(7);
   DmaxHndl = fhndlList(8);
   PwrHndl = fhndlList(9);
   AWGNHndl = fhndlList(10);
   CSAHndl = fhndlList(11);
   viewHndl = fhndlList(12);
   helpHndl = fhndlList(13);
   closeHndl = fhndlList(14);
   methodHndl = fhndlList(15);
   methRateHndl = fhndlList(16);
   methSNRHndl = fhndlList(17);
   methSSNRHndl = fhndlList(18);
   methMSEHndl = fhndlList(19);
   methDelayHndl = fhndlList(20);
   methMaxHndl = fhndlList(21);
   
   colors = get(gca,'colororder'); 
  
   % init variables
   Nb = get(NbHndl,'userdata');
   Nw = get(NwHndl,'userdata');
   N = get(NHndl,'userdata');
   codingGain = get(CGHndl,'userdata');
   margin = get(MarginHndl,'userdata');
   Dmin = get(DminHndl,'userdata');
   Dmax = get(DmaxHndl,'userdata');
   totalInputPower = get(PwrHndl,'userdata');
   AWGNpower = get(AWGNHndl,'userdata');
   loopNum = get(CSAHndl,'userdata');
   method = get(methodHndl,'value');
   barflag = 1;
   %number of iteration per TEQ tap
   numIter = 1000;

   
   %generate signal,noise, and channel data
   [recNoisySig,receivedSignal,noise,channel,inputSignal,gamma,fs] = ...
      siggen(N,AWGNpower,loopNum,totalInputPower,codingGain,margin,barflag);
  
  %save tmp1.mat inputSignal, receivedSignal, noise;

   %estimate the power spectrums
   [inputSpecAll, noiseSpecAll, channelGainAll]=...
      specestim(inputSignal,noise,channel,N,barflag);
   
   %calculate SNRs and usable channels
   [subMFBall,subMFB,usedSubs,noiseSpec,channelGain,inputSpec] = ...
   calcsnrs(inputSpecAll,channelGainAll,noiseSpecAll,gamma);

   % chose method and design TEQ
   [B, W, D, MSE, Dv, I,title_str] = selmeth(method,inputSignal,...
      receivedSignal,noise,channel(1:N),...
   Nb,Nw,Dmin,Dmax,barflag,N,numIter,inputSpec,noiseSpec,channelGain,...
   usedSubs,gamma);

   % get performance results
   [SSNR, SNR, subSNR,geoSNRfinal,geoSNRmfb,bDMTfinal,bDMTmfb,...
         RDMTfinal,RDMTmfb,hw,Fh,Fw,colorNoiseaft,Fhw] = ...
      perform(W,B,channel,D,Nb,N,inputSignal,noise,...
      channelGainAll,inputSpecAll,noiseSpecAll,margin,codingGain,fs,subMFBall,usedSubs);
   
   % put results into table
   str = sprintf('%1.3fe6',RDMTfinal/1e6);
   set(methRateHndl,'string',str);
   str = sprintf('%3.1f',SNR);
   set(methSNRHndl,'string',str);
   str = sprintf('%3.1f',SSNR);
   set(methSSNRHndl,'string',str);
   str = sprintf('%3.0f',D);
   set(methDelayHndl,'string',str);
   if MSE == 0, 
      str = '-';
   else
      str = sprintf('%3.1e',MSE);
   end
   set(methMSEHndl,'string',str);
   str = sprintf('%1.3fe6',RDMTmfb/1e6);
   set(methMaxHndl,'string',str);
   
   
   % save resuls to disk in order to use it in other functions 
   save tmpresults title_str N hw B D Nb W Nw subSNR subMFB channel Fh ...
	   noiseSpec colorNoiseaft Dv Dmin Dmax Fhw Fw usedSubs ... 
       inputSignal receivedSignal noise;
   
   TEQequalizerTaps = W;
   TEQtargetImpulseResp = B;
   TEQoptimalDelay = D;
   TEQoriginalChannel = channel;
   TEQequalizedChannel = hw;
   TEQmatchedFilterBoundPerChannel = subMFB;
   TEQSNRperChannel = subSNR;
   TEQoriginalChannelFreqResp = Fh;
   TEQequalizerFreqResp = Fw;
   TEQequalizedChannelFreqResp = Fhw;
   TEQchanNoisePowSpecAfterEqual = colorNoiseaft;
   TEQchanNoisePowSpecBeforeEqual = noiseSpec;
   TEQperformanceVersusDelay = Dv;
   TEQDmin = Dmin;
   TEQDmax = Dmax;
   TEQNb = Nb;
   TEQNw = Nw;
   TEQN = N;
   
   save teqresults TEQequalizerTaps TEQtargetImpulseResp  ...
	   TEQoptimalDelay TEQoriginalChannel TEQequalizedChannel ...
	   TEQmatchedFilterBoundPerChannel TEQSNRperChannel ...
	   TEQoriginalChannelFreqResp TEQequalizerFreqResp ...
	   TEQequalizedChannelFreqResp TEQchanNoisePowSpecAfterEqual...
	   TEQchanNoisePowSpecBeforeEqual TEQDmin TEQDmax TEQNb TEQNw TEQN;
	   
   % plot result
   teqdemo('plotGraph');
   set(gcf,'Pointer','arrow');
   return
   
   % plot graph
elseif strcmp(action,'plotGraph'),
  legend off
  set(gcf,'Pointer','watch');
   % get handles
   fhndlList = get(gcf,'Userdata');
   viewHndl = fhndlList(12);
   freqzHndl = fhndlList(1);
   % set axis
   axHndl = gca;
   axes(freqzHndl);
   % which graph is going to be ploted
   graphNum = get(viewHndl,'value'); 
   %1= TIR/SIR, 2= TEQ impulse res., 3= TEQ freq. res., 4= SNR/MFB
   %5=original and short channel impulse res, 6 = noise power spec., 
   %7= delay plot, 8= equalized channel freq
   
   % load results for current TEQ
   load tmpresults;
   freqAxis = 1:N/2+1;
   
   % plot chosen graph
   if graphNum == 1      
 	  legend off
      hw100 = hw(1:100)/norm(hw);
      axhndlList = get(freqzHndl,'UserData');
      hndl = plot(hw100); 
      hold on; 
      plot(D+1:D+Nb,B,'r','linewidth',2); hold off;
      set(gcf,'handlevisibility','callback')
      grid on
      xlabel('tap number')
      ylabel('amplitude')
      title([title_str,' TIR and SIR'])
      legend('SIR','TIR');
      
   elseif graphNum == 2
      legend off
      hndl = stem(W); 
      if Nw > 3
         hold on;
         plot(W,':');
         hold off
      end
      axis([0 Nw+1 1.1*min(W)*(min(W)<=0) 1.1*max(W)*(max(W)>0)])
      xlabel('tap number')
      ylabel('amplitude')
      title([title_str,' TEQ impulse response'])
      
   elseif graphNum == 3
      legend off
      plot(freqAxis(usedSubs),20*log10(abs(Fw(usedSubs))+eps));
      grid on
      axisSet = axis;
      axisSet(1) = min(freqAxis(usedSubs))-2;
      axisSet(2) = max(freqAxis(usedSubs))+2;
      axis(axisSet);
      xlabel('frequency bin number')
      ylabel('magnitude (dB)')
      title([title_str,' TEQ frequency response'])
      
   elseif graphNum == 4
      legend off
      plot(freqAxis(usedSubs),10*log10(subSNR(usedSubs)+eps))
      hold on;
      plot(freqAxis(usedSubs),10*log10(subMFB+eps),'r')
      hold off
      grid on
      axisSet = axis;
      axisSet(1) = min(freqAxis(usedSubs))-2;
      axisSet(2) = max(freqAxis(usedSubs))+2;
      axis(axisSet);
      xlabel('frequency')
      ylabel('magnitude (dB)')
      title([title_str,' achieved and upper-bound SNR'])
      legend('SNR','MFB',0);
      
   elseif graphNum == 5
      legend off
      maxAmp = max(abs(hw));
      plot(channel(1:500)/max(abs(channel))*maxAmp);
      hold on
      plot(hw(1:500),'r');
      hold off
      grid on	
      xlabel('discrete time (n)')
      ylabel('amplitude')
      title([title_str,' Original & Shortened Channel'])
      legend('original','shortened',0);         
      
   elseif graphNum == 6
      legend off
      plot(freqAxis(usedSubs),20*log10(noiseSpec+eps),'r')
      hold on
      plot(freqAxis(usedSubs),20*log10(colorNoiseaft(usedSubs)+eps))
      hold off
      grid on	
      axisSet = axis;
      axisSet(1) = min(freqAxis(usedSubs))-2;
      axisSet(2) = max(freqAxis(usedSubs))+2;
      axis(axisSet);
      xlabel('frequency')
      ylabel('magnitude (dB)')
      title([title_str,' Noise Power Spectrum'])
      legend('before TEQ','after TEQ',0);         
      
   elseif graphNum == 7
      legend off
      plot(Dmin:Dmax,10*log10(real(Dv(Dmin:Dmax))+eps))
      grid on	
      ylabel('objective function (e.g. MSE for MMSE method)')
      xlabel('delay')
      title([title_str,' Cost vs delay plot'])
      
   elseif graphNum == 8
      legend off
      plot(freqAxis(usedSubs),10*log10(abs(Fhw(usedSubs)).^2+eps));
      grid on
      axisSet = axis;
      axisSet(1) = min(freqAxis(usedSubs))-2;
      axisSet(2) = max(freqAxis(usedSubs))+2;
      axis(axisSet);
      ylabel('magnitude (dB)')
      xlabel('frequency')
      title([title_str,' Equalized Channel Frequency Response'])
   end
   
   set(gcf,'Pointer','arrow');
   return
   
   % get help
elseif strcmp(action,'info'),
   ttlStr = get(gcf,'name');
   myFig = gcf;
   
   topic1 =  ['Single-path TEQ design demo'];
   helptop1 = [...
     'This demo lets you design a TEQ and analyze its performance.          '
     '                                                                      '
     'Usage: Set parameters to desired values and hit the Calculate button. '
     '                                                                      '
     'The following design methods are supported:                           '
	 '  Minimum mean squared error -- unit energy constraint,               '		  
	 '  Minimum mean squared error -- unit tap constraint,                  ' 
     '  Maximum shortening signal to noise ratio method,                    '
     '  Maximum geometric SNR method,                                       '
     '  Divide and conquer -- cancellation,                                 '
     '  Divide and conquer -- minimization,                                 '
     '  Maximum channel capacity method,                                    '
     '  Minimum intersymbol interference method,                            '
     '  Matrix pencil design method,                                        '
     '  Modified matrix pencil design method.                               '
     '  ARMA modeling design method                                         '];
   
   topic2 =  ['Design and Simulation Paramters'];
   helptop2 = [...
	 'On the top of he control window on the right is a pulldown menu       '
     'from which a design method can be chosen.                             '
	 'Below this pulldown menu are the following editable text windows      '
     'which are used to set the desired parameters,                         '
     '  Shortened impulse response (SIR) length. This is the desired length '
     '   of the channel after equalization. For example, it should be set   '
     '   to 33 (one plus the cyclic prefix length) for the ADSL standard.   '
     '  Time domain equalizer (TEQ) length. Defines the number of taps of   ' 
	 '   the TEQ.                                                           '
     '  Fast Fourier transform (FFT) size. Sets the FFT size used in DMT    '
     '   modulation. It is twice the number of subchannels.                 '
     '  Coding gain (dB). Defines a coding in dB used during capacity       '
     '   calculations.                                                      '
     '  Margin (dB). Sets the desired system margin in dB. This is also     '
     '   used in capacity calculations.                                     '
     '  Dmin and Dmax. The optimal delay is search in the interval between  '
     '   Dmin and Dmax.                                                     '
	 '  Input power (dBm). Defines the input signal power in dBm.           '
	 '  AWGN pow (dBm/Hz). Sets the amount of additive white Gaussian noise '
     '   in dBm/Hz. AWGN is added to the near-end crosstalk noise.          '
     '  CSA loop \# (1-8). Selects the desired channel to run the simulation'
     '   on. Currently 8 standard CSA loops are supported.                  '];
   
     topic3 =  ['Graphics'];
   helptop3 = [...
	 'Below the editable text windows is another pull-down menu which       '
     'is used to select the desired graph to be displayed.                  '
     'The following graphics can be selected:                               '
     '  Target & shortened channel. Displays the shortened channel impulse  '	 
	 '   response and the target channel impulse response for the MMSE and  ' 
	 '   geometric SNR methods. For all other methods the location of the   ' 
     '   target window is displayed instead of a target impulse response.   '
     '  TEQ impulse response. Shows the impulse response of the TEQ.        '
     '  TEQ frequency response. Shows the TEQ in frequency domain.          '
     '  SNR & MFB. The SNR and matched filter bound to the SNR is displayed '
     '   as a function of frequency (subchannels).                          '
     '  Original & shortened channel. Displays the channel impulse response '
     '   before and after equalization.                                     '
     '  Noise Power Spectrum. Shows the power spectrum of the noise which   '
     '   consists of NEXT noise plus AWGN.                                  '
     '  Delay Plot. Displays the performance measure of the method with     '
     '   respect to the delay.                                              '
     '  Equalized channel freq resp. Displays the frequency response of the '
     '   channel after equalization.                                        '];
 
   topic4 =  ['Performance Measures'];
   helptop4 = [...
'The two remaining buttons in the control frame are                    '
'   Info. Displays this help.                                          '
'   Calculate. Starts the calculation and performance evaluation of the'
'    TEQ.                                                              '
'The following performance measures are calculated and listed in the   '
'table:                                                                '
'   Rate. Gives the achievable bit rate with the given channel and TEQ '
'    settings.                                                         '
'   SNR. Shows the SNR at the output of the equalizer in dB.           '
'   SSNR. Shows the shortening SNR in dB. This is defined as the ratio '
'    of the energy of the shortened channel impulse response in the    '
'    target window the the energy outside the target window.           '
'   MSE. Gives the MSE for the MMSE and geometric SNR methods.         '
'   Delay. Shows the optimal delay for the system.                     '
'   Max Rate. Shows the absolute maximum achievable bit rate given the '
'    channel and equalizer settings. It is basically the capacity      '
'   calculated from the MFB.                                           '];
		 
   str =  { topic1 helptop1; topic2 helptop2; topic3 helptop3; topic4 helptop4 };
   
   helpwin(str,'Topic 1','Signal-path TEQ design demo')                   
   return  
   
else
   disp(sprintf( 'Internal ERROR: ''%s'' not recognized',action))
end