Path: news.mathworks.com!newsfeed-00.mathworks.com!newsfeed2.dallas1.level3.net!news.level3.com!postnews.google.com!k37g2000hsf.googlegroups.com!not-for-mail
From: Greg Heath <heath@alumni.brown.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: solution/fit for f(x)=A*cos(B*x)
Date: Thu, 10 Jul 2008 01:26:59 -0700 (PDT)
Organization: http://groups.google.com
Lines: 82
Message-ID: <7439521d-9f0f-409b-ab60-f7f662b0d019@k37g2000hsf.googlegroups.com>
References: <24427446.1215576592399.JavaMail.jakarta@nitrogen.mathforum.org> 
NNTP-Posting-Host: 69.141.173.117
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
X-Trace: posting.google.com 1215678419 11023 127.0.0.1 (10 Jul 2008 08:26:59 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Thu, 10 Jul 2008 08:26:59 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: k37g2000hsf.googlegroups.com; posting-host=69.141.173.117; 
User-Agent: G2/1.0
X-HTTP-UserAgent: Mozilla/4.0 (compatible; MSIE 7.0; Windows NT 
Xref: news.mathworks.com comp.soft-sys.matlab:478550



On Jul 10, 12:24=A0am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jul 9, 12:09=A0am, Leaf Man <imille...@live.com> wrote:
>
> > Hello,
>
> > I'm trying to solve f(x)=3DA*cos(B*x) to get estimates of the parameter=
s A and B where I have a vector for f(x),x. I'm also looking to give a conf=
idence interval for f(x)=3DAcos(Bx) given f(x),x and report errors for A an=
d B. Any help would be much appreciated!
>
> > Thanks.
>
> F =3D fft(f);

Example:

close all, clear all, clc, k =3D 0

a0  =3D 10;   b0  =3D 3;   t   =3D (0:0.125:5)';
N  =3D length(t)                         % N =3D 41 (odd),
dt =3D 0.125,   T  =3D 5+dt            % T =3D  5.125

y0          =3D a0*cos(b0*t);
meany0  =3D mean(y0)              % 0.4473
stdy0     =3D std(y0)                   % 7.0850

randn('state',0)
n        =3D 0.1*stdy0*randn(N,1);
meann =3D mean(n)                  % 0.0788
stdn    =3D std(n)                       % 0.6770

y1        =3D y0 + n;                   % cosine with noise
SNRdB =3D 10*log10(var(y0)/var(n))  % 20.3956

k =3D k+1, figure(k),hold on
plot(t,y0),
plot(t,y1,'ro')

meany1 =3D mean(y1)                  % 0.5261
y          =3D y1 - meany1;
y01      =3D y0 - meany0;

k =3D k+1, figure(k), hold on
plot(t,y01),hold on
plot(t,y,'ro'),hold on

% Solution:
df  =3D 1/T                               % 0.1951
dfM =3D b0/(2*pi)/10                 % 0.0477
M   =3D ceil(N*df/dfM)               % 168

f   =3D dfM*(0:M-1)';
Q  =3D ceil((M+1)/2)                 % 85
fQ =3D (Q-1)*dfM                     % 4.0107
fb =3D f-fQ;

Y  =3D dt*fft(y,M);

Yb         =3D fftshift(Y);
ReYb     =3D real(Yb);
ImYb      =3D imag(Yb);
absYb    =3D abs(Yb);
phaseYb =3D angle(Yb);

[maxabsYb indx] =3D max(absYb)       % [25.1987  75]
a =3D 2*maxabsYb*df                 % 9.8336
b =3D 2*pi*abs(fb(indx))             % 3.0000

k=3Dk+1,figure(k)
subplot(2,2,1)
plot(fb,ReYb)
subplot(2,2,2)
plot(fb,ImYb)
subplot(2,2,3)
plot(fb,absYb)
subplot(2,2,4)
plot(fb,phaseYb)


Hope this helps,

Greg