www.gusucode.com > matlab 案例源码 matlab代码程序 > matlab/DataInterpolationUsingFFTExample.m

    %% Interpolate Asteroid Data
% In a paper by Gauss, he describes an approach to estimating the orbit of
% the Pallas asteroid.  He starts with the following twelve 2-D positional
% data points |x| and |y|.
x = 0:30:330;
y = [408 89 -66 10 338 807 1238 1511 1583 1462 1183 804];
plot(x,y,'ro')
xlim([0 360])

%%
% Gauss models the asteroid's orbit with a trigonometric polynomial of the
% following form.
%
% $$\begin{array}{rcl}
% y = a_{0} &+& a_{1}\cos(2\pi(x/360))+b_{1}\sin(2\pi(x/360))\\
% &&{}a_{2}\cos(2\pi (2x/360))+b_{2}\sin(2\pi(2x/360))\\
% &&{}\cdots \\
% &&{}a_{5}\cos(2\pi(5x/360))+b_{5}\sin(2\pi(5x/360))\\
% &&{}a_{6}\cos(2\pi(6x/360))
% \end{array}$$

%%
% Use |fft| to compute the coefficients of the polynomial.
m = length(y);
n = floor((m+1)/2);
z = fft(y)/m;

a0 = z(1); 
an = 2*real(z(2:n));
a6 = z(n+1);
bn = -2*imag(z(2:n));

%%
% Plot the interpolating polynomial over the original data points.
hold on
px = 0:0.01:360;
k = 1:length(an);
py = a0 + an*cos(2*pi*k'*px/360) ...
        + bn*sin(2*pi*k'*px/360) ...
        + a6*cos(2*pi*6*px/360); 

plot(px,py)