Guest User

Untitled

a guest
Aug 8th, 2012
31
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. clear, clc    
  2.      
  3.     % Orders must be odd
  4.     order = 3;
  5.    
  6.     % *******************************************************************
  7.     % Set up signals
  8.     % *******************************************************************    
  9.        
  10.     fun = inline('sin(2*pi*x*3/40) + sin(2*pi*x*5/40)', 'x');
  11.    
  12.     % Input points and samples
  13.     nIn = 401;
  14.     xIn = 1:nIn;
  15.     inData = fun(xIn);
  16.    
  17.     % Output points and reference samples
  18.     ratio = 53/40;
  19.     xOut = 1 : 1/ratio : nIn;
  20.     interpolatorDelay = order/2-1;
  21.     ref = fun(xOut - interpolatorDelay);
  22.  
  23.     nOut = numel(xOut);
  24.     outData = zeros(1, nOut);
  25.  
  26.     % *******************************************************************
  27.     % Apply Lagrange interpolation to input data
  28.     % *******************************************************************    
  29.  
  30.     % Lagrange polynomial coefficients for each tap
  31.     for ixTap = 0:order
  32.         c(:,ixTap+1) = lagrangeBasisPolynomial(order, ixTap);
  33.     end
  34.  
  35.     % Calculate samples at the new points xOut
  36.     outData = [];
  37.     i=1;
  38.     for point = xOut
  39.  
  40.         % Calculate the fractional delay of 'point'
  41.         xOutFrac = point - floor(point);
  42.  
  43.         outTemp = 0;
  44.         for ixTap = 0:order
  45.             % inSample is the input sample that appears at FIR tap 'ixTap'
  46.             % to contribute to the output sample
  47.            
  48.             % Use end samples for the out-of-range samples
  49.             % inSample = mod(floor(point) + ixTap - order, nIn) + 1;
  50.            
  51.             % Just skip the out-of-range samples
  52.             inSample = floor(point) + ixTap - order + 1;
  53.  
  54.             if inSample < 1 || inSample > nIn
  55.                 continue;
  56.             end
  57.  
  58.             % The value of said input sample
  59.             d = inData( inSample );
  60.            
  61.             % Evaluate the Lagrange polynomial, resulting in the
  62.             % time-varying FIR tap weight
  63.             cTap = polyval(c(:,ixTap+1), xOutFrac);
  64.            
  65.             % FIR operation: multiply FIR tap weight with input sample and
  66.             % add to output sample
  67.             outTemp = outTemp + d.*cTap;
  68.         end
  69.                
  70.         outData(end+1) = outTemp;
  71.     end
  72.    
  73.     % *******************************************************************
  74.     % Trim ends to ignore filter response time
  75.     % *******************************************************************
  76.     trim = ceil(length(xOut) * 0.05); % trim first and last 5%
  77.     xOut = xOut(trim+1:end-trim);
  78.     outData = outData(trim+1:end-trim);
  79.     ref = ref(trim+1:end-trim);
  80.    
  81.     % *******************************************************************
  82.     % Plot
  83.     % *******************************************************************    
  84.     figure(1), clf, hold on
  85.     plot(xIn, inData, 'k+-');
  86.     plot(xOut, outData, 'b.-');
  87.     stem(xOut, ref, 'r');
  88.     legend('input', 'interpolation result', 'ideal result');
  89.     title('Lagrange interpolation, impulse response');
  90.     axis([25 45 -2 2])
  91.    
  92.     % *******************************************************************
  93.     % signal-to-noise ratio
  94.     % *******************************************************************    
  95.     SNR_dB = snr(ref, outData)
RAW Paste Data