# 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