Zobrazení zdroje stránky Trigonometrická interpolace

'''Trigonometrická interpolace''' je v [[Matematika|matematice]] interpolační metoda zvlášť vhodná pro [[Interpolace|interpolaci]] [[Periodická funkce|periodických funkcí]]. Cílem interpolace je nalezení funkce, která prochází zadanými [[Datový bod|datovými body]]. U trigonometrické interpolace se funkce interpoluje [[Trigonometrický polynom|trigonometrickými polynomy]], což je suma [[Goniometrická funkce|sinů a kosinů]] daných period.

Důležitým speciálním případem je, když dané datové body jsou stejně vzdálené; pak je řešení dáno [[Fourierova transformace#Diskrétní Fourierova transformace|diskrétní Fourierovou transformací]].

== Formulace problému interpolace ==
Trigonometrický polynom stupně ''K'' má tvar
{{Vzorec|<math> p(x) = a_0 + \sum_{k=1}^K a_k \cos(kx) + \sum_{k=1}^K b_k \sin(kx). \,  </math>|1}}
Tento výraz obsahuje 2''K'' + 1 koeficientů, ''a''<sub>0</sub>, ''a''<sub>1</sub>, … ''a''<sub>''K''</sub>, ''b''<sub>1</sub>, …, ''b''<sub>''K''</sub>, které je třeba vypočítat tak, aby funkce procházela ''N'' body:
:<math> p(x_n) = y_n, \quad n=0, \ldots, N-1. \, </math>
Protože trigonometrický polynom je periodický s periodou 2π, lze těchto ''N'' bodů rozdělit a uspořádat do jedné periody:
:<math> 0 \leq x_0 < x_1 < x_2 < \ldots  < x_{N-1} < 2  \pi. \, </math>
Přičemž obecně ''není'' požadováno, aby tyto body byly stejně vzdálené. Úkolem interpolace je pak nalézt takové koeficienty, aby trigonometrický polynom ''p'' vyhovoval interpolačním podmínkám.

== Formulace v komplexní rovině ==
Nejpřirozenější je formulace trigonometrické interpolace v [[Komplexní rovina|komplexní rovině]]. Vzorec pro trigonometrický polynom můžeme přepsat na tvar
<math> p(x) = \sum_{k=-K}^K c_k e^{ikx}, \, </math>
kde ''i'' je [[Imaginární jednotka]]. Pokud položíme ''z'' = ''e''<sup>''ix''</sup>, dostaneme tvar
:<math> q(z) = \sum_{k=-K}^K c_k z^{k}, \, </math>
kde
:<math> q(e^{ix}) \triangleq p(x). \, </math>
Tím je problém trigonometrické interpolace převeden na polynomiální interpolaci na [[Jednotková kružnice|jednotkové kružnici]]. Existence a jednoznačnost trigonometrické interpolace nyní okamžitě vyplývá z odpovídajících výsledků pro polynomiální interpolaci.

Více informací o formulaci trigonometrických interpolačních polynomů v komplexní rovině obsahuje text „Interpolation using Fourier Polynomials“ od strany 156.{{Sfn|Interpolation using Fourier Polynomials}}

== Řešení problému ==
Za výše uvedených podmínek existuje řešení problému pro ''jakoukoli'' danou množinu datových bodů {''x''<sub>''k''</sub>, ''y''<sub>''k''</sub>}, pokud počet datových bodů ''N'' není větší než počet koeficientů v polynomu, tj. ''N''&nbsp;≤&nbsp;2''K''+1 (i pro ''N''>2''K''+1 může existovat řešení, podle toho, jaká je konkrétní množina datových bodů). Interpolační polynom je jednoznačný právě tehdy, když počet nastavitelných koeficientů je roven počtu datových bodů, tj. ''N''&nbsp;=&nbsp;2''K''&nbsp;+&nbsp;1. Ve zbytku tohoto článku budeme předpokládat, že je tato podmínka splněna.

=== Lichý počet bodů ===
Pokud je počet bodů ''N'' lichý (''N=2K+1''), použití [[Lagrangeova interpolace|Lagrangeva vzorce pro polynomiální interpolaci]] na polynomiální formulaci v komplexní rovině dostaneme řešení ve tvaru
{{Vzorec|<math> p(x) = \sum_{k=0}^{2K} y_k\,t_k(x),</math>|5}}
kde
:<math> t_k(x) = e^{-iKx+iKx_k} \prod_{m=0,m\ne k}^{2K} \frac{e^{ix}-e^{ix_m}}{e^{ix_k}-e^{ix_m}}.</math>
Faktor <math>e^{-iKx+iKx_k}</math> je v tomto vzorci kvůli tomu, že formulace pro komplexní rovinu obsahuje i záporné mocniny <math>e^{ix}</math> a proto není polynomiálním výrazem pro <math>e^{ix}</math>. Korektnost tohoto výrazu lze snadno ověřit díky pozorování, že <math>t_k(x_k)=1</math> a že <math>t_k(x)</math> je lineární kombinací vhodných mocnin <math>e^{ix}</math>.
Využitím identity
{{Vzorec|<math>e^{iz_1}-e^{iz_2}=2i\sin\left(\frac{z_1-z_2}{2}\right)e^{i\frac12 z_1+i\frac12 z_2},</math>|2}}
lze koeficient <math>t_k(x)</math> zapsat ve tvaru
{{Vzorec|<math> t_k(x) = \prod_{m=0,m\ne k}^{2K} \frac{\sin\frac12(x-x_m)}{\sin\frac12(x_k-x_m)}.</math>|4}}

=== Sudý počet bodů ===
Pokud počet bodů ''N'' je sudý (''N=2K''), použití [[Lagrangeova interpolace|Lagrangeova vzorce pro polynomiální interpolaci]] na polynomiální formulaci v komplexní rovině dostaneme řešení ve tvaru
{{Vzorec|<math> p(x) = \sum_{k=0}^{2K-1} y_k\,t_k(x),</math>|6}}
kde
{{Vzorec|<math> t_k(x) = e^{-iKx+iKx_k} \frac{e^{ix}-e^{i\alpha_k}}{e^{ix_k}-e^{i\alpha_k}} \prod_{m=0,m\ne k}^{2K-1} \frac{e^{ix}-e^{ix_m}}{e^{ix_k}-e^{ix_m}}.</math>|3}}
Přičemž konstanty <math>\alpha_k</math> lze libovolně zvolit. To je způsobené faktem, že interpolační funkce ({{Odkaz na vzorec|1}}) obsahuje lichý počet neznámých konstant. Běžnou volbou je vyžadovat, že nejvyšší frekvence měla tvar konstanta krát <math>\cos(Kx)</math>, tj. aby člen <math>\sin(Kx)</math> měl nulový koeficient, ale fázi nejvyšší frekvence lze obecně zvolit tak, aby byla <math>\varphi_K</math>. Pro výpočet <math>\alpha_k</math> podle ({{Odkaz na vzorec|2}}) je třeba ({{Odkaz na vzorec|3}}) zapsat ve tvaru
:<math> t_k(x) = \frac{\cos\left(Kx-\frac12\alpha_k+\frac12\sum\limits_{m=0,m\ne k}^{2K-1} x_m\right)+\sum\limits_{m=-(K-1)}^{K-1}c_k e^{imx}}{2^N\sin(\frac{x_k-\alpha_k}{2})\prod\limits_{m=0,m\ne k}^{2K-1}\sin(\frac{x_k-x_m}{2})}.</math>
Což dává
:<math>\alpha_k=\sum_{m=0,m\ne k}^{2K-1} x_m - 2 \varphi_K</math>
a
:<math> t_k(x) = \frac{\sin\frac12(x-\alpha_k)}{\sin\frac12(x_k-\alpha_k)}\prod_{m=0,m\ne k}^{2K-1} \frac{\sin\frac12(x-x_m)}{\sin\frac12(x_k-x_m)}.</math>

Přitom je třeba se vyhnout nekonečnům způsobeným nulovými hodnotami ve jmenovateli.

== Ekvidistantní uzly ==
Další zjednodušení problému je možné,{{Sfn|Zygmund|1988}} pokud uzly <math>x_m</math> jsou ekvidistantní, tj.
:<math>x_m=\frac{2\pi m}{N}.</math>

=== Lichý počet bodů ===
Přestože je možné další zjednodušení pomocí ({{Odkaz na vzorec|4}}), mnohem jednodušší přístup je uvažovat [[Dirichletovo jádro]]
:<math>D(x,N)=\frac{1}{N} +\frac{2}{N} \sum_{k=1}^{\frac12 (N-1)}\cos(kx) = \frac{\sin\frac12 Nx}{N\sin\frac12 x},</math>
kde <math>N>0</math> je liché. Je zřejmé, že <math>D(x,N)</math> je lineární kombinací vhodných mocnin <math>e^{ix}</math> a platí
:<math>D(x_m,N)=\begin{cases}0\text{ pro } m\neq0 \\1\text{ pro } m=0\end{cases}.</math>
Protože tyto dvě vlastnosti definují koeficienty <math>t_k(x)</math> v ({{Odkaz na vzorec|5}}) jednoznačně, z toho plyne, že
:<math>\begin{align}
t_k(x) &= D(x-x_k,N)=\begin{cases}
\frac{\sin\frac12 N(x-x_k)}{N\sin\frac12 (x-x_k)} \text{ pro } x\neq x_k\\
\lim\limits_{x\to 0} \frac{\sin\frac12 Nx}{N\sin\frac12 x}=1 \text{ pro } x= x_k
\end{cases}\\&= \frac{\mathrm{sinc}\,\frac12 N(x-x_k)}{\mathrm{sinc}\,\frac12 (x-x_k)}.
\end{align}</math>
Funkce [[Sinc]] definovaná vztahem
:<math> \mathrm{sinc}\,x=\frac{\sin x}{x}.</math>
zabrání singularitě.

=== Sudý počet bodů ===
Pro <math>N</math> sudé definujeme Dirichletovo jádro jako
:<math>D(x,N)=\frac{1}{N} +\frac{1}{N}\cos \frac12 Nx + \frac{2}{N} \sum_{k=1}^{\frac12 N-1}\cos(kx) = \frac{\sin\frac12 Nx}{N\operatorname{tg}\frac12 x}.</math>
Opět je zřejmé, že <math>D(x,N)</math> je lineární kombinací vhodných mocnin <math>e^{ix}</math>, neobsahuje člen <math> \sin \frac12 Nx </math> a vyhovuje
:<math>D(x_m,N)=\begin{cases}0\text{ pro } m\neq0 \\1\text{ pro } m=0\end{cases}.</math>
Z těchto vlastností vyplývá, že koeficienty <math>t_k(x)</math> v ({{Odkaz na vzorec|6}}) jsou dány vzorcem
:<math>\begin{align}
t_k(x) &= D(x-x_k,N)=\begin{cases}
\frac{\sin\frac12 N(x-x_k)}{N\operatorname{tg}\frac12 (x-x_k)}\text{ pro } x\neq x_k\\
\lim\limits_{x\to 0} \frac{\sin\frac12 Nx}{N\operatorname{tg}\frac12 x}=1 \text{ pro } x= x_k.
\end{cases}\\&= \frac{\mathrm{sinc}\,\frac12 N(x-x_k)}{ \mathrm{sinc}\,\frac12 (x-x_k)}\cos\frac12 (x-x_k)
\end{align}</math>
Přitom <math>t_k(x)</math> neobsahuje <math> \sin \frac12 Nx </math>. Funkce <math> \sin \frac12 Nx </math> bude mít nulovou hodnotu také pro všechny body <math>x_m</math>, proto se násobky tohoto členu obvykle vynechávají.

=== Implementace ===
Na webu University of Delaware byla zveřejněna tato implementace trigonometrické interpolace pro [[MATLAB]]:{{Sfn|triginterp}}
<syntaxhighlight lang="matlab">
function P = triginterp(xi,x,y)
% TRIGINTERP Trigonometrická interpolace.
% Vstup:
%   xi  vektor bodů, v nichž se vyhodnocuje funkce
%   x   ekvidistantní interpolační body (vektor délky N)
%   y   interpolační hodnoty (vektor délky N)
% Výstup:
%   P   vektor hodnot trigonometrického interpolantu
N = length(x);
% Nastavení odstup nezávislé proměnné.
h = 2/N;
scale = (x(2)-x(1)) / h;
x = x/scale;  xi = xi/scale;
% Vyhodnocení interpolantu.
P = zeros(size(xi));
for k = 1:N
  P = P + y(k)*trigcardinal(xi-x(k),N);
end

function tau = trigcardinal(x,N)
ws = warning('off','MATLAB:divideByZero');
% Výpočet je různý pro sudá a lichá N.
if rem(N,2)==1   % liché
  tau = sin(N*pi*x/2) ./ (N*sin(pi*x/2));
else             % sudé
  tau = sin(N*pi*x/2) ./ (N*tan(pi*x/2));
end
warning(ws)
tau(x==0) = 1;     % oprava hodnotu pro x=0
</syntaxhighlight>

=== Spojitost s diskrétní Fourierovou transformací ===
Zvláště důležitý je případ, kdy body ''x''<sub>''n''</sub> jsou stejně vzdálené. V tomto případě platí
:<math> x_n = 2 \pi \frac{n}{N}, \qquad 0 \leq n < N.</math>

Transformace převádějící datové body ''y''<sub>''n''</sub> na koeficienty ''a''<sub>''k''</sub>, ''b''<sub>''k''</sub> je získaná z [[Fourierova transformace#Diskrétní Fourierova transformace|diskrétní Fourierovy transformace]] (DFT) řádu&nbsp;N.

:<math> Y_k = \sum_{n=0}^{N-1} y_n \ e^{-i 2 \pi \frac{nk}{N}} \, </math>

:<math> y_n = p(x_n) = \frac{1}{N} \sum_{k=0}^{N-1} Y_k \ e^{i 2 \pi \frac{nk}{N}} \, </math>

(Kvůli tomu, jak byl výše problém formulován, je třeba se omezit na lichý počet bodů. To není nezbytně nutné, ale pro sudý počet bodů je třeba zahrnout další kosinový člen odpovídající [[Nyquistova frekvence|Nyquistově frekvenci]].)

Případ interpolace používající pouze kosiny pro stejně vzdálené body odpovídající trigonometrické interpolaci se [[Parita funkce|sudou symetrie]] bodů zkoumal [[Alexis Clairaut]] v roce 1754. Řešení je v tomto případě ekvivalentní s [[Diskrétní kosinová transformace|diskrétní kosinovou transformací]]. Rozvoj pro stejně vzdálené body obsahující pouze siny odpovídající liché symetrii vyřešil v roce 1762 [[Joseph-Louis Lagrange]]; řešením je v tomto případě [[diskrétní sinová transformace]]. Úplný interpolační polynom s kosiny i siny, který odpovídá DFT, vyřešil okolo roku 1805 [[Carl Friedrich Gauss]] v nepublikované práci, v níž také odvodil algoritmus [[Cooley–Tukey FFT algoritmus|rychlé Fourierovy transformace]]. Clairaut, Lagrange a Gauss pracovali na výpočtech parametrů [[Oběžná dráha|oběžných drah]] [[Planeta|planet]], [[asteroid]]ů, atd. z konečného počtu pozorování; protože oběžné dráhy jsou periodické, byla trigonometrická interpolace přirozenou volbou.{{Sfn|Heideman|Johnson|Burrus|1984}}

== Použití v numerických výpočtech ==
[[Chebfun]] je plně integrovaný softwarový systém pro výpočty s funkcemi napsaný v [[MATLAB]]u používající trigonometrickou interpolaci a Fourierovy rozvoje pro výpočty s periodickými funkcemi. V [[Chebfun]] je okamžitě dostupné množství algoritmů týkajících se trigonometrické interpolace; několik příkladů je dostupných v popisu výpočtů s periodickými funkcemi.{{Sfn|Chebfun11|loc=Section 11}}

== Odkazy ==

=== Reference ===
{{Překlad|en|Trigonometric interpolation|1072278093}}<references />

=== Literatura ===
* {{Citace monografie
 | jméno = Kendall E.
 | příjmení = Atkinson
 | titul = An Introduction to Numerical Analysis
 | url = https://archive.org/details/introductiontonu0000atki
 | vydání = 2
 | rok = 1988
 | vydavatel = John Wiley & Sons
 | místo = New York
 | isbn = 0-471-50023-2
 | ref = harv
}}
* {{Citace periodika
 | jméno = M. T.
 | příjmení = Heideman
 | jméno2 = D. H.
 | příjmení2 = Johnson
 | jméno3 = C. S.
 | příjmení3 = Burrus
 | url = https://web.archive.org/web/20170914220621/https://pdfs.semanticscholar.org/1790/fe007bc1ab161a1ea814748b42e3acbdc958.pdf
 | titul =  Gauss and the history of the fast Fourier transform
 | periodikum = IEEE ASSP Magazine
 | ročník = 1
 | číslo = 4
 | strany = 14–21
 | rok = 1984
 | ref = harv
}}
* {{Citace periodika
 | jméno = G.B.
 | příjmení = Wright
 | jméno2 = M.
 | příjmení2 = Javed
 | jméno3 = H.
 | příjmení3 = Montanelli
 | jméno4 = L.N.
 | příjmení4 = Trefethen
 | url = http://eprints.maths.ox.ac.uk/1876/1/NA_15_02.pdf
 | titul = Extension of Chebfun to periodic functions
 | periodikum = SIAM. J. Sci. Comput.
 | ričník = 37
 | rok = 2015
 | strany = C554-C573
 | ref = harv
 }} {{Wayback|url=http://eprints.maths.ox.ac.uk/1876/1/NA_15_02.pdf |date=20150915211700 }}
* {{Citace monografie
 | příjmení = Zygmund
 | jméno = A.
 | titul = Trigonometric Series
 | svazek = II
 | kapitola = X
 | vydavatel = Cambridge University Press
 | rok = 1988
 | ref = harv
}}
* {{Citace elektronické monografie
 | url = http://www.physics.arizona.edu/~restrepo/475A/Notes/sourcea.pdf
 | titul = Interpolation using Fourier Polynomials
 | vydavatel = University of Arizona
 | ref = Interpolation using Fourier Polynomials
}}
* {{Citace elektronické monografie
 | url = http://www.chebfun.org/docs/guide/guide11.html
 | titul = Periodic Chebfuns
 | jméno = Grady B.
 | příjmení = Wright
 | jméno2 = Lloyd N.
 | příjmení2 = Trefethen
 | datum = prosinec 2014, červenec 2019
 | rok = 2014
 | ref = Chebfun11
}}
* {{Citace elektronické monografie
 | url = http://www.math.udel.edu/~braun/M428/Matlab/interpolation/triginterp.m
 | ref = triginterp
}}

=== Externí odkazy ===
* [http://www.chebfun.org www.chebfun.org]
{{Autoritní data}}

[[Kategorie:Interpolace]]
[[Kategorie:Goniometrie]]