Thread Subject: Trigonometric polynomial approximation

Hi,

I've a specific problem at hand and wondered whether MATLAB provides already this functionality or if I have to implement it by myself.

I have a (discrete) time series over some periods of length L and want to approximate this series with trigonometric polynomials via the least squares method. Thus, the function is of the form a_0 + \sum_{k = 1}^n a_k * sin(2* \pi * k/L) + b_k * cos(2* \pi * k/L) with the time series at hand given over some periods, e.g. N * L, and the coefficients a_0, a_k, b_k to be determined.

I've already tried a workaround with the Gaussian fft Asteroid interpolation example, but it gets mixed up since my time series is given over several periods and I want only a best fitting trigonometric polynom, not an interpolation.

Does MATLAB provide the desired functionality (where I can choose the degree n freely) or does anyone know how to easiest implement it?

Thanks for your help.

Thomas

> Hi,
>
> I've a specific problem at hand and wondered whether
> MATLAB provides already this functionality or if I
> have to implement it by myself.
>
> I have a (discrete) time series over some periods of
> length L and want to approximate this series with
> trigonometric polynomials via the least squares
> method. Thus, the function is of the form a_0 +
> \sum_{k = 1}^n a_k * sin(2* \pi * k/L) + b_k * cos(2*
> \pi * k/L) with the time series at hand given over
> some periods, e.g. N * L, and the coefficients a_0,
> a_k, b_k to be determined.
>

I think, the function is of the form
a_0 +
\sum_{k = 1}^n a_k * sin(2* \pi * k* t/L) + b_k * cos(2*
\pi * k * t/L)

> I've already tried a workaround with the Gaussian fft
> Asteroid interpolation example, but it gets mixed up
> since my time series is given over several periods
> and I want only a best fitting trigonometric polynom,
> not an interpolation.
>
> Does MATLAB provide the desired functionality (where
> I can choose the degree n freely) or does anyone know
> how to easiest implement it?
>
> Thanks for your help.
>
> Thomas

This is a linear least squares problem.
Formulate the problem as A*x = b with x the vector
of unknown coefficients and b the vector of your
measurement. Then use x = A\b to solve.

Best wishes
Torsten.

> I think, the function is of the form
> a_0 +
> \sum_{k = 1}^n a_k * sin(2* \pi * k* t/L) + b_k * cos(2*
> \pi * k * t/L)

Yes, I forgot t, thanks.

> This is a linear least squares problem.
> Formulate the problem as A*x = b with x the vector
> of unknown coefficients and b the vector of your
> measurement. Then use x = A\b to solve.

Thus, do I have to set up Matrix A manually (or by an algorithm) every time I change n? Is there no function already implemented in MATLAB providing above functionality? I thought this is a very common used function.

"Thomas " <meinl@iism.uni-karlsruhe.de> wrote in message <gk76ec$ett$1@fred.mathworks.com>...
> > I think, the function is of the form
> > a_0 +
> > \sum_{k = 1}^n a_k * sin(2* \pi * k* t/L) + b_k * cos(2*
> > \pi * k * t/L)
>
> Yes, I forgot t, thanks.
>
> > This is a linear least squares problem.
> > Formulate the problem as A*x = b with x the vector
> > of unknown coefficients and b the vector of your
> > measurement. Then use x = A\b to solve.
>
> Thus, do I have to set up Matrix A manually (or by an algorithm) every time I change n? Is there no function already implemented in MATLAB providing above functionality? I thought this is a very common used function.

No, it is apparently not terribly common,
although Fourier series are used often enough.
And this is just that - a Fourier series, built
from scattered data. The coefficients you want
are solved for using least squares, so it is not
an interpolant.

I'll argue that the reason why people don't use
this more often is that better, more useful tools
exist - splines. Least squares splines are easily
built and easily enough controlled to have a
desired shape.

But nothing stops you from using the ability of
matlab to write a function to solve this problem.
This is what makes a language like matlab so
valuable - that we can essentially customize it to
our own purposes. If I frequently do something
that the rest of the world finds less useful than
I do, I can write a function to solve it. Give the
function an easily remembered name and
document it so you can use it whenever you
need it.

If the code is good enough, then put it on
the file exchange, so that maybe someone
else can augment their own version of matlab
too.

John

> > I think, the function is of the form
> > a_0 +
> > \sum_{k = 1}^n a_k * sin(2* \pi * k* t/L) + b_k *
> cos(2*
> > \pi * k * t/L)
>
> Yes, I forgot t, thanks.
>
> > This is a linear least squares problem.
> > Formulate the problem as A*x = b with x the vector
> > of unknown coefficients and b the vector of your
> > measurement. Then use x = A\b to solve.
>
> Thus, do I have to set up Matrix A manually (or by an
> algorithm) every time I change n? Is there no
> function already implemented in MATLAB providing
> above functionality? I thought this is a very common
> used function.

I think you will have to do it manually, but that's
not that hard:

A(i,1) = 1 (i=1,...,n)
A(i,2*j) = sin(2*pi*t(i)/L*j)
A(i,2*j+1) = cos(2*pi*t(i)/L*j) (i,j = 1,...,n)
b(i) = f(t(i)) (i=1,...,n)

Ordering of solution vector:
x = (a_0,a_1,b_1,a_2,b_2,...,a_n,b_n)

Best wishes
Torsten.

Thanks a lot. Meanwhile, I did it myself like this:

function x = fourierintp(Data,n,period)

k = 1:n;

A = zeros(length(Data),1+2*n);
for j = 1:length(Data)
    A(j,:) = [1 sin(2*pi*k*j/period) cos(2*pi*k*j/period)];
end

x = A\Data;

fourierintpfun = zeros(1,length(Data));
for j = 1:length(Data)
    fourierintpfun(j) = x(1) + sum(sin(2*pi*k*j/period)*x(2:2+n-1)) + sum(cos(2*pi*k*j/period)*x(2+n:end));
end

I know it's not very nice because of the loops, but I don't work that long with MATLAB. Any comments on how to improve the code are appreciated.

Thomas

> Thanks a lot. Meanwhile, I did it myself like this:
>
> function x = fourierintp(Data,n,period)
>
> k = 1:n;
>
> A = zeros(length(Data),1+2*n);
> for j = 1:length(Data)
> A(j,:) = [1 sin(2*pi*k*j/period)
> cos(2*pi*k*j/period)];

A(j,:) = [1 sin(2*pi*k*t(j)/period)cos(2*pi*k*t(j)/period)];

> end
>
> x = A\Data;
>
> fourierintpfun = zeros(1,length(Data));

fourierintpfun = A*x ;
The rest is superfluos.

> for j = 1:length(Data)
> fourierintpfun(j) = x(1) +
> 1) + sum(sin(2*pi*k*j/period)*x(2:2+n-1)) +
> sum(cos(2*pi*k*j/period)*x(2+n:end));
> end
>
> I know it's not very nice because of the loops, but I
> don't work that long with MATLAB. Any comments on how
> to improve the code are appreciated.
>
> Thomas

Best wishes
Torsten.

Just another question:

What (or is there any) difference between using

x = A\Data;

and

x = lsqr(A,Data);

?

> Just another question:
>
> What (or is there any) difference between using
>
> x = A\Data;
>
> and
>
> x = lsqr(A,Data);
>
> ?

Different numerical methods are used, but the results
should be the same.
With x = A\data, x is computed by a direct method for
the linear system and should be used if A is not
too large.
With x = lsqr(A,Data), x is computed iteratively
and should be used if A is large and sparse.

Best wishes
Torsten.

On Jan 9, 3:45=A0am, "Thomas " <me...@iism.uni-karlsruhe.de> wrote:
> Hi,
>
> I've a specific problem at hand and wondered whether MATLAB provides alre=
ady this functionality or if I have to implement it by myself.
>
> I have a (discrete) time series over some periods of length L and want to=
 approximate this series with trigonometric polynomials via the least squar=
es method. Thus, the function is of the form a_0 + \sum_{k =3D 1}^n a_k * s=
in(2* \pi * k/L) + b_k * cos(2* \pi * k/L) with the time series at hand giv=
en over some periods, e.g. N * L, and the coefficients a_0, a_k, b_k to be =
determined.
>
> I've already tried a workaround with the GaussianfftAsteroid interpolatio=
n example, but it gets mixed up since my time series is given over several =
periods and I want only a best fitting trigonometric polynom, not an interp=
olation.
>
> Does MATLAB provide the desired functionality (where I can choose the deg=
ree n freely) or does anyone know how to easiest implement it?
>
> Thanks for your help.
>
> Thomas

Search the srchives with

greg-heath DFT

Hope this helps.

Greg

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <8478126.1231510981344.JavaMail.jakarta@nitrogen.mathforum.org>...
> > Just another question:
> >
> > What (or is there any) difference between using
> >
> > x = A\Data;
> >
> > and
> >
> > x = lsqr(A,Data);
> >
> > ?
>
> Different numerical methods are used, but the results
> should be the same.
> With x = A\data, x is computed by a direct method for
> the linear system and should be used if A is not
> too large.
> With x = lsqr(A,Data), x is computed iteratively
> and should be used if A is large and sparse.
>
> Best wishes
> Torsten.

As A is neither large and sparse, I used the first method. But now I get warnings like

Warning: Rank deficient, rank = 405, tol = 1.1631e-011. and
Warning: Rank deficient, rank = 11, tol = 1.1631e-011


Usage of

x = pinv(A)*data;

provided the same solution as

x = A\data.

I do not really know how to interpret the above warnings, and what methods (the first one or the pseudoinverse) should be used.

Thanks for your help.

Thomas

> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> in message
> <8478126.1231510981344.JavaMail.jakarta@nitrogen.mathf
> orum.org>...
> > > Just another question:
> > >
> > > What (or is there any) difference between using
> > >
> > > x = A\Data;
> > >
> > > and
> > >
> > > x = lsqr(A,Data);
> > >
> > > ?
> >
> > Different numerical methods are used, but the
> results
> > should be the same.
> > With x = A\data, x is computed by a direct method
> for
> > the linear system and should be used if A is not
> > too large.
> > With x = lsqr(A,Data), x is computed iteratively
> > and should be used if A is large and sparse.
> >
> > Best wishes
> > Torsten.
>
> As A is neither large and sparse, I used the first
> method. But now I get warnings like
>
> Warning: Rank deficient, rank = 405, tol =
> 1.1631e-011. and
> Warning: Rank deficient, rank = 11, tol =
> 1.1631e-011
>
>
> Usage of
>
> x = pinv(A)*data;
>
> provided the same solution as
>
> x = A\data.
>
> I do not really know how to interpret the above
> warnings, and what methods (the first one or the
> pseudoinverse) should be used.
>
> Thanks for your help.
>
> Thomas

Your problem reads
A*x = b
with A having (2*n+1) colums and say m >> 2n+1 rows.
The warning message means that the rank of the matrix
is less than (2*n+1) so that the solution to the
associated least-squares problem to determine x is
no longer unique.
x = A\data provides gives you an arbitrary solution while
x=pinv(A)*data provides the x of minimum norm.

Best wishes
Torsten.