Guest User

Untitled

a guest
Aug 3rd, 2012
29
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. % Lagrange interpolation for resampling
  2. % References:
  3. % [1] A digital signal processing approach to Interpolation
  4. %     Ronald W. Schafer and Lawrence R. Rabiner
  5. %     Proc. IEEE vol 61, pp.692-702, June 1973
  6. % [2] https://ccrma.stanford.edu/~jos/Interpolation/Lagrange_Interpolation.html
  7. % [3] Splitting the Unit delay: Tools for fractional delay filter design
  8. %     T. I. Laakso, V. Valimaki, M. Karjalainen, and U. K. Laine
  9. %     IEEE Signal Processing Magazine, vol. 13, no. 1, pp. 30..60, January 1996
  10. function lagrangeResamplingDemo()
  11.     originDefinition = 0; % see comment for LagrangeBasisPolynomial() below
  12.    
  13.     % Regarding order: From [1]
  14.     % "Indeed, it is easy to show that whenever Q is odd, none of the
  15.     % impulse responses corresponding to Lagrange interpolation can have
  16.     % linear phase"
  17.     % Here, order = Q+1 => odd orders are preferable
  18.     order = 5;
  19.    
  20.     % *******************************************************************
  21.     % Set up signals
  22.     % *******************************************************************    
  23.     fun = inline('sin(2*pi*x*3) + sin(2*pi*x*5)');
  24.     fs = 40; % sampling freq
  25.     x = 0:1/fs:10-1/fs;
  26.     inData = fun(x);
  27.  
  28.     nIn = length(inData);
  29.     nOut = 500;
  30.    
  31.     % Get Lagrange polynomial coefficients
  32.     c0 = lagrangeBasisPolynomial(order, 0, originDefinition);
  33.     c1 = lagrangeBasisPolynomial(order, 1, originDefinition);
  34.     c2 = lagrangeBasisPolynomial(order, 2, originDefinition);
  35.     c3 = lagrangeBasisPolynomial(order, 3, originDefinition);
  36.     c4 = lagrangeBasisPolynomial(order, 4, originDefinition);
  37.     c5 = lagrangeBasisPolynomial(order, 5, originDefinition);
  38.    
  39.     % *******************************************************************
  40.     % Plot
  41.     % *******************************************************************    
  42.     figure(2), clf, hold on  
  43.     xr = [0:10/nOut:10-10/nOut]; % fsn = 50Hz
  44.     yr = fun(xr);
  45.     plot(xr,yr,'g')
  46.    
  47.     % *******************************************************************
  48.     % Apply Lagrange interpolation to input data
  49.     % *******************************************************************  
  50.     outData = zeros(1,nOut);
  51.     for i=0:nOut-1
  52.        
  53.         % outTimeAtInputInteger
  54.         temp = floor( (i/nOut)*nIn );
  55.        
  56.         % the value of the input sample that appears at FIR tap 'ixTap'
  57.         d0 = inData( mod(temp + 0 - order, nIn) + 1 );
  58.         d1 = inData( mod(temp + 1 - order, nIn) + 1 );
  59.         d2 = inData( mod(temp + 2 - order, nIn) + 1 );
  60.         d3 = inData( mod(temp + 3 - order, nIn) + 1 );
  61.         d4 = inData( mod(temp + 4 - order, nIn) + 1 );
  62.         d5 = inData( mod(temp + 5 - order, nIn) + 1 );
  63.        
  64.         % Evaluate the Lagrange polynomials, resulting in the time-varying FIR tap weight
  65.         evalFracTime = -0.5 + mod(i,nOut/nIn) * (nIn/nOut);
  66.        
  67.         cTap0 = polyval(c0, evalFracTime);
  68.         cTap1 = polyval(c1, evalFracTime);
  69.         cTap2 = polyval(c2, evalFracTime);
  70.         cTap3 = polyval(c3, evalFracTime);
  71.         cTap4 = polyval(c4, evalFracTime);
  72.         cTap5 = polyval(c5, evalFracTime);
  73.        
  74.         % FIR operation: multiply FIR tap weight with input sample
  75.         outData(i+1) = d0*cTap0 + d1*cTap1 + d2*cTap2 + d3*cTap3 + d4*cTap4 + d5*cTap5;
  76.                
  77.         %plot( i/nOut, outData(i+1), 'b.');
  78.         %pause(0.01)
  79.     end
  80.  
  81.     plot( xr(1:end-4), outData(5:end), 'b.-');
  82.     legend('Reference', 'Interpolated')
  83.     axis([-0.1 1.2 -2 2])
  84.    
  85.     snr(yr(1:end-4), outData(5:end))
  86. end
  87.  
  88. % Returns the coefficients of a Lagrange basis polynomial
  89. % 1 <= order: Polynomial order
  90. % 0 <= evalIx <= order: index of basis function.
  91. %
  92. % At the set of support points, the basis polynomials evaluate as follows:
  93. % evalIx = 1: [1 0 0 ...]
  94. % evalIx = 2: [0 1 0 ...]
  95. % evalIx = 3: [0 0 1 ...]
  96. %
  97. % The support point are equally spaced.
  98. % Example, using originDefinition=0:
  99. % order = 1: [-0.5 0.5]
  100. % order = 2: [-1 0 1]
  101. % order = 3: [-1.5 -0.5 0.5 1.5]
  102. %
  103. % The area around the mid-point corresponds to -0.5 <= x <= 0.5.
  104. % If a resampler implementation uses by convention 0 <= x <= 1 instead, set
  105. % originDefinition=0.5 to offset
  106. % the polynomial.
  107. function polyCoeff = lagrangeBasisPolynomial(order, evalIx, originDefinition)
  108.     assert(evalIx >= 0);
  109.     assert(evalIx <= order);
  110.    
  111.     tapLocations = -0.5*(order) + (0:order) + originDefinition;
  112.  
  113.     polyCoeff = [1];
  114.     for loopIx = 0:order
  115.         if loopIx ~= evalIx        
  116.             % numerator: places a zero in the polynomial via (x-xTap(k)), with k != evalIx
  117.             % denominator: scales to 1 amplitude at x=xTap(evalIx)
  118.             term = [1 -tapLocations(loopIx+1)] / (tapLocations(evalIx+1)-tapLocations(loopIx+1));
  119.  
  120.             % multiply polynomials => convolve coefficients
  121.             polyCoeff = conv(polyCoeff, term);
  122.         end
  123.     end
  124.    
  125. end
RAW Paste Data