# 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