%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gc5016_filterdesign_0711
%
% filter generator for GC5016, 2 stage multirate DDC or DUC
% low pass filter, DDC and DUC 
% Texas Instruments - JQ 
% modified 2/28/2011 - to allow zeros at end of PFIR for repeater or low
% latency
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gc5016_filterdesign_0711( )
   
   % filter design for GC5016  
   % 4 channel computes both I and Q interleaved
   % split IQ has separate I and Q in parallel across 2 channels
   % ftype 1, DDC - clock rate 122.88Mhz, 4 chan dec2CIC, dec4PFIR - 8Mhz 
   % ftype 2,                             2 chan splitIQ  
   % ftype 3, DUC - clock rate 122.88Mhz, 4 chan int2PFIR, int4CIC - 8Mhz
   % ftype 4                              2 chan splitIQ    

   ftype = 4;

   outfolder = 'C:\Temp';
   DDCtapfile = 'gc5016_ddc_channel.taps';
   DUCtapfile = 'gc5016_duc_channel.taps';
   addtrailzeros = 0;

   if( ftype == 1)
   % ftype1 - 8Mhz BW, 12.5Msps IQ rate 
      clock_rate = 122.88e6;
      pfir_oversample_factor = 1;   
      pfir_ratio = 4;  % 4 channel mode limit 2 -> 8, splitIQ 1 -> 16
            
      cic_ratio = 2;
      cic_stages = 5;              % DDC
      cic_m = 1;
      symbol_rate = clock_rate / (cic_ratio * pfir_ratio * pfir_oversample_factor);

      cic_customcorr = 1;
      cic_corr_res = 100;
      ciccorr_numtaps = 5;  % value is typically 3, 5, if cic_customcorre == 1, must match size(CIC_sparray)
      CIC_sparray(1:ciccorr_numtaps) = [ .005, -.22, 1.57, -.22, .005];
   
      forcePFIRnum = 127; % 255 taps for splitIQ mode

      method = 1;    % 0 - LP, 1 - LP w/Window, 2 - RRC, 3- RRC w/Window, 4 - Hybrid, if PFIR taps > 160, 5 - manual file input

      Beta = 1;    % pulse shaping filter excess bandwidth (not used)
      hybrid_Fp = 1; % (not used)
      hybrid_Fs = 1; % (not used)
      
      windowtype = 1.002; % 0 - none, 1 - Hamming, >1 - KaiserBessel Alpha
      Fp = 4.1e6;        % used with Least Squares as end of passband
      Fs = 5.8e6;     % used with Least Squares as start of stopband
      Passbd = .05;   % used with Least Squares as passband ripple in db
      Stopbd = -90; % used with Least Squared as stopband rejection in db
      
      PFIRcoef_bits = 16; % fixed point resolution
      PFIRtap_infilename = 'GC5016ddc_pfirin.taps';
      PFIRtap_filename = 'GC5016ddc_pfir1.taps';
      Convtap_filename = DDCtapfile;
   end

   if( ftype == 2)
   % ftype1 - 8Mhz BW, 12.5Msps IQ rate 
      clock_rate = 122.88e6;
      pfir_oversample_factor = 1;   
      pfir_ratio = 4;  % 4 channel mode limit 2 -> 8, splitIQ 1 -> 16
            
      cic_ratio = 2;
      cic_stages = 5;              % DDC
      cic_m = 1;
      symbol_rate = clock_rate / (cic_ratio * pfir_ratio * pfir_oversample_factor);

      cic_customcorr = 1;
      cic_corr_res = 100;
      ciccorr_numtaps = 5;  % value is typically 3, 5, if cic_customcorre == 1, must match size(CIC_sparray)
      CIC_sparray(1:ciccorr_numtaps) = [ .005, -.22, 1.57, -.22, .005];
   
      forcePFIRnum = 255; % 255 taps for splitIQ mode

      method = 1;    % 0 - LP, 1 - LP w/Window, 2 - RRC, 3- RRC w/Window, 4 - Hybrid, if PFIR taps > 160, 5 - manual file input

      Beta = 1;    % pulse shaping filter excess bandwidth (not used)
      hybrid_Fp = 1; % (not used)
      hybrid_Fs = 1; % (not used)
      
      windowtype = 1.002; % 0 - none, 1 - Hamming, >1 - KaiserBessel Alpha
      Fp = 4.1e6;        % used with Least Squares as end of passband
      Fs = 5.4e6;     % used with Least Squares as start of stopband
      Passbd = .05;   % used with Least Squares as passband ripple in db
      Stopbd = -90; % used with Least Squared as stopband rejection in db
      
      PFIRcoef_bits = 16; % fixed point resolution
      PFIRtap_infilename = 'GC5016ddc_pfirin.taps';
      PFIRtap_filename = 'GC5016ddc_pfir2.taps';
      Convtap_filename = DDCtapfile;
   end
   

   if( ftype == 3)
   % ftype3 - DUC PFIR int2, CIC int4 - 8Mhz BW
      clock_rate = 122.88e6;
      pfir_oversample_factor = 1;   
      pfir_ratio = 2;  % 4 channel mode limit 2 -> 8, splitIQ 1 -> 16
            
      cic_ratio = 4;
      cic_stages = 6;              % DUC
      cic_m = 1;
      symbol_rate = clock_rate / (cic_ratio * pfir_ratio * pfir_oversample_factor);

      cic_customcorr = 1;
      cic_corr_res = 100;
      ciccorr_numtaps = 5;  % value is typically 3, 5, if cic_customcorre == 1, must match size(CIC_sparray)
      CIC_sparray(1:ciccorr_numtaps) = [ .006, -.25, 1.3, -.25, .006];
   
      totaltaps  = 127;
      addtrailzeros = 0;
      forcePFIRnum = totaltaps - addtrailzeros; 

      method = 1;    % 0 - LP, 1 - LP w/Window, 2 - RRC, 3- RRC w/Window, 4 - Hybrid, if PFIR taps > 160, 5 - manual file input

      Beta = 1;    % pulse shaping filter excess bandwidth (not used)
      hybrid_Fp = 1; % (not used)
      hybrid_Fs = 1; % (not used)
      
      windowtype = 1.002; % 0 - none, 1 - Hamming, >1 - KaiserBessel Alpha
      Fp = 4.1e6;        % used with Least Squares as end of passband
      Fs = 5.8e6;     % used with Least Squares as start of stopband
      Passbd = .05;   % used with Least Squares as passband ripple in db
      Stopbd = -90; % used with Least Squared as stopband rejection in db
      
      PFIRcoef_bits = 16; % fixed point resolution
      PFIRtap_infilename = 'GC5016duc_pfirin.taps';
      PFIRtap_filename = 'GC5016duc_pfir3.taps';
      Convtap_filename = DUCtapfile;
   end
  
   if( ftype == 4)
   % ftype4 - DUC PFIR int2, CIC int4 - 8Mhz BW
      clock_rate = 122.88e6;
      pfir_oversample_factor = 1;   
      pfir_ratio = 2;  % 4 channel mode limit 2 -> 8, splitIQ 1 -> 16
            
      cic_ratio = 4;
      cic_stages = 6;              % DUC
      cic_m = 1;
      symbol_rate = clock_rate / (cic_ratio * pfir_ratio * pfir_oversample_factor);

      cic_customcorr = 1;
      cic_corr_res = 100;
      ciccorr_numtaps = 5;  % value is typically 3, 5, if cic_customcorre == 1, must match size(CIC_sparray)
      CIC_sparray(1:ciccorr_numtaps) = [ .006, -.25, 1.3, -.25, .006];
   
      totaltaps  = 255;  %splitIQ mode
      addtrailzeros = 0;
      forcePFIRnum = totaltaps - addtrailzeros; 

      method = 1;    % 0 - LP, 1 - LP w/Window, 2 - RRC, 3- RRC w/Window, 4 - Hybrid, if PFIR taps > 160, 5 - manual file input

      Beta = 1;    % pulse shaping filter excess bandwidth (not used)
      hybrid_Fp = 1; % (not used)
      hybrid_Fs = 1; % (not used)
      
      windowtype = 1.002; % 0 - none, 1 - Hamming, >1 - KaiserBessel Alpha
      Fp = 4.1e6;        % used with Least Squares as end of passband
      Fs = 5.1e6;     % used with Least Squares as start of stopband
      Passbd = .05;   % used with Least Squares as passband ripple in db
      Stopbd = -90; % used with Least Squared as stopband rejection in db
      
      PFIRcoef_bits = 16; % fixed point resolution
      PFIRtap_infilename = 'GC5016duc_pfirin.taps';
      PFIRtap_filename = 'GC5016duc_pfir4.taps';
      Convtap_filename = DUCtapfile;
   end

   
   %%%%%%%%%%%%%%%%%%% start of code section %%%%%%%%%%%%%%%%%%%%%%%%
   % constants and fudge factors 
   Remez_fudge = 80;  % fudge factor for Remez versus Least Squares calc
   
   LPFilter = 0;
   LPFilter_Window = 1;
   RRC = 2;
   RRC_Window = 3;
   Hybrid = 4;
   fromfile = 5;

   PFIRtap_filename = strcat(outfolder, '\', PFIRtap_filename);
   Convtap_filename = strcat(outfolder, '\', Convtap_filename);
   CICtap_filename = strcat(outfolder, '\', 'cictaps.taps');
      
   clock_rate = cic_ratio * pfir_ratio * pfir_oversample_factor * symbol_rate;
   
   Firres_bits = PFIRcoef_bits;
   FIRmax = 2^ (Firres_bits - 1) - 1;
   
   HybridPFIR_StopRej = -85;
   
   numPFIRtaps = forcePFIRnum;
   
   CIC_Comparray(1) = [1];
   if( cic_ratio > 1) 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % determine number of CIC correction coefficients
   % if custom coefficients, load the CIC correction coefficients
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      if(cic_customcorr == 1)
          ciccorr_numtaps = max(size(CIC_sparray));
          clear CIC_Comparray;
          CIC_Comparray = CIC_sparray;
      end    
      numPFIRtaps = numPFIRtaps - ciccorr_numtaps + 1;
   end      
         
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   % determine PFIR clock rate
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   PFIRClk = clock_rate / cic_ratio;
   fracBWtap = (PFIRClk/(Fs - Fp))*255/64
      
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % normalize frequency and amplitude
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   freq(1:4) =    [0, Fp, Fs, (PFIRClk /2)];
   ampldbV(1:4) = [0, -.5 * Passbd, Stopbd, -150];

   % convert freq[] to normalized freq
   % convert amplitude in dbV to ratio, assumes same number of points as freq[]
   for n = 1: 4
      norm_freq(n) = freq(n) / (PFIRClk /2);  % normalize to filter samplerate
      if(n < 4)
         ratio_ampl(n) = 10^( ampldbV(n) / 20);
      else
         ratio_ampl(n) = 0;
      end   
   end
         
   % weighting vector uses passband ripple for quality factor
   % w(1: 2) = [1/Passbd, 10];
   w(1: 2) = [10, 2];

   % RRC filter parameters
   sampleratio = symbol_rate / PFIRClk;
      
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % run specific filter mode
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(method == LPFilter)
            
      % matlab wants odd symmetry, cmd5016 will add a 0 for non-symmetric mode
            
      % generate taps  
      if( numPFIRtaps < Remez_fudge)
         FIRtaps = remez( numPFIRtaps-1, norm_freq, ratio_ampl, w); % remez exchange   
      else    
         FIRtaps = firls( numPFIRtaps-1, norm_freq, ratio_ampl, w); % Least Square error
      end 
      bandwidth = 2 * Fp * 1.02;
   end 
      
   if(method == LPFilter_Window) 

      % generate taps  
      if( numPFIRtaps < Remez_fudge)
         FIRtaps = remez( numPFIRtaps-1, norm_freq, ratio_ampl, w); % remez exchange   
      else    
         FIRtaps = firls( numPFIRtaps-1, norm_freq, ratio_ampl, w); % Least Square error
      end  
 
      % windowing
      % 0 - none, 1 - Hamming, >1 - KaiserBessel Alpha
      flen = length(FIRtaps);
      w(1: flen) = 1;

      if(windowtype == 1)
         w = hamming(flen);                      
      end       
      if( windowtype > 1)
         w = kaiser( flen, windowtype);                      
      end    
      for n = 1: flen
         FIRtaps(n) = FIRtaps(n) * w(n);    
      end
      bandwidth = 2 * Fp * 1.02;
   end
   
   if(method == RRC)
      % generate closed form taps 
      FIRtaps = Calc_RRCTaps( sampleratio, Beta, FIRmax, numPFIRtaps);
      bandwidth = 2 * (1-Beta) * symbol_rate;
   end
      
   if( method == RRC_Window)
      % windowing excessBW compensation
      if(windowtype >= 1)
         Beta = Beta - (windowtype - 1);
      end
               
      % generate closed form taps 
      FIRtaps = Calc_RRCTaps( sampleratio, Beta, FIRmax, numPFIRtaps);
             
      % windowing
      % 0 - none, 1 - Hamming, >1 - KaiserBessel Alpha
               
      flen = length(FIRtaps);
      w(1: flen) = 1;
      if(windowtype == 1)
         w = hamming(flen);                      
      end       
      if(windowtype > 1)
         w = kaiser( flen, windowtype);                      
      end
      for n = 1: flen
         FIRtaps(n) = FIRtaps(n) * w(n);    
      end
      bandwidth = 2 * (1-Beta) * symbol_rate  * 1.02;
   end
   
   if(method == Hybrid) 
      % divide the number of taps between the RRC filter and the Low Pass Filter
      numRRCfil_taps = 32 * pfir_ratio * pfir_oversample_factor;

      numLPfil_taps  = numPFIRtaps - numRRCfil_taps + 1;
                  
      % generate closed form taps
      RRCtaps = Calc_RRCTaps( sampleratio, Beta, FIRmax, numRRCfil_taps);
            
      % generate the Low Pass filter mask Fp as 115% of RRC Fc, and the 1.2x of Fp for Fs
      symbolrate = PFIRClk * sampleratio;
            
      freq(1:4) =    [0, hybrid_Fp, hybrid_Fs, (PFIRClk /2)];
      ampldbV(1:4) = [0,  -.1, HybridPFIR_StopRej, -150]; 

      % convert freq[] to normalized freq
      % convert amplitude in dbV to ratio, assumes same number of points as freq[]
      for n = 1: 4
         norm_freq(n) = freq(n) / (PFIRClk /2);  % normalize to filter samplerate
         if(n < 4)
            ratio_ampl(n) = 10^( ampldbV(n) / 20);
         else
            ratio_ampl(n) = 0;
         end   
      end
      w(1: 2) = [10, 10];
             
      if( numPFIRtaps < Remez_fudge)
         LPtaps = remez( numLPfil_taps -1, norm_freq, ratio_ampl, w); % remez exchange   
      else    
         LPtaps = firls( numLPfil_taps -1, norm_freq, ratio_ampl, w); % Least Square error
      end 
            
      % combine the 32symbol RRC + Low Pass Stopband suppression
      FIRtaps = conv(LPtaps, RRCtaps);
      bandwidth = 2 * (1-Beta) * symbol_rate  * 1.02;
   end
  
   if(method == fromfile)
      filename = strcat(outfolder,'\', PFIRtap_infilename); 
      FIRtaps = read_filterdat( filename);
   end
   

   % either calculate the cic coefficients or generate them, them convolve
   % with FIRtaps
   cic_active = 0;
   if( cic_ratio > 1)
      if(cic_customcorr == 0)
         CIC_Comparray = ciccomp3or5taps(clock_rate, bandwidth, cic_stages, cic_m, cic_ratio, ciccorr_numtaps, cic_corr_res, 0, 0, 'auto_cic_corr', FIRtaps);
      end 
  
      % Convolve CIC Correction if needed
      FIRtaps = conv( FIRtaps, CIC_Comparray);
      cic_active = 1;
   end
         
   % round and scale
   sigstat(1:2) = [max(FIRtaps) (min(FIRtaps) * -1) ];
   FIRtaps = convrnd( FIRtaps * FIRmax / max(sigstat), Firres_bits );
   
   if(addtrailzeros > 0)
       firlen = length(FIRtaps);
       newtaplen = addtrailzeros + firlen;
       newtaps(1: newtaplen) = 0;
       newtaps(1: firlen) = FIRtaps(1: firlen);
       clear FIRtaps;
       FIRtaps = newtaps;
       clear newtaps;
   end    
      
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % plot filter results
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   flen = length(FIRtaps);
      
   X(1: flen) = 1: 1: flen;
   figure(1);
   plot(X, FIRtaps);
   titletext = sprintf('PFIR filter, %d taps', flen);
   title( titletext);
   zoom on;
      
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % generate equiv CIC taps, generate Convolved Response
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if( cic_active == 1) 
      CICtaps = GenCICTaps( cic_stages, cic_ratio, cic_m);
            
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % interpolate PFIR taps by CIC_ratio
      % Convolve Interp_PFIR with CIC
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      InterpPFIRtaps= zeros(1,length(FIRtaps) * cic_ratio);
      InterpPFIRtaps( 1 : cic_ratio : length(FIRtaps) * cic_ratio)= FIRtaps;
      Convtaps = conv(CICtaps, InterpPFIRtaps); 
      
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % generate equivalent PFIR taps by decimating by cic_ratio
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      dectaps = polyphase_dec_max(Convtaps, cic_ratio);

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % real FFT of taps
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      figure(2);
      titletext = sprintf('PFIR filter, convolved with CIC, and decimated back to PFIR rate');
      dectaps = scalereal( dectaps, FIRmax); 

      FFT64KNoFig( dectaps, PFIRClk, titletext); 
      zoom on;

   else
      CICtaps(1) = 1;
      Convtaps = FIRtaps;

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % real FFT of taps
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      figure(2);
      titletext = sprintf('PFIR filter');
      FFT64KNoFig( FIRtaps, PFIRClk, titletext); 
      zoom on;

   end   
            
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % plot and FFT Convolved Taps
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   figure(3);
   titletext = 'Convolved GC5016 channel taps, estimated delay';
   clen = length(Convtaps);
   ck_period = 1/ clock_rate;
   X(1: clen) = 0: ck_period: (clen - 1) * ck_period;
   plot( X, Convtaps);
   title( titletext);
   zoom on;
      
   figure(4);
   if(cic_ratio > 1)
      titletext = 'Convolved GC5016 channel CIC + PFIR';
   else
      titletext = 'Convolved GC5016 channel PFIR';
   end  
   FFTC64KNoWinNoFig( Convtaps, clock_rate, titletext);
   % FFT64KNoFig( Convtaps, clock_rate, titletext); 
   zoom on;
      
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % export PFIR coefficients to file
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   write_dcolumn(FIRtaps, PFIRtap_filename)
   write_filter( Convtap_filename, Convtaps);
  
return;
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write filter taps, single column vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function write_filter( filename, filtertaps)
   fid= fopen(filename,'w');
   DataLen = length(filtertaps);

   % print to file a line at a time
   for n = 1: DataLen
      fprintf( fid, '%12.5e \n', filtertaps(n) );
   end
   fclose(fid);
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computes Root Raised Cosine array, windows, scales
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function RRCtaps = Calc_RRCTaps( sampleratio, alpha, centertap, numtaps)
   
   titletext = sprintf('RRC FIR Response %d taps', numtaps);
   
   % RRC calculation
   % formula is:
   %   
   % tap(k) = [CT/((pi(1-A)+4A)]*
   %          [1/(1-(4AkR)**2)]* 
   %          [sin(kRpi(1-A))/kR + 4Acos(kRpi(1+A))]
   %    

   gain = centertap / ( pi*(1.0 - alpha) + (4.0 * alpha) );
 
   RRCtaps(1: numtaps) = 0.0;
   Fudge = 1.0e-10;

   totaltaps = 0.0;
   for k = 0 : (numtaps - 1)
      
      % second term 
      a = Fudge + sampleratio * (numtaps - (2*k) - 1.0)/2.0;
      t = 4.0 * alpha * a;
      tg = gain / (1.0-(t*t)+Fudge);
      s= sin(pi * (1.0 - alpha) * a);
      c = cos(pi * (1.0 + alpha) * a);
      
      % third term and combination
      RRCtaps(k+1) = tg * ((s/a) + (4.0*alpha*c));
   end  % next, completes loop calculation

   % re-scale for centertap value
   maxtap = max(RRCtaps);
   tapscale = centertap / maxtap;
   RRCtaps(1: numtaps) = RRCtaps(1: numtaps) * tapscale;

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% use GC5016 datasheet equations to determine ntaps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function numPFIRtaps = GC5016_CalcPFIRtaps( mode, cic_ratio, pfir_ratio, splitIQ, odd_symmetry) 
   numPFIRtaps = 0;
   % DDC mode
   if(mode == 1)
      if(cic_ratio < (16/pfir_ratio) )
         temp = 2 ^ fix( log2(cic_ratio / (2 - splitIQ) ) );
      else
         temp = 2 ^ fix( log2( (16/pfir_ratio) / (2 - splitIQ) ) );
      end    
      numPFIRtaps = (pfir_ratio * 16 * temp * (1 + odd_symmetry) ) - odd_symmetry;
   end    
   % DUC mode
   if(mode == 2)
      temp = 16 * pfir_ratio * (2 ^ fix( cic_ratio / (2 - splitIQ) ) ) ;
      if( temp > 256)
         temp = 256;
      end   
      numPFIRtaps = (temp * (odd_symmetry + 1) ) - odd_symmetry;
   end   
return;   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calc CIC Taps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CICtaps = GenCICTaps( cic_stages, cic_ratio, cic_m)

   % determine number of loops
   numConv = cic_stages - 1;
   
   % determine boxcar equiv
   boxtaps = ones((cic_ratio * cic_m), 1);
   
   % convolve n - (1,2) stages
   a = boxtaps;
   for n = 1 : numConv
      a = conv(a, boxtaps);
   end   
   
   CICtaps = a;
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% rescale complex array to max value
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Cscaled_array = scalecomplex( inputC_array, newmaxvalue)
   maxI = real( max( inputC_array) );
   minI = -1 * real( min( inputC_array) );
   maxQ = imag( max( inputC_array) );
   minQ = -1 * imag( min( inputC_array) );
   
   arraystat(1: 4) = [maxI, minI, maxQ, minQ];
   
   scalef = newmaxvalue / max(arraystat);
   Cscaled_array = scalef * inputC_array;
return; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RealFFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FFT64KNoFig( d, FClk, titletext)
   
   FFTsize = 65536;
   lend = length(d);
   b(1: FFTsize) = 0;                       % zero fill the array
   b(1: lend) = d(1: lend);

   F= fft(b, FFTsize);                      % FFT h(t)
   F= F(1: FFTsize/2);                       % We only need the first half of the FFT

   Xmax = FClk/2;
   X= linspace(0, Xmax, FFTsize/2);

   F= F + (realmin + realmin*i);             % This prevents log(0) warnings
   FF= 20*log10(abs(F));                     % Compute log of magnitude

   maxFF = max(FF);
   FF = FF - maxFF;                          % dBc 

   plot(X, FF);
   title( titletext);
   axis( [0 Xmax (max(FF)-120) max(FF)]); grid;            % Plot the data
return;

function rnd_array = convrnd(in_array, numbits)

   dlen = length(in_array);
   rnd_array(1: dlen) = 0;
   
   maxpos = 2^(numbits-1) - 1;
   maxneg = (-1 * maxpos) - 1;
   halfpos = 2^(numbits-2);
   halfneg = -1 * halfpos;
   pos_status = 0;
   neg_status = 0;
   
   for n = 1: dlen
      if (in_array(n) >= maxpos)
         rnd_array(n) = maxpos;
      elseif (in_array(n) <= maxneg)
         rnd_array(n) = maxneg;
      elseif (in_array(n) == halfpos)
         if(pos_status == 1)
             rnd_array(n) = fix(in_array(n));
             pos_status = 0;
         else
             rnd_array(n) = fix(in_array(n) + 0.5);
             pos_status = 1;
         end
      elseif (in_array(n) == halfneg)
          if(neg_status == 1)
             rnd_array(n) = fix(in_array(n));
             neg_status = 0;
         else
             rnd_array(n) = fix(in_array(n) + 0.5);
             neg_status = 1;
          end
      else
          rnd_array(n) = fix(in_array(n) + 0.5);
      end
   end   

return;

function FFTC64KNoWinNoFig( d, samplerate, titletext)

windowflag = 0;

lendata = length(d);               % size of data
sizeFFT = 2^( ( fix(log2(lendata)) ) + 1);       % set FFT size to next larger power of 2

if( sizeFFT < (2^16) )                     % at least 64K FFT 
   sizeFFT = (2^16);
end

if( lendata < sizeFFT)
   temp = lendata;
else
   temp = sizeFFT;
end   
b(1: sizeFFT) = 0 + i*0;
b(1: temp) = d(1: temp);  % put in input data

if( windowflag == 1)
   numbins = temp;                     % apply window 
   w(1: sizeFFT) = 0;
   w(1: numbins) = hamming(numbins); 
   
   bi(1: sizeFFT) = real(b(1: sizeFFT) );
   bq(1: sizeFFT) = imag(b(1: sizeFFT) );
   bi = bi .* w;
   bq = bq .* w;
   b = bi + (i*bq);
end

xmin = -.5 * samplerate;                   % plotted at the input rate, shifted for complex representation
xmax = .5 * samplerate;
X= linspace(xmin, xmax , sizeFFT);               % Create the X scale frequency axis

% b= finalTaps/sum(finalTaps);             % fraction of sum(? needed)
% F = fft(b, M);

F(1: sizeFFT) = 0 + i*0;
F = fft(b, sizeFFT);                    % FFT of convolved filter response
F= F + (realmin + realmin*i);             % This prevents log(0) warnings

FF= 20*log10(abs(F));                     % Compute log of magnitude

maxFF = max(FF);
FF = FF - maxFF;                          % dBc 

outFFT(1:sizeFFT) = 0 + 0*i;
hsizeFFT = fix(sizeFFT/2);

outFFT(1: hsizeFFT ) = FF( (hsizeFFT+1): sizeFFT);         % lower half of FFT plot is neg freq
outFFT( (hsizeFFT+1) : sizeFFT) = FF(1: hsizeFFT);          % upper half of FFT plot is pos freq

% figure(1);

plot( X, outFFT);
   
title( titletext);
axis( [xmin xmax -139 1]); 
grid;            % Plot the data

grid on;
zoom on;
return;

function dectaps = polyphase_dec_max(Convtaps, decimation_ratio) 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % decimate the convolved response of the PFIR and CIC back to the PFIR
   % rate, using the decimation phase, select the highest sum of abs(taps)
   % as the reference phase for the test
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   dlen = length(Convtaps);
   polyphaselen = fix(dlen / decimation_ratio);
   endarray = polyphaselen * decimation_ratio;
   
   % there are cic_interpolation polyphases, only one will be processed for
   % testing, the max sum(abs(taps)) 
   index = 1;
   maxsum = 0;
      
   for n = 1: decimation_ratio
      Decphase_array = CDecimate(Convtaps(1: endarray), decimation_ratio, n);
      dlen = length(Decphase_array);
      temp = 0;
      for m = 1: dlen
         temp = temp + abs( Decphase_array(m) );   
      end    

      if(temp > maxsum)
         index = n;
         maxsum = temp;
      end   
   end   
   dectaps = CDecimate(Convtaps(1: endarray), decimation_ratio, index);
return;

function [sq_err, sqrt_err] = mse(array1, array2)
   dlen1 = length(array1);
   dlen2 = length(array2);
   
   error = 0;
   if( dlen1 == dlen2)
      for n = 1: dlen1    
         error = error + (array1(n) - array2(n))^2; 
      end 
      sq_err = error / dlen1;
      sqrt_err = sqrt(error) / dlen1;
   else
      sq_err = -999;
      sqrt_err = -999;
   end
return; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CDecimate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CDecData = CDecimate(CData, decimation_ratio, decimation_phase) 
   
   %output data array size
   filtlength = length(real(CData));
   declength = fix(length(real(CData)) / decimation_ratio) - 1;
   
   % matrix decimation method?
   CDecData(1: declength) = CData( decimation_phase: decimation_ratio: ((declength-1) * decimation_ratio) + decimation_phase);
   clear CData;
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% rescale real array to max value
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function scaled_array = scalereal( input_array, newmaxvalue)
   maxI = max( input_array);
   minI = -1 * min( input_array);
   
   arraystat(1: 2) = [maxI, minI];
   
   scalef = newmaxvalue / max(arraystat);
   scaled_array = scalef * input_array;
return;

function sel_tap_array = LargestAbsSumTaps( Taparray, num_extractedtaps)
  % loop through the taps, looking for the largest abs sum
   testloops = length(Taparray) - num_extractedtaps;
   test_sum(1: testloops) = 0;
   testindex = 1;
   testmax = 0;
   for n = 1: testloops
      for m = 1: num_extractedtaps
         test_sum(n) = test_sum(n) + sqrt( (Taparray( n + m) )^2);       
      end
  
      if( test_sum(n) > testmax)
         testindex = n;
         testmax = test_sum(n);
      end   
   end
  
   % provide the subset indices
   start_index = testindex;
   stop_index = testindex + num_extractedtaps - 1;
   sel_tap_array(1: num_extractedtaps) =  Taparray(start_index: stop_index);
return;

function compfilter = ciccomp3or5taps(CICclockrate,bandwidth,CICstage,m,CICrate,numtaps,resolution,plotripple,savetaps,filepath, PFIRtaps);

    % CICclockrate = 160e6;    %Rate of PFIR clock, i.e. after interpolation or before decimation = CIC I/O rate
    % bandwidth = 10e6;           %Passband bandwidth
    % CICstage = 5;               %5 for DDC, 6 for DUC
    % m=1;                        %CIC delay
    % CICrate = 6;                %Interpolation/Decimation rate
    % numtaps = 5;                %number of compensation taps 3 or 5 taps
    % resolution = 200;           %resolution of search matrix
    % plotripple = 1;             %1 = plot ripple across bandwidth
    % savetaps = 1;               %1 = save text file of taps
    % filepath='C:\Documents and Settings\A0188431\My Documents\My GC Studio Projects\GC5016\TSW4100\DDC Experiments\GC5016_RX_CIC6_PFIR2_4CH\GC5016_RX_CIC6_PFIR2_4CH_10MHz\CICcomp.txt';
    % PFIRtaps = []               % array of PFIR taps
    
    PFIRclockrate = CICclockrate/CICrate;
    
    cictaps=cicgeneratorlocal(CICstage,CICrate,m);
    
    firlen = max(size(PFIRtaps));
    InterpPFIRtaps= zeros(1, firlen * CICrate);
    InterpPFIRtaps( 1 : CICrate : firlen * CICrate)= PFIRtaps;
    Convtaps = conv(cictaps, InterpPFIRtaps); 

    BW = bandwidth/(2*PFIRclockrate*CICrate);

    % xmin=.1;
    xmin=.02;
    % xmax=.35;
    xmax=.4;
    xresolution=resolution;

    if numtaps == 5
        ymin=.002;
        ymax=.1;
        yresolution=resolution;
    elseif numtaps == 3
        ymin=0;
        ymax=0;
        yresolution=1;  
    else
        sprintf('numtaps should be 3 or 5')
    end

    numfft=round(25/BW);

    BWminmax=zeros(xresolution,yresolution);

    numtaps_int = CICrate*5 - 1; 

    BWpoints=round(BW*numfft)+1;

    for xindex = 1:xresolution
        for yindex = 1:yresolution
        xtap = -(xmin+(xmax-xmin)*xindex/xresolution);
        ytap = (ymin+(ymax-ymin)*yindex/yresolution);
        comparray=[ytap xtap 1 xtap ytap];
        comparray_int=zeros(numtaps_int,1);
        comparray_int(1:CICrate:numtaps_int)=comparray;
        % convolved=conv(comparray_int,cictaps);
        convolved=conv(comparray_int,Convtaps);
        convolved_fft=fftextendlocal(convolved,numfft);
        BWarray=20*log10(abs(convolved_fft(1:BWpoints)));
        BWminmax(xindex,yindex)=max(BWarray)-min(BWarray);
        end
    end

    [minxindex,minyindex]=find(BWminmax==min(min(BWminmax)));
    xtap = -(xmin+(xmax-xmin)*minxindex/xresolution);
    ytap = (ymin+(ymax-ymin)*minyindex/yresolution);
    comparray=[ytap xtap 1 xtap ytap];
    comparray_int=zeros(numtaps_int,1);
    comparray_int(1:CICrate:numtaps_int)=comparray;
    
    if plotripple == 1
        % convolved=conv(comparray_int,cictaps);
        convolved=conv(comparray_int,Convtaps);
        convolved_fft=fftextendlocal(convolved,numfft);
        BWarray=20*log10(abs(convolved_fft(1:BWpoints)));
        x=(0:BWpoints-1)*bandwidth/((BWpoints-1)*2);
        figure(1);plot(x,BWarray-mean(BWarray));
        grid on;
        xlabel('Frequency (Hz)');
        ylabel('Ripple (dB)')
    end

    ripple=BWminmax(minxindex,minyindex)
    comparray = comparray_int(1:CICrate:length(comparray_int));
    
    if numtaps == 3
        comparray(5)=[];
        comparray(1)=[];
    end
    
    comparray
    
    if savetaps == 1
        filepath
        write_column(comparray, filepath);
    end

    compfilter = comparray;

return;

function cicfilter = cicgeneratorlocal(stages,rate,delay)
    a=ones(1,rate*delay);
    b=a;
    for index = 1:(stages-1)
        b=conv(b,a);
    end
    cicfilter = b/sum(b);
return;

 function outfft=fftextendlocal(inarray,numfft);
     inarray_extend=zeros(numfft,1);
     inarray_extend(1:length(inarray))=inarray;
     outfft=fft(inarray_extend);
 return;
 
 function write_column(data, outfilename)
   fout = fopen(outfilename,'wt');      %open a file
   lend = length(data);
   for k=1: lend
      fprintf(fout,'%g \r',data(k) );
   end
   fclose(fout);
 return;
 
  function write_dcolumn(data, outfilename)
   fout = fopen(outfilename,'wt');      %open a file
   lend = length(data);
   for k=1: lend
      fprintf(fout,'%d \n',round(data(k)) );
   end
   fclose(fout);
 return;

 function taps = read_filterdat( filename)

   fid= fopen(filename,'r');

   taps = fscanf(fid, '%g');
   fclose(fid);
return;