blob: c639ca256860bf8666fb4babdcd8e96c7a4ff5fb (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
|
function d = pchips(x,y,delta)
//Piecewise Cubic Hermite Interpolating Polynomial (PCHIP)
//Parameters
// x: a vector
// y: is Y is vector then it must have the same length as x and Y is matrix then the last dimension of Y must equal
//length(X).
// delta: Points for interpolation
// d: vector of interpolantant at delta
//Examples:
//x = -3:3;
//y = [-1 -1 -1 0 1 1 1];
//xq1 = -3:.01:3;
//v=pchips(x,y,xq1)
//
n = length(x);
if n==2
d = repmat(delta(1),size(y));
else
d = zeros(size(y));
k = find(sign(delta(1:n-2)).*sign(delta(2:n-1)) > 0);
h = diff(x);
hs = h(k)+h(k+1);
w1 = (h(k)+hs)./(3*hs);
w2 = (hs+h(k+1))./(3*hs);
if ~isempty (k) then
del_mx = max(abs(delta(k)), abs(delta(k+1)));
del_mn = min(abs(delta(k)), abs(delta(k+1)));
d(k+1) = del_mn./conj(w1.*(delta(k)./del_mx) + w2.*(delta(k+1)./del_mx));
else
d(1) = ((2*h(1)+h(2))*delta(1) - h(1)*delta(2))/(h(1)+h(2));
if sign(d(1)) ~= sign(delta(1))
d(1) = 0;
elseif (sign(delta(1)) ~= sign(delta(2))) & (abs(d(1)) > abs(3*delta(1)))
d(1) = 3*delta(1);
end
d(n) = ((2*h(n-1)+h(n-2))*delta(n-1) - h(n-1)*delta(n-2))/(h(n-1)+h(n-2));
if sign(d(n)) ~= sign(delta(n-1))
d(n) = 0;
elseif (sign(delta(n-1)) ~= sign(delta(n-2))) & (abs(d(n)) > abs(3*delta(n-1)))
d(n) = 3*delta(n-1);
end
end
end
endfunction
|
