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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
|
<!DOCTYPE html>
<html>
<!-- Created by GNU Texinfo 7.1.1, https://www.gnu.org/software/texinfo/ -->
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<title>Polynomial Interpolation (GNU Octave (version 10.3.0))</title>
<meta name="description" content="Polynomial Interpolation (GNU Octave (version 10.3.0))">
<meta name="keywords" content="Polynomial Interpolation (GNU Octave (version 10.3.0))">
<meta name="resource-type" content="document">
<meta name="distribution" content="global">
<meta name="Generator" content="makeinfo">
<meta name="viewport" content="width=device-width,initial-scale=1">
<link href="index.html" rel="start" title="Top">
<link href="Concept-Index.html" rel="index" title="Concept Index">
<link href="index.html#SEC_Contents" rel="contents" title="Table of Contents">
<link href="Polynomial-Manipulations.html" rel="up" title="Polynomial Manipulations">
<link href="Miscellaneous-Functions.html" rel="next" title="Miscellaneous Functions">
<link href="Derivatives-_002f-Integrals-_002f-Transforms.html" rel="prev" title="Derivatives / Integrals / Transforms">
<style type="text/css">
<!--
a.copiable-link {visibility: hidden; text-decoration: none; line-height: 0em}
div.center {text-align:center}
div.example {margin-left: 3.2em}
span:hover a.copiable-link {visibility: visible}
strong.def-name {font-family: monospace; font-weight: bold; font-size: larger}
-->
</style>
<link rel="stylesheet" type="text/css" href="octave.css">
</head>
<body lang="en">
<div class="section-level-extent" id="Polynomial-Interpolation">
<div class="nav-panel">
<p>
Next: <a href="Miscellaneous-Functions.html" accesskey="n" rel="next">Miscellaneous Functions</a>, Previous: <a href="Derivatives-_002f-Integrals-_002f-Transforms.html" accesskey="p" rel="prev">Derivatives / Integrals / Transforms</a>, Up: <a href="Polynomial-Manipulations.html" accesskey="u" rel="up">Polynomial Manipulations</a> [<a href="index.html#SEC_Contents" title="Table of contents" rel="contents">Contents</a>][<a href="Concept-Index.html" title="Index" rel="index">Index</a>]</p>
</div>
<hr>
<h3 class="section" id="Polynomial-Interpolation-1"><span>28.5 Polynomial Interpolation<a class="copiable-link" href="#Polynomial-Interpolation-1"> ¶</a></span></h3>
<p>Octave comes with good support for various kinds of interpolation,
most of which are described in <a class="ref" href="Interpolation.html">Interpolation</a>. One simple alternative
to the functions described in the aforementioned chapter, is to fit
a single polynomial, or a piecewise polynomial (spline) to some given
data points. To avoid a highly fluctuating polynomial, one most often
wants to fit a low-order polynomial to data. This usually means that it
is necessary to fit the polynomial in a least-squares sense, which just
is what the <code class="code">polyfit</code> function does.
</p>
<a class="anchor" id="XREFpolyfit"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-polyfit"><span><code class="def-type"><var class="var">p</var> =</code> <strong class="def-name">polyfit</strong> <code class="def-code-arguments">(<var class="var">x</var>, <var class="var">y</var>, <var class="var">n</var>)</code><a class="copiable-link" href="#index-polyfit"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-polyfit-1"><span><code class="def-type">[<var class="var">p</var>, <var class="var">s</var>] =</code> <strong class="def-name">polyfit</strong> <code class="def-code-arguments">(<var class="var">x</var>, <var class="var">y</var>, <var class="var">n</var>)</code><a class="copiable-link" href="#index-polyfit-1"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-polyfit-2"><span><code class="def-type">[<var class="var">p</var>, <var class="var">s</var>, <var class="var">mu</var>] =</code> <strong class="def-name">polyfit</strong> <code class="def-code-arguments">(<var class="var">x</var>, <var class="var">y</var>, <var class="var">n</var>)</code><a class="copiable-link" href="#index-polyfit-2"> ¶</a></span></dt>
<dd><p>Return the coefficients of a polynomial <var class="var">p</var>(<var class="var">x</var>) of degree <var class="var">n</var>
that minimizes the least-squares-error of the fit to the points
<code class="code">[<var class="var">x</var>(:), <var class="var">y</var>(:)]</code>.
</p>
<p><var class="var">n</var> is typically an integer ≥ 0 specifying the degree of the
approximating polynomial. If <var class="var">n</var> is a logical vector, it is used as a
mask to selectively force the corresponding polynomial coefficients to be
used or ignored.
</p>
<p>The polynomial coefficients are returned in the row vector <var class="var">p</var>. The
output <var class="var">p</var> may be directly used with <code class="code">polyval</code> to estimate values
using the fitted polynomial.
</p>
<p>The optional output <var class="var">s</var> is a structure containing the following fields:
</p>
<dl class="table">
<dt>‘<samp class="samp">yf</samp>’</dt>
<dd><p>The values of the polynomial for each value of <var class="var">x</var>.
</p>
</dd>
<dt>‘<samp class="samp">X</samp>’</dt>
<dd><p>The Vandermonde matrix used to compute the polynomial
coefficients.
</p>
</dd>
<dt>‘<samp class="samp">R</samp>’</dt>
<dd><p>Triangular factor R from the QR decomposition.
</p>
</dd>
<dt>‘<samp class="samp">C</samp>’</dt>
<dd><p>The unscaled covariance matrix, formally equal to the inverse of
<var class="var">x’</var>*<var class="var">x</var>, but computed in a way minimizing roundoff error
propagation.
</p>
</dd>
<dt>‘<samp class="samp">df</samp>’</dt>
<dd><p>The degrees of freedom.
</p>
</dd>
<dt>‘<samp class="samp">normr</samp>’</dt>
<dd><p>The norm of the residuals.
</p></dd>
</dl>
<p>The second output may be used by <code class="code">polyval</code> to calculate the statistical
error limits of the predicted values. In particular, the standard deviation
of <var class="var">p</var> coefficients is given by
</p>
<p><code class="code">sqrt (diag (<var class="var">s.C</var>)/<var class="var">s.df</var>) * <var class="var">s.normr</var></code>.
</p>
<p>When the third output, <var class="var">mu</var>, is present the original data is centered
and scaled which can improve the numerical stability of the fit. The
coefficients <var class="var">p</var> are associated with a polynomial in
</p>
<p><code class="code"><var class="var">xhat</var> = (<var class="var">x</var> - <var class="var">mu</var>(1)) / <var class="var">mu</var>(2)</code> <br>
where <var class="var">mu</var>(1) = mean (<var class="var">x</var>), and <var class="var">mu</var>(2) = std (<var class="var">x</var>).
</p>
<p>Example 1 : logical <var class="var">n</var> and integer <var class="var">n</var>
</p>
<div class="example">
<div class="group"><pre class="example-preformatted">f = @(x) x.^2 + 5; # data-generating function
x = 0:5;
y = f (x);
## Fit data to polynomial A*x^3 + B*x^1
p = polyfit (x, y, logical ([1, 0, 1, 0]))
⇒ p = [ 0.0680, 0, 4.2444, 0 ]
## Fit data to polynomial using all terms up to x^3
p = polyfit (x, y, 3)
⇒ p = [ -4.9608e-17, 1.0000e+00, -3.3813e-15, 5.0000e+00 ]
</pre></div></div>
<p><strong class="strong">See also:</strong> <a class="ref" href="Evaluating-Polynomials.html#XREFpolyval">polyval</a>, <a class="ref" href="Derivatives-_002f-Integrals-_002f-Transforms.html#XREFpolyaffine">polyaffine</a>, <a class="ref" href="Finding-Roots.html#XREFroots">roots</a>, <a class="ref" href="Famous-Matrices.html#XREFvander">vander</a>, <a class="ref" href="Basic-Statistical-Functions.html#XREFzscore">zscore</a>.
</p></dd></dl>
<p>In situations where a single polynomial isn’t good enough, a solution
is to use several polynomials pieced together. The function
<code class="code">splinefit</code> fits a piecewise polynomial (spline) to a set of
data.
</p>
<a class="anchor" id="XREFsplinefit"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-splinefit"><span><code class="def-type"><var class="var">pp</var> =</code> <strong class="def-name">splinefit</strong> <code class="def-code-arguments">(<var class="var">x</var>, <var class="var">y</var>, <var class="var">breaks</var>)</code><a class="copiable-link" href="#index-splinefit"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-splinefit-1"><span><code class="def-type"><var class="var">pp</var> =</code> <strong class="def-name">splinefit</strong> <code class="def-code-arguments">(<var class="var">x</var>, <var class="var">y</var>, <var class="var">p</var>)</code><a class="copiable-link" href="#index-splinefit-1"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-splinefit-2"><span><code class="def-type"><var class="var">pp</var> =</code> <strong class="def-name">splinefit</strong> <code class="def-code-arguments">(…, "periodic", <var class="var">periodic</var>)</code><a class="copiable-link" href="#index-splinefit-2"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-splinefit-3"><span><code class="def-type"><var class="var">pp</var> =</code> <strong class="def-name">splinefit</strong> <code class="def-code-arguments">(…, "robust", <var class="var">robust</var>)</code><a class="copiable-link" href="#index-splinefit-3"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-splinefit-4"><span><code class="def-type"><var class="var">pp</var> =</code> <strong class="def-name">splinefit</strong> <code class="def-code-arguments">(…, "beta", <var class="var">beta</var>)</code><a class="copiable-link" href="#index-splinefit-4"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-splinefit-5"><span><code class="def-type"><var class="var">pp</var> =</code> <strong class="def-name">splinefit</strong> <code class="def-code-arguments">(…, "order", <var class="var">order</var>)</code><a class="copiable-link" href="#index-splinefit-5"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-splinefit-6"><span><code class="def-type"><var class="var">pp</var> =</code> <strong class="def-name">splinefit</strong> <code class="def-code-arguments">(…, "constraints", <var class="var">constraints</var>)</code><a class="copiable-link" href="#index-splinefit-6"> ¶</a></span></dt>
<dd>
<p>Fit a piecewise cubic spline with breaks (knots) <var class="var">breaks</var> to the
noisy data, <var class="var">x</var> and <var class="var">y</var>.
</p>
<p><var class="var">x</var> is a vector, and <var class="var">y</var> is a vector or N-D array. If <var class="var">y</var> is an
N-D array, then <var class="var">x</var>(j) is matched to <var class="var">y</var>(:,…,:,j).
</p>
<p><var class="var">p</var> is a positive integer defining the number of intervals along
<var class="var">x</var>, and <var class="var">p</var>+1 is the number of breaks. The number of points in
each interval differ by no more than 1.
</p>
<p>The optional property <var class="var">periodic</var> is a logical value which specifies
whether a periodic boundary condition is applied to the spline. The
length of the period is <code class="code">max (<var class="var">breaks</var>) - min (<var class="var">breaks</var>)</code>.
The default value is <code class="code">false</code>.
</p>
<p>The optional property <var class="var">robust</var> is a logical value which specifies
if robust fitting is to be applied to reduce the influence of outlying
data points. Three iterations of weighted least squares are performed.
Weights are computed from previous residuals. The sensitivity of outlier
identification is controlled by the property <var class="var">beta</var>. The value of
<var class="var">beta</var> is restricted to the range, 0 < <var class="var">beta</var> < 1. The default
value is <var class="var">beta</var> = 1/2. Values close to 0 give all data equal
weighting. Increasing values of <var class="var">beta</var> reduce the influence of
outlying data. Values close to unity may cause instability or rank
deficiency.
</p>
<p>The fitted spline is returned as a piecewise polynomial, <var class="var">pp</var>, and
may be evaluated using <code class="code">ppval</code>.
</p>
<p>The splines are constructed of polynomials with degree <var class="var">order</var>.
The default is a cubic, <var class="var">order</var>=3. A spline with P pieces has
P+<var class="var">order</var> degrees of freedom. With periodic boundary conditions
the degrees of freedom are reduced to P.
</p>
<p>The optional property, <var class="var">constraints</var>, is a structure specifying linear
constraints on the fit. The structure has three fields, <code class="code">"xc"</code>,
<code class="code">"yc"</code>, and <code class="code">"cc"</code>.
</p>
<dl class="table">
<dt><code class="code">"xc"</code></dt>
<dd><p>Vector of the x-locations of the constraints.
</p>
</dd>
<dt><code class="code">"yc"</code></dt>
<dd><p>Constraining values at the locations <var class="var">xc</var>.
The default is an array of zeros.
</p>
</dd>
<dt><code class="code">"cc"</code></dt>
<dd><p>Coefficients (matrix). The default is an array of ones. The number of
rows is limited to the order of the piecewise polynomials, <var class="var">order</var>.
</p></dd>
</dl>
<p>Constraints are linear combinations of derivatives of order 0 to
<var class="var">order</var>-1 according to
</p>
<div class="example">
<div class="group"><pre class="example-preformatted">cc(1,j) * y(xc(j)) + cc(2,j) * y'(xc(j)) + ... = yc(:,...,:,j).
</pre></div></div>
<p><strong class="strong">See also:</strong> <a class="ref" href="One_002ddimensional-Interpolation.html#XREFinterp1">interp1</a>, <a class="ref" href="#XREFunmkpp">unmkpp</a>, <a class="ref" href="#XREFppval">ppval</a>, <a class="ref" href="One_002ddimensional-Interpolation.html#XREFspline">spline</a>, <a class="ref" href="Signal-Processing.html#XREFpchip">pchip</a>, <a class="ref" href="#XREFppder">ppder</a>, <a class="ref" href="#XREFppint">ppint</a>, <a class="ref" href="#XREFppjumps">ppjumps</a>.
</p></dd></dl>
<p>The number of <var class="var">breaks</var> (or knots) used to construct the piecewise
polynomial is a significant factor in suppressing the noise present in
the input data, <var class="var">x</var> and <var class="var">y</var>. This is demonstrated by the example
below.
</p>
<div class="example">
<div class="group"><pre class="example-preformatted">x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.2 * randn (size (x));
## Uniform breaks
breaks = linspace (0, 2 * pi, 41); % 41 breaks, 40 pieces
pp1 = splinefit (x, y, breaks);
## Breaks interpolated from data
pp2 = splinefit (x, y, 10); % 11 breaks, 10 pieces
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "41 breaks, 40 pieces", "11 breaks, 10 pieces"})
</pre></div></div>
<p>The result of which can be seen in <a class="ref" href="#fig_003asplinefit1">Figure 28.1</a>.
</p>
<div class="float" id="fig_003asplinefit1">
<div class="center"><img class="image" src="splinefit1.png" alt="splinefit1">
</div><div class="caption"><p><strong class="strong">Figure 28.1: </strong>Comparison of a fitting a piecewise polynomial with 41 breaks to one
with 11 breaks. The fit with the large number of breaks exhibits a fast ripple
that is not present in the underlying function.</p></div></div>
<p>The piecewise polynomial fit, provided by <code class="code">splinefit</code>, has
continuous derivatives up to the <var class="var">order</var>-1. For example, a cubic fit
has continuous first and second derivatives. This is demonstrated by
the code
</p>
<div class="example">
<pre class="example-preformatted">## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.1 * randn (size (x));
## Piecewise constant
pp1 = splinefit (x, y, 8, "order", 0);
## Piecewise linear
pp2 = splinefit (x, y, 8, "order", 1);
## Piecewise quadratic
pp3 = splinefit (x, y, 8, "order", 2);
## Piecewise cubic
pp4 = splinefit (x, y, 8, "order", 3);
## Piecewise quartic
pp5 = splinefit (x, y, 8, "order", 4);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
y4 = ppval (pp4, xx);
y5 = ppval (pp5, xx);
plot (x, y, ".", xx, [y1; y2; y3; y4; y5])
axis tight
ylim auto
legend ({"data", "order 0", "order 1", "order 2", "order 3", "order 4"})
</pre></div>
<p>The result of which can be seen in <a class="ref" href="#fig_003asplinefit2">Figure 28.2</a>.
</p>
<div class="float" id="fig_003asplinefit2">
<div class="center"><img class="image" src="splinefit2.png" alt="splinefit2">
</div><div class="caption"><p><strong class="strong">Figure 28.2: </strong>Comparison of a piecewise constant, linear, quadratic, cubic, and
quartic polynomials with 8 breaks to noisy data. The higher order solutions
more accurately represent the underlying function, but come with the
expense of computational complexity.</p></div></div>
<p>When the underlying function to provide a fit to is periodic, <code class="code">splinefit</code>
is able to apply the boundary conditions needed to manifest a periodic fit.
This is demonstrated by the code below.
</p>
<div class="example">
<div class="group"><pre class="example-preformatted">## Data (100 points)
x = 2 * pi * [0, (rand (1, 98)), 1];
y = sin (x) - cos (2 * x) + 0.2 * randn (size (x));
## No constraints
pp1 = splinefit (x, y, 10, "order", 5);
## Periodic boundaries
pp2 = splinefit (x, y, 10, "order", 5, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "no constraints", "periodic"})
</pre></div></div>
<p>The result of which can be seen in <a class="ref" href="#fig_003asplinefit3">Figure 28.3</a>.
</p>
<div class="float" id="fig_003asplinefit3">
<div class="center"><img class="image" src="splinefit3.png" alt="splinefit3">
</div><div class="caption"><p><strong class="strong">Figure 28.3: </strong>Comparison of piecewise polynomial fits to a noisy periodic
function with, and without, periodic boundary conditions.</p></div></div>
<p>More complex constraints may be added as well. For example, the code below
illustrates a periodic fit with values that have been clamped at the endpoints,
and a second periodic fit which is hinged at the endpoints.
</p>
<div class="example">
<pre class="example-preformatted">## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (2 * x) + 0.1 * randn (size (x));
## Breaks
breaks = linspace (0, 2 * pi, 10);
## Clamped endpoints, y = y' = 0
xc = [0, 0, 2*pi, 2*pi];
cc = [(eye (2)), (eye (2))];
con = struct ("xc", xc, "cc", cc);
pp1 = splinefit (x, y, breaks, "constraints", con);
## Hinged periodic endpoints, y = 0
con = struct ("xc", 0);
pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "clamped", "hinged periodic"})
</pre></div>
<p>The result of which can be seen in <a class="ref" href="#fig_003asplinefit4">Figure 28.4</a>.
</p>
<div class="float" id="fig_003asplinefit4">
<div class="center"><img class="image" src="splinefit4.png" alt="splinefit4">
</div><div class="caption"><p><strong class="strong">Figure 28.4: </strong>Comparison of two periodic piecewise cubic fits to a noisy periodic
signal. One fit has its endpoints clamped and the second has its endpoints
hinged.</p></div></div>
<p>The <code class="code">splinefit</code> function also provides the convenience of a <var class="var">robust</var>
fitting, where the effect of outlying data is reduced. In the example below,
three different fits are provided. Two with differing levels of outlier
suppression and a third illustrating the non-robust solution.
</p>
<div class="example">
<pre class="example-preformatted">## Data
x = linspace (0, 2*pi, 200);
y = sin (x) + sin (2 * x) + 0.05 * randn (size (x));
## Add outliers
x = [x, linspace(0,2*pi,60)];
y = [y, -ones(1,60)];
## Fit splines with hinged conditions
con = struct ("xc", [0, 2*pi]);
## Robust fitting, beta = 0.25
pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25);
## Robust fitting, beta = 0.75
pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75);
## No robust fitting
pp3 = splinefit (x, y, 8, "constraints", con);
## Plot
xx = linspace (0, 2*pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
plot (x, y, ".", xx, [y1; y2; y3])
legend ({"data with outliers","robust, beta = 0.25", ...
"robust, beta = 0.75", "no robust fitting"})
axis tight
ylim auto
</pre></div>
<p>The result of which can be seen in <a class="ref" href="#fig_003asplinefit6">Figure 28.5</a>.
</p>
<div class="float" id="fig_003asplinefit6">
<div class="center"><img class="image" src="splinefit6.png" alt="splinefit6">
</div><div class="caption"><p><strong class="strong">Figure 28.5: </strong>Comparison of two different levels of robust fitting (<var class="var">beta</var> = 0.25 and 0.75) to noisy data combined with outlying data. A conventional fit, without
robust fitting (<var class="var">beta</var> = 0) is also included.</p></div></div>
<p>A very specific form of polynomial interpretation is the Padé approximant.
For control systems, a continuous-time delay can be modeled very simply with
the approximant.
</p>
<a class="anchor" id="XREFpadecoef"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-padecoef"><span><code class="def-type">[<var class="var">num</var>, <var class="var">den</var>] =</code> <strong class="def-name">padecoef</strong> <code class="def-code-arguments">(<var class="var">T</var>)</code><a class="copiable-link" href="#index-padecoef"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-padecoef-1"><span><code class="def-type">[<var class="var">num</var>, <var class="var">den</var>] =</code> <strong class="def-name">padecoef</strong> <code class="def-code-arguments">(<var class="var">T</var>, <var class="var">N</var>)</code><a class="copiable-link" href="#index-padecoef-1"> ¶</a></span></dt>
<dd><p>Compute the <var class="var">N</var>th-order Padé approximant of the continuous-time
delay <var class="var">T</var> in transfer function form.
The Padé approximant of <code class="code">exp (-sT)</code> is defined by the
following equation
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"> Pn(s)
exp (-sT) ~ -------
Qn(s)
</pre></div></div>
<p>Where both Pn(s) and Qn(s) are <var class="var">N</var>th-order rational functions
defined by the following expressions
</p>
<div class="example">
<div class="group"><pre class="example-preformatted"> N (2N - k)!N! k
Pn(s) = SUM --------------- (-sT)
k=0 (2N)!k!(N - k)!
Qn(s) = Pn(-s)
</pre></div></div>
<p>The inputs <var class="var">T</var> and <var class="var">N</var> must be non-negative numeric scalars. If
<var class="var">N</var> is unspecified it defaults to 1.
</p>
<p>The output row vectors <var class="var">num</var> and <var class="var">den</var> contain the numerator and
denominator coefficients in descending powers of s. Both are
<var class="var">N</var>th-order polynomials.
</p>
<p>For example:
</p>
<div class="example smallexample">
<div class="group"><pre class="example-preformatted">t = 0.1;
n = 4;
[num, den] = padecoef (t, n)
⇒ num =
1.0000e-04 -2.0000e-02 1.8000e+00 -8.4000e+01 1.6800e+03
⇒ den =
1.0000e-04 2.0000e-02 1.8000e+00 8.4000e+01 1.6800e+03
</pre></div></div>
</dd></dl>
<p>The function, <code class="code">ppval</code>, evaluates the piecewise polynomials, created
by <code class="code">mkpp</code> or other means, and <code class="code">unmkpp</code> returns detailed
information about the piecewise polynomial.
</p>
<p>The following example shows how to combine two linear functions and a
quadratic into one function. Each of these functions is expressed
on adjoined intervals.
</p>
<div class="example">
<div class="group"><pre class="example-preformatted">x = [-2, -1, 1, 2];
p = [ 0, 1, 0;
1, -2, 1;
0, -1, 1 ];
pp = mkpp (x, p);
xi = linspace (-2, 2, 50);
yi = ppval (pp, xi);
plot (xi, yi);
</pre></div></div>
<a class="anchor" id="XREFmkpp"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-mkpp"><span><code class="def-type"><var class="var">pp</var> =</code> <strong class="def-name">mkpp</strong> <code class="def-code-arguments">(<var class="var">breaks</var>, <var class="var">coefs</var>)</code><a class="copiable-link" href="#index-mkpp"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-mkpp-1"><span><code class="def-type"><var class="var">pp</var> =</code> <strong class="def-name">mkpp</strong> <code class="def-code-arguments">(<var class="var">breaks</var>, <var class="var">coefs</var>, <var class="var">d</var>)</code><a class="copiable-link" href="#index-mkpp-1"> ¶</a></span></dt>
<dd>
<p>Construct a piecewise polynomial (pp) structure from sample points
<var class="var">breaks</var> and coefficients <var class="var">coefs</var>.
</p>
<p><var class="var">breaks</var> must be a vector of strictly increasing values. The number of
intervals is given by <code class="code"><var class="var">ni</var> = length (<var class="var">breaks</var>) - 1</code>.
</p>
<p>When <var class="var">m</var> is the polynomial order <var class="var">coefs</var> must be of size:
<var class="var">ni</var>-by-(<var class="var">m</var> + 1)<!-- /@w -->.
</p>
<p>The i-th row of <var class="var">coefs</var>, <code class="code"><var class="var">coefs</var>(<var class="var">i</var>,:)</code>, contains the
coefficients for the polynomial over the <var class="var">i</var>-th interval, ordered from
highest (<var class="var">m</var>) to lowest (<var class="var">0</var>) degree.
</p>
<p><var class="var">coefs</var> may also be a multi-dimensional array, specifying a
vector-valued or array-valued polynomial. In that case the polynomial
order <var class="var">m</var> is defined by the length of the last dimension of <var class="var">coefs</var>.
The size of first dimension(s) are given by the scalar or vector <var class="var">d</var>.
If <var class="var">d</var> is not given it is set to <code class="code">1</code>. In this case
<code class="code"><var class="var">p</var>(<var class="var">r</var>, <var class="var">i</var>, :)</code> contains the coefficients for the
<var class="var">r</var>-th polynomial defined on interval <var class="var">i</var>. In any case <var class="var">coefs</var>
is reshaped to a 2-D matrix of size <code class="code">[<var class="var">ni</var>*prod(<var class="var">d</var>) <var class="var">m</var>]</code>.
</p>
<p>Programming Note: <code class="code">ppval</code> evaluates polynomials at
<code class="code"><var class="var">xi</var> - <var class="var">breaks</var>(i)</code>, i.e., it subtracts the lower endpoint of
the current interval from <var class="var">xi</var>. This must be taken into account when
creating piecewise polynomials objects with <code class="code">mkpp</code>.
</p>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFunmkpp">unmkpp</a>, <a class="ref" href="#XREFppval">ppval</a>, <a class="ref" href="One_002ddimensional-Interpolation.html#XREFspline">spline</a>, <a class="ref" href="Signal-Processing.html#XREFpchip">pchip</a>, <a class="ref" href="#XREFppder">ppder</a>, <a class="ref" href="#XREFppint">ppint</a>, <a class="ref" href="#XREFppjumps">ppjumps</a>.
</p></dd></dl>
<a class="anchor" id="XREFunmkpp"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-unmkpp"><span><code class="def-type">[<var class="var">x</var>, <var class="var">p</var>, <var class="var">n</var>, <var class="var">k</var>, <var class="var">d</var>] =</code> <strong class="def-name">unmkpp</strong> <code class="def-code-arguments">(<var class="var">pp</var>)</code><a class="copiable-link" href="#index-unmkpp"> ¶</a></span></dt>
<dd>
<p>Extract the components of a piecewise polynomial structure <var class="var">pp</var>.
</p>
<p>This function is the inverse of <code class="code">mkpp</code>: it extracts the inputs to
<code class="code">mkpp</code> needed to create the piecewise polynomial structure <var class="var">pp</var>.
The code below makes this relation explicit:
</p>
<div class="example">
<div class="group"><pre class="example-preformatted">[breaks, coefs, numinter, order, dim] = unmkpp (pp);
pp2 = mkpp (breaks, coefs, dim);
</pre></div></div>
<p>The piecewise polynomial structure <code class="code">pp2</code> obtained in this way, is
identical to the original <code class="code">pp</code>. The same can be obtained by directly
accessing the fields of the structure <code class="code">pp</code>.
</p>
<p>The components are:
</p>
<dl class="table">
<dt><var class="var">x</var></dt>
<dd><p>Sample points or breaks.
</p>
</dd>
<dt><var class="var">p</var></dt>
<dd><p>Polynomial coefficients for points in sample interval.
<code class="code"><var class="var">p</var>(<var class="var">i</var>, :)</code> contains the coefficients for the polynomial
over interval <var class="var">i</var> ordered from highest to lowest degree.
If <code class="code"><var class="var">d</var> > 1</code>, then <var class="var">p</var> is a matrix of size
<code class="code">[<var class="var">n</var>*prod(<var class="var">d</var>) <var class="var">m</var>]</code>, where the
<code class="code"><var class="var">i</var> + (1:<var class="var">d</var>)</code> rows are the coefficients of all the <var class="var">d</var>
polynomials in the interval <var class="var">i</var>.
</p>
</dd>
<dt><var class="var">n</var></dt>
<dd><p>Number of polynomial pieces or intervals,
<code class="code"><var class="var">n</var> = length (<var class="var">x</var>) - 1</code>.
</p>
</dd>
<dt><var class="var">k</var></dt>
<dd><p>Order of the polynomial plus 1.
</p>
</dd>
<dt><var class="var">d</var></dt>
<dd><p>Number of polynomials defined for each interval.
</p></dd>
</dl>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFmkpp">mkpp</a>, <a class="ref" href="#XREFppval">ppval</a>, <a class="ref" href="One_002ddimensional-Interpolation.html#XREFspline">spline</a>, <a class="ref" href="Signal-Processing.html#XREFpchip">pchip</a>.
</p></dd></dl>
<a class="anchor" id="XREFppval"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-ppval"><span><code class="def-type"><var class="var">yi</var> =</code> <strong class="def-name">ppval</strong> <code class="def-code-arguments">(<var class="var">pp</var>, <var class="var">xi</var>)</code><a class="copiable-link" href="#index-ppval"> ¶</a></span></dt>
<dd><p>Evaluate the piecewise polynomial structure <var class="var">pp</var> at the points <var class="var">xi</var>.
</p>
<p>If <var class="var">pp</var> describes a scalar polynomial function, the result is an array
of the same shape as <var class="var">xi</var>. Otherwise, the size of the result is
<code class="code">[pp.dim, length(<var class="var">xi</var>)]</code> if <var class="var">xi</var> is a vector, or
<code class="code">[pp.dim, size(<var class="var">xi</var>)]</code> if it is a multi-dimensional array.
</p>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFmkpp">mkpp</a>, <a class="ref" href="#XREFunmkpp">unmkpp</a>, <a class="ref" href="One_002ddimensional-Interpolation.html#XREFspline">spline</a>, <a class="ref" href="Signal-Processing.html#XREFpchip">pchip</a>.
</p></dd></dl>
<a class="anchor" id="XREFppder"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-ppder"><span><code class="def-type">ppd =</code> <strong class="def-name">ppder</strong> <code class="def-code-arguments">(pp)</code><a class="copiable-link" href="#index-ppder"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-ppder-1"><span><code class="def-type">ppd =</code> <strong class="def-name">ppder</strong> <code class="def-code-arguments">(pp, m)</code><a class="copiable-link" href="#index-ppder-1"> ¶</a></span></dt>
<dd><p>Compute the piecewise <var class="var">m</var>-th derivative of a piecewise polynomial
struct <var class="var">pp</var>.
</p>
<p>If <var class="var">m</var> is omitted the first derivative is calculated.
</p>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFmkpp">mkpp</a>, <a class="ref" href="#XREFppval">ppval</a>, <a class="ref" href="#XREFppint">ppint</a>.
</p></dd></dl>
<a class="anchor" id="XREFppint"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-ppint"><span><code class="def-type"><var class="var">ppi</var> =</code> <strong class="def-name">ppint</strong> <code class="def-code-arguments">(<var class="var">pp</var>)</code><a class="copiable-link" href="#index-ppint"> ¶</a></span></dt>
<dt class="deftypefnx def-cmd-deftypefn" id="index-ppint-1"><span><code class="def-type"><var class="var">ppi</var> =</code> <strong class="def-name">ppint</strong> <code class="def-code-arguments">(<var class="var">pp</var>, <var class="var">c</var>)</code><a class="copiable-link" href="#index-ppint-1"> ¶</a></span></dt>
<dd><p>Compute the integral of the piecewise polynomial struct <var class="var">pp</var>.
</p>
<p><var class="var">c</var>, if given, is the constant of integration.
</p>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFmkpp">mkpp</a>, <a class="ref" href="#XREFppval">ppval</a>, <a class="ref" href="#XREFppder">ppder</a>.
</p></dd></dl>
<a class="anchor" id="XREFppjumps"></a><span style="display:block; margin-top:-4.5ex;"> </span>
<dl class="first-deftypefn">
<dt class="deftypefn" id="index-ppjumps"><span><code class="def-type"><var class="var">jumps</var> =</code> <strong class="def-name">ppjumps</strong> <code class="def-code-arguments">(<var class="var">pp</var>)</code><a class="copiable-link" href="#index-ppjumps"> ¶</a></span></dt>
<dd><p>Evaluate the boundary jumps of a piecewise polynomial.
</p>
<p>If there are <em class="math">n</em> intervals, and the dimensionality of <var class="var">pp</var> is
<em class="math">d</em>, the resulting array has dimensions <code class="code">[d, n-1]</code>.
</p>
<p><strong class="strong">See also:</strong> <a class="ref" href="#XREFmkpp">mkpp</a>.
</p></dd></dl>
</div>
<hr>
<div class="nav-panel">
<p>
Next: <a href="Miscellaneous-Functions.html">Miscellaneous Functions</a>, Previous: <a href="Derivatives-_002f-Integrals-_002f-Transforms.html">Derivatives / Integrals / Transforms</a>, Up: <a href="Polynomial-Manipulations.html">Polynomial Manipulations</a> [<a href="index.html#SEC_Contents" title="Table of contents" rel="contents">Contents</a>][<a href="Concept-Index.html" title="Index" rel="index">Index</a>]</p>
</div>
</body>
</html>
|
