Thread Subject: solution/fit for f(x)=A*cos(B*x)

Hello,

I'm trying to solve f(x)=A*cos(B*x) to get estimates of the parameters A and B where I have a vector for f(x),x. I'm also looking to give a confidence interval for f(x)=Acos(Bx) given f(x),x and report errors for A and B. Any help would be much appreciated!

Thanks.

Leaf Man <imiller25@live.com> wrote in message
<24427446.1215576592399.JavaMail.jakarta@nitrogen.mathforum.org>...
> Hello,
>
> I'm trying to solve f(x)=A*cos(B*x) to get estimates of
the parameters A and B where I have a vector for f(x),x. I'm
also looking to give a confidence interval for f(x)=Acos(Bx)
given f(x),x and report errors for A and B. Any help would
be much appreciated!
>
> Thanks.

Hi

If nobody answers, I am sending you a dirty solution of a
simulated example:

% acosbx.m
% Model parameters:
a = 10;
b = 3;
x = (0:.125:5)';
y = a*cos(b*x) + randn(size(x)); % cosine with noise
plot(x,y,'or')
% Solution:
ae = max(y); % Initial guess of a
w = acos(y/ae)./x;
I = find(diff(w)<0);
J = find(diff(I)==1);
be = mean(w(I(2:round(end/3)))); % Initial guess of b
fx = @(c) c(1)*cos(c(2)*x) - y;
[c,ssq,cnt] = LMFnlsq(fx,[ae,be]) % from FEX, Id 17534
c = real(c);
hold on
plot(x,fx(c)+y)
hold off

The solution c = [a; b] has been obtained in 5 iterations:

>> acosbx
ae =
   11.0337
be =
    3.0964
c =
    9.7358
    2.9971
ssq =
   45.5168
cnt =
     5

Mira

"Miroslav Balda" <balda.nospam@cdm.it.cas.cz> wrote in
message <g51pbu$5ln$1@fred.mathworks.com>...

Hi

I analyzed the problem in more detail and found a clean
solution in initial guess. Here is the final code:

% acosbx.m
% Model parameters:
a = 10;
b = 3;
x = (0:.125:5)';
y = a*cos(b*x) + randn(size(x)); % cosine with noise
plot(x,y,'or')

% Solution:
L = extr(abs(y)); % from FEX, Id 10272
if L{1}(end)==true, L{1}(end)=false; end
ymax = y(L{1});
ae = mean(abs(ymax)); % mean of all extremes
be = 2*length(xmax)/(xmax(end)-1);% Initial guess of b
fx = @(c) c(1)*cos(c(2)*x) - y;
[c,ssq,cnt] = LMFnlsq(fx,[ae,be]) % from FEX, Id 17534
c = real(c);
hold on
plot(x,fx(c)+y), grid
hold off

The solution has been found as follows:

>> acosbx
c =
    9.5286
    3.0080
ssq =
   42.9551
cnt =
     5
Hope it helps.

Mira

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 parameters =
A and B where I have a vector for f(x),x. I'm also looking to give a confid=
ence interval for f(x)=3DAcos(Bx) given f(x),x and report errors for A and =
B. Any help would be much appreciated!
>
> Thanks.

F =3D fft(f);

etc

Hope this helps.

Greg

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