dsp_main.m

PHS多天线基带调制与解调MATLAB程序

function [results] = DSP_Main(desc, env_in, alg_in) 
global env alg sys;
global n;

env = struct(...
          'nBursts', 500, ...             
          'SNR', 40, ...                  
          'nInt', 1, ...          % support 1 interfer for now!        
          'CIR', 20, ...                  
          'intDelay', 10, ...             
          'freqDrift', 0, ...             
          'freqOffset', 1500, ...           
          'timingDrift', 1, ...         
          'timingOffset', 2, ...            
          'SO_I', 1, ...                   
          'spatSig', 'rayleigh', ...       
          'dopplerS', 7, ...               
          'dopplerI', 7, ...               
          'refSigBer', 0 ...              
           );

alg = struct(...
          'mode', 'MRC', ...                   
          'FreqAlignEstMode', 'normal', ...    
          ...                            
          ...                            
          'nFreqEstSym', 80, ...          
          'beginFreqEstSym', 20, ...      
          'nRzzSym', 75, ...              
          'beginRzzSym', 30, ...          
          'alphaScale', 0, ...          
          'nIterations', 0, ...           
          'nCma', 2, ...                  
          'feedback', 1, ...              
          'UWWeights', 1 ...             
          ); 

sys = struct(...
          'frame', 200, ...               
          'over', 8, ...                   
          'rate', 192e3, ...              
          'nBit', 240, ...                
          'nSym', 120, ...                
          'nSamp', 960, ...               
          'nAnt', 4, ...                  
          'tracking', 0.5, ...            
          'keyS', 501, ...                
          'keyI', 547 ...                 
          );

%--- Initialize data structure -----------------------------------------------% 

 
  data = struct(...
          'fo_S', [], ...                 
          'filt_fo_S', [], ...            
          'new_fo_S', [], ...             
          'to_S', [], ...                 
          'filt_to_S', [], ...            
          'new_to_S', [], ...             
          'fo_I', [], ...                 
          'to_I', [], ...                 
          'be_S', [], ...                  
                      ...                 
          'be_I', [], ...                 
                      ...                 
          'w', [], ...                    
	       'a_I', [], ... 
          'w_src', [], ...                
          'w_int', [], ...                
          'w_ideal', [], ...              
          'w_mrc', [], ...                
          'mrc_src',[], ...               
          'mrc_int',[], ...               
          'mrc_ideal', [], ...            
          'src_int', [], ...                   
          'src_ideal', [], ...            
          'int_ideal', [], ...            
          'lock_detect', [], ...           
          'trackstate', [] ...		  
          );

  global FAD_S;                 
  global FAD_I;                 

  global FO_S;                  
  global FO_I;                  

  global TO_S;                  
  global TO_I;                  

  global N_TCH_S;               
  global TCH_S_ARRAY;           

  global N_TCH_I;               
  global TCH_I_ARRAY;           

  global N_NOISE;               
  global NOISE;                 

  global ENV_PREV;              
  global ALG_PREV;              

  global Z

  timestamp = date;
  user = getenv('USER');

  readme = [];
  readme = [readme sprintf('Started for %s on %s\n', user, timestamp)];

  readme = [readme sprintf('Description: %s\n', desc)];

  readme = [readme sprintf('\nENV structure modifications:\n')];
  if ((nargin >= 2) & ~isempty(env_in))
    names = fieldnames(env_in);
    for i=1:length(names),
      
      if (~isfield(env, names{i}))
        error(sprintf('env field (%s) is invalid', names{i}));
      end;

      new = getfield(env_in, names{i});
      old = getfield(env, names{i});
      if (ischar(new))
        readme = [readme sprintf('%s: %s --> %s\n', names{i}, old, new)];
      else 
        readme = [readme sprintf('%s: %f --> %f\n', names{i}, old, new)];
      end

      env = setfield(env, names{i}, new);
    end
  else
    readme = [readme sprintf('None.\n')];
  end

  readme = [readme sprintf('\nALG structure modifications:\n')];
  if ((nargin >= 3) & ~isempty(alg_in))
    names = fieldnames(alg_in);
    for i=1:length(names)
      
      if (~isfield(alg, names{i}))
        error(sprintf('alg field (%s) is invalid', names{i}));
      end

      new = getfield(alg_in, names{i});
      old = getfield(alg, names{i});
      if (ischar(new))
        readme = [readme sprintf('%s: %s --> %s\n', names{i}, old, new)];
      else 
        readme = [readme sprintf('%s: %f --> %f\n', names{i}, old, new)];
      end

      alg = setfield(alg, names{i}, new);
    end
  else
    readme = [readme sprintf('None.\n')];
 end
 

  j = sqrt(-1);

  start_seed = 1;

  maxFreqOffset = 5e3;

  maxFreqDrift = 50;
  
  minPrime = 23;

% primes = nextPrimes(max(env.nBursts,minPrime^2),3);  
NNN=max(env.nBursts,minPrime^2);
nnn=3;
x = ceil(sqrt(NNN)) ;
primes = [] ;
  k = x ;

  while length(primes) < nnn

     r = rem(k,[2:floor(sqrt(k))]) ;
     index = find(r == 0) ;
     if length(index) == 0
        primes = [ primes k ] ;
     end
     k = k + 1 ;

  end



noiseSymShift = 40;

  rand('seed',start_seed);  
  start_seed = start_seed + 1;
  
  randn('seed',start_seed); 
  start_seed = start_seed + 1;

  gen_fad = 0;
  if (isempty(ENV_PREV))
    gen_fad = 1;
  elseif (ENV_PREV.nBursts < env.nBursts)
    gen_fad = 1;
  elseif (~strcmp(ENV_PREV.spatSig, env.spatSig))
    gen_fad = 1;
  elseif (ENV_PREV.dopplerS ~= env.dopplerS)
    gen_fad = 1;
  elseif (ENV_PREV.dopplerI ~= env.dopplerI)
    gen_fad = 1;
  elseif (ENV_PREV.nInt ~= env.nInt)
    gen_fad = 1;
  end;


  if (gen_fad)
    if (strcmp (env.spatSig, 'rayleigh'))
      fprintf(1, 'Building Rayleigh fading coeffs...\n');

      nFFT = 2^(ceil(log(env.nBursts)/log(2))) ;

      FAD_S = [];
      FAD_I = [];


      for n = 1:sys.nAnt
        r = rayleigh(env.dopplerS,sys.frame,nFFT,env.nBursts) ;
        FAD_S = [ FAD_S ; r ; ] ;
      end

      for i = 1:env.nInt
        for n = 1:sys.nAnt
          r = rayleigh(env.dopplerI,sys.frame,nFFT,env.nBursts) ;
          FAD_I = [ FAD_I ; r];
        end
      end
      
    else
      fprintf(1, 'Building constant amplitude spatial signatures...\n');

      theta_S = 2*pi*rand(env.nBursts, sys.nAnt)*env.dopplerS/sys.frame;
      theta_S = cumsum(theta_S)';

      theta_I = 2*pi*rand(env.nBursts, sys.nAnt)*env.dopplerI/sys.frame;
      theta_I = cumsum(theta_I)';

      FAD_S = exp(j*theta_S);
      FAD_I = exp(j*theta_I);
    end
 end
 
  rand('seed',start_seed);  
  start_seed = start_seed + 1;
  
  randn('seed',start_seed); 
  start_seed = start_seed + 1;
  
  gen_freq = 0;
  if (isempty(ENV_PREV))
    gen_freq = 1;
  elseif (ENV_PREV.nBursts < env.nBursts)
    gen_freq = 1;
  elseif (ENV_PREV.nInt ~= env.nInt)
    gen_freq = 1;
  elseif (ENV_PREV.freqDrift ~= env.freqDrift)
    gen_freq = 1;
  elseif (ENV_PREV.timingDrift ~= env.timingDrift)
    gen_freq = 1;
  elseif ((env.freqDrift == 0) & (ENV_PREV.freqOffset ~= env.freqOffset))
    gen_freq = 1;
  elseif ((env.timingDrift == 0) & (ENV_PREV.timingOffset ~= env.timingOffset))
    gen_freq = 1;
  end

  if (gen_freq)
    fprintf(1, 'Generating frequency/timing offsets between Source and Interferer...\n');

    FO_S = [];
    FO_I = [];

    TO_S = [];
    TO_I = [];

    if (env.freqDrift ~= 0)

        FO_S(1) = min([maxFreqOffset, (env.nBursts-1)*env.freqDrift]);
        freqDriftNew = (FO_S(1)/(env.nBursts-1));
        
        if (freqDriftNew ~= env.freqDrift)
          readme = [readme sprintf('Scaling freqDrift from %f to %f\n', ...
                               env.freqDrift, freqDriftNew)];
        end

        FO_I(1) = max([-maxFreqOffset, -(env.nBursts-1)*freqDriftNew]);

        for n = 2:env.nBursts

	    FO_S(n) = FO_S(1) - (n-1)*freqDriftNew + ((2*rand-1)*maxFreqDrift/2);  
	    FO_I(n) = FO_I(1) + (n-1)*freqDriftNew + ((2*rand-1)*maxFreqDrift/2);

	end
    else
	FO_S = zeros(1,env.nBursts);
	FO_I = env.freqOffset * ones(1,env.nBursts);
    end;

    if (env.timingDrift ~= 0)

	TO_S(1) = sys.over * rand;
	TO_I(1) = sys.over * rand;
	for n = 2:env.nBursts

	  TO_S(n) = TO_S(n-1) + env.timingDrift * 2 * (rand-0.5) ;
	  TO_I(n) = TO_I(n-1) + env.timingDrift * 2 * (rand-0.5) ;

	end


	TO_S = mod (round(TO_S), sys.over);
	TO_I = mod (round(TO_I), sys.over);
    else
	TO_S = zeros(1, env.nBursts);
	TO_I = env.timingOffset * ones(1, env.nBursts);
    end;
  end


  H = RC_filt_mex(0,8,8) ;
  nH = length(H) ;

  h = RC_filt_mex(1,8,8) ;
  nh = length(h) ;
  
  rand('seed',start_seed);  
  start_seed = start_seed + 1;
  
  randn('seed',start_seed); 
  start_seed = start_seed + 1; 

  gen_tch = 0;
  if (isempty(ENV_PREV))
    gen_tch = 1;
  elseif (ENV_PREV.nBursts < env.nBursts)
    gen_tch = 1;
  elseif (ENV_PREV.nInt ~= env.nInt)
    gen_tch = 1;
  elseif (ENV_PREV.refSigBer ~= env.refSigBer)
    gen_tch = 1;
  end

  if (gen_tch)
    fprintf(1,'Generating random TCH bursts...\n');

    N_TCH_S = primes(1);
    N_TCH_I = primes(2);
    
    for n=1:N_TCH_S
      [TCH_S_ARRAY(n).s_S,TCH_S_ARRAY(n).r_S,TCH_S_ARRAY(n).TCH_S,TCH_S_ARRAY(n).ref_S] = ...
      buildTCHburst(h,H,sys.over,sys.keyS,env.refSigBer) ;

    end

    if (env.nInt ~= 0)
      for n=1:N_TCH_I
        [TCH_I_ARRAY(n).s_I,TCH_I_ARRAY(n).r_I,TCH_I_ARRAY(n).TCH_I,TCH_I_ARRAY(n).ref_I] = ...
        buildTCHburst(h,H,sys.over,sys.keyI,env.refSigBer) ;
      end
    end
  end

  rand('seed',start_seed);  
  start_seed = start_seed + 1;
  
  randn('seed',start_seed); 
  start_seed = start_seed + 1;


  gen_noise = 0;
  if (isempty(ENV_PREV))
    gen_noise = 1;
  elseif (ENV_PREV.nBursts < env.nBursts)
    gen_noise = 1;
  end

  if (gen_noise)
    fprintf(1,'Generating filtered AWGN channel noise...\n');

    N_NOISE = primes(3);

    noise_len = sys.nSamp + (N_NOISE-1)*noiseSymShift*sys.over;

    NOISE = randn(noise_len,sys.nAnt) + j * randn(noise_len,sys.nAnt);
    NOISE = 1 / sqrt(2) * NOISE;
    
    NOISE = conv2(NOISE,h.','same');

    NOISE = NOISE / sqrt(mean(mean(NOISE .* conj(NOISE))));
  end

  results.readme = readme;
  results.env = env;
  results.alg = alg;
  results.sys = sys;

%------------------------------------------------------------------

  
  w = ones(sys.nAnt,1)/sqrt(sys.nAnt);

  fprintf(1,'\nSNR: %3d CIR: %3d, Burst: %5.5d,', env.SNR, env.CIR, 0);
  for n = 1:env.nBursts
    if( mod(n,10) == 0 )
        fprintf(1,'\b\b\b\b\b%5.5d,', env.nBursts-n);
    end

      
    idx_S = rem(n-1, N_TCH_S) + 1;
    idx_I = rem(n-1, N_TCH_I) + 1;
    idx_N = rem(n-1, N_NOISE) + 1;

    a_S = FAD_S(:,n).' ;   

    if (env.nInt == 1)
      a_I = FAD_I(:,n).' ; 
    end

    
    s_S = TCH_S_ARRAY(idx_S).s_S;
    r_S = TCH_S_ARRAY(idx_S).r_S;
    TCH_S = TCH_S_ARRAY(idx_S).TCH_S;
    perfect_ref_S = TCH_S_ARRAY(idx_S).ref_S;

    if (env.nInt == 1)
      s_I = TCH_I_ARRAY(idx_I).s_I;
      r_I = TCH_I_ARRAY(idx_I).r_I;
      TCH_I = TCH_I_ARRAY(idx_I).TCH_I;
      perfect_ref_I = TCH_I_ARRAY(idx_I).ref_I;
    end

    cisoid_S = exp(j*2*pi*[1:sys.over*sys.nSym]*FO_S(n)/(sys.over*sys.rate)).' ;
    cisoid_I = exp(j*2*pi*[1:sys.over*sys.nSym]*FO_I(n)/(sys.over*sys.rate)).' ;

    Z_S = 10^(env.SNR/20) * (cisoid_S .* r_S) * a_S ;

    if (env.nInt == 1)
      Z_I = 10^((env.SNR-env.CIR)/20) * (cisoid_I .* r_I) * a_I ;
    else
      Z_I = zeros(size(Z_S));
    end

    if (n <= env.intDelay)
      scale = 10^(-((env.intDelay-n)*(40)/(env.intDelay-1))/20);
      Z_I = Z_I * scale;
    end

    lower = (idx_N-1)*noiseSymShift*sys.over+1;
    upper = lower + sys.nSamp -1;
    Z_N = NOISE(lower:upper,:);

    Z = [zeros(TO_S(n), sys.nAnt);Z_S(1:sys.nSamp-TO_S(n), :)] + ...
        [zeros(TO_I(n)+8*env.SO_I, sys.nAnt);Z_I(1:sys.nSamp-TO_I(n)-8*env.SO_I, :)] + ...    
    	Z_N;

%        [zeros(TO_I(n)+8*3, sys.nAnt);0.1*Z_S(1:sys.nSamp-TO_I(n)-8*3, :)] + ...
    

    if( n == 1 )
	filt_to_S(n) = TO_S(n);
        filt_fo_S(n) = FO_S(n);
    else
        filt_to_S(n) = mod(filt_to_S(n-1) + 1/16 * (new_to_S - filt_to_S(n-1)),sys.over);
        filt_fo_S(n) = filt_fo_S(n-1) + 1/32 * (new_fo_S - filt_fo_S(n-1));
%        filt_fo_S(n) = filt_fo_S(n-1) + 1/32 * new_fo_S;
    end

    BitParam.PerfectRef = perfect_ref_S(1+alg.beginRzzSym : ...
                          alg.beginRzzSym+alg.nRzzSym);
    BitParam.uwref = perfect_ref_S(8:17);
    if( n == 1)
        AirParam.w = [0 0 0 0]'; 
        BitParam.UsePerfectRef = 1;
    else
        BitParam.UsePerfectRef = 0;
	AirParam.w = new_w;
    end

    if( strcmp(alg.FreqAlignEstMode,'perfect'))
        AirParam.TO = TO_S(n);
        AirParam.FO = FO_S(n);
    else
        AirParam.TO = filt_to_S(n);
        AirParam.FO = filt_fo_S(n);
    end


    switch alg.mode

      case 'MRC'
		  %disp('MRC Mode ...');
        [AirData, BitData] = MRC(AirParam);
        mrc = AirData.w;
        trackstate = 0 ;
     
      case 'ModifiedWiener'
		  %disp('ModifiedWiener Mode ...');
        [AirData, BitData] = ModifiedWiener(AirParam, BitParam);
        trackstate = 1 ;
        
      otherwise

        error (sprintf('Don''t recognize %s mode of operation', env.mode));

    end
   

    if (env.nInt == 0)
      w_ideal = a_S';
    elseif (env.nInt == 1)
      w_ideal = (a_S - (a_S * (a_I'/norm(a_I)) * (a_I/norm(a_I))))';
    end
   

    if length(BitData.S_info) == 160
      data.be_S(n) = sum(xor(BitData.S_info,TCH_S.')) ;
    else
      data.be_S(n) = -1 ;
    end

    if length(BitData.I_info) == 160
      data.be_I(n) = sum(xor(BitData.I_info,TCH_I.')) ;
    else
      data.be_I(n) = -1 ;
    end

    new_w = AirData.w; 
    new_to_S = AirData.TO;
    new_fo_S = AirData.FO;
   

    AirData.w = new_w;
    AirParam = AirData;
   

    data.w(:,n) = AirData.w;
    data.a_I(:,n) = a_I;


  end

  ENV_PREV = env;
  ALG_PREV = alg;

  results.data = data;

  fprintf(1, '\nAll done with SNR: %3d, CIR: %3d\n\n', env.SNR, env.CIR);

return;