calc_gc1115coef.m

一个滤波器设计程序

   end    
       
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % 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) * carrierBW;
   end
      
   if( method == RRC_Window)
      % windowing excessBW compensation
      if(windowtype > 1)
         Beta = Beta - (windowtype - 1);
      else
         if(windowtype == 1)
            Beta = Beta - .025;
         end
      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) * carrierBW  * 1.02;
   end

   if(chopmode > 0)
      dlen = max(size(FIRtaps));
      numtaps = dlen - (2 * chopmode);
      cFIRtaps(1: numtaps) = FIRtaps(chopmode: dlen - chopmode - 1);
      clear FIRtaps;
      FIRtaps = cFIRtaps;
   end   
   
   % round and scale
   sigstat(1:2) = [max(FIRtaps) (min(FIRtaps) * -1) ];
   FIRtaps = convrnd( FIRtaps * FIRmax / max(sigstat), Firres_bits );
      
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % plot filter results
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   flen = length(FIRtaps);
      
   X(1: flen) = 1: 1: flen;
   figure(1);
   plot(X, FIRtaps);
   titletext = sprintf('GC1115 TimeSeries Cancellation filter, %d taps', flen);
   title( titletext);
   zoom on;

   Convtaps = FIRtaps;

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % real FFT of taps
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   figure(2);
   if(method == 0)
      titletext = sprintf('GC1115 cancellation Low Pass filter');
   end   
   if(method == 1)
      titletext = sprintf('GC1115 cancellation Low Pass with window filter');
   end   
   if(method == 2)
      titletext = sprintf('GC1115 cancellation Root Raised Cosine filter');
   end   
   if(method == 3)
      titletext = sprintf('GC1115 cancellation Root Raised Cosine with window filter');
   end   
   FFT64KNoFig( FIRtaps, clock_rate, titletext); 
   zoom on;

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % export PFIR coefficients to file
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   write_filter( strcat(outcoef_foldername, '\',outcoef_filename) , FIRtaps);

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write filter taps, single column vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function write_filter( filename, filtertaps)
   fid= fopen(filename,'wt');
   DataLen = length(filtertaps);

   % print to file a line at a time
   for n = 1: DataLen
      fprintf( fid, '%d \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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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 FFTdata = 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;