Hello 

The output of our frequency estimation methods gives a frequency
function f(n), where 'n' is the discrete-time index.

We know that this frequency function f(n) can be nonlinear. 
That is f(n) = ax^3 + bx^2 + cx + d.

I am looking for line/curve fitting or any other method/algorithm to
estimate the coefficients a,b,c & d given f(n).  The method is
required to be used with a Digital Signal Processor(DSP).

Please comment.
Also suggest some references/newsgroups/links/ideas etc.


Awaiting, eagerly for your reply.


With regards
Parthasarathy

On Thu, 08 Jul 2004 21:07:40 -0700, Parthasarathy wrote:

> Hello 
> 
> The output of our frequency estimation methods gives a frequency
> function f(n), where 'n' is the discrete-time index.
> 
> We know that this frequency function f(n) can be nonlinear. 
> That is f(n) = ax^3 + bx^2 + cx + d.
> 
> I am looking for line/curve fitting or any other method/algorithm to
> estimate the coefficients a,b,c & d given f(n).  The method is
> required to be used with a Digital Signal Processor(DSP).
> 
> Please comment.
> Also suggest some references/newsgroups/links/ideas etc.
> 
> 
> Awaiting, eagerly for your reply.
> 
> 
> With regards
> Parthasarathy
 

Look for Numerical Recipes in C and download it. It describes several
curve fitting algorithms.

D

Parthasarathy wrote:

> The output of our frequency estimation methods gives a frequency
> function f(n), where 'n' is the discrete-time index.
> 
> We know that this frequency function f(n) can be nonlinear. 
> That is f(n) = ax^3 + bx^2 + cx + d.
> 
> I am looking for line/curve fitting or any other method/algorithm to
> estimate the coefficients a,b,c & d given f(n).  The method is
> required to be used with a Digital Signal Processor(DSP).
> 
> Please comment.
> Also suggest some references/newsgroups/links/ideas etc.

Take a look at
The C++ Scalar, Vector, Matrix and Tensor class library

	http://www.netwood.net/~edwin/svmtl/

cat svmtl/examples/polynomial.cc
/*
The C++ Scalar, Vector, Matrix and Tensor classes.
Copyright (C) 1998 E. Robert Tisdale

This file is part of The C++ Scalar, Vector, Matrix and Tensor classes.
This library is free software which you can redistribute and/or modify
under the terms of the GNU Library General Public License
as published by the Free Software Foundation;
either version 2, or (at your option) any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty
of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.

You should have received a copy of
the GNU Library General Public License along with this library.
If not, write to the Free Software Foundation,
Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

Written by E. Robert Tisdale
*/
/*

This little program fits a polynomial

         y(x) = w_0 + w_1*x + w_2*x^2 + ... + w_n*x^n

of order n to a set of m > n data pairs {x_i, y_i}
from standard input and prints the coefficients to standard output.

         $ make polynomial
         $ polynomial 2 < data
         w =
         -3.44668 1.01307 0.499563
         $ gnuplot
         gnuplot> y(x) = -3.44668 + 1.01307*x + 0.499563*x*x
         gnuplot> plot "data", y(x)

*/

#include<doubleSquare.h>

inline
doubleVector
         solve(const doubleSubVector& b, doubleSubSquare& A) {
   // solve b = xA^T assuming that A is a positive definite symmetric matrix
   offsetVector  p = A.lld();    // Cholesky decomposition PAP^T = 
(LD)(LD)^T
   return b.pld(p, A).dup(A.t(), p);
   // triangular system solvers b = y(P^T(LD))^T and y = x((DU)P)^T
  }

int
main(int argc, char* argv[]) {

   // The algorithm:

   // For an overdetermined set of equations
   //
   //    Y = wX
   //
   // where
   //
   //    Y = [y_1, y_2, ..., y_i, ..., y_m],
   //
   //    X = [(u_1)^T, (u_2)^T, ..., (u_i)^T, ..., (u_m)^T]
   //
   // and
   //
   //    u_i = [1, x_i, (x_i)^2, ..., (x_i)^j, ..., (x_i)^n],
   //
   // find the set of coefficients
   //
   //    w = [w_0, w_1, ..., w_j, ..., w_n]
   //
   // which represents a least squares best fit
   // of a polynomial of order n to the data pairs {x_i, y_i}.
   // First, post multiply both sides of the equation by X^T
   //
   //    v = YX^T = wXX^T = wS
   //
   // then solve v = wS for w.
   //
   // The square symmetric matrix
   //
   //    S = S_1 + S_2 + ... + S_i + ... + S_m
   //
   // where the outer products
   //
   //    S_i = ((u_i)^T)u_i
   //
   // and
   //
   //    v = y_1*u_1 + y_2*u_2 + ... + y_i*u_i + ... + y_m*u_m.

   // The implementation:

   using std::cin; using std::cout;

   Extent n = (1 < argc)? atoi(argv[1]): 1;      // polynomial order

   doubleSquare          S(1+n, 0.0);            // symmetric matrix
   doubleVector          u(1+n, 1.0);
   doubleVector          v(1+n, 0.0);
   double                x;
   double                y;
   while (cin >> x && cin >> y) {                // {x_i, y_i}
     for (Offset j = 0; j < n; ++j)
       u[j+1] = x*u[j];                          // u_ij = (x_i)^j
     S += u.submatrix().t().dot();               // S += ((u_i)^T)u_i
     v += y*u;                                   // v += y_i*u_i
     }

   doubleVector          w = solve(v, S);        // solve v = wS for w
   cout << "w =\n" << w;

   return 0;
   }