Question

I have a certain power spectrum array saved as test, depending on a frequency array f.

This power spectrum generally looks like in the following figure enter image description here

It is to be remarked that the above power spectrum stemms from a simulated time series.

What I want to do, in principle, is to fit a curve approaching the simulated power spectrum as shown below:

enter image description here

I know that the theoretical power spectrum can be defined as follows:

function ps_theo = ps_th(L,Uhub,f)

const = L/Uhub;
f_x = 6.*f.*(L/Uhub);
exp = 5/3;
ps_theo = (4*const)./((1 + f_x).^exp);

end

where L is a constant length scale, Uhub, a constant speed and f the frequency vector.

The question is: I don't know the value of 'L', therefore I was thinking of an optimized non-linear resolution by using lsqcurvefit.

I have been proceeding as follows:

xdata = f;
ydata = test;
Uhub = 10;

fit_func = @(L) ps_th(L,Uhub,f);

L_opt = lsqcurvefit(@fit_func,330.2,xdata,ydata)

which retrieves an error msg on the number of inputs variables for the fit_func function.

Would you mind to shed a light?

Was it helpful?

Solution

The function you're fitting can only take in two arguments. You can rewrite ps_th like this:

function ps_theo = ps_th(x0,f)
L = x0(1);
Uhub = x0(2);

const = L/Uhub;
f_x = 6.*f.*(L/Uhub);
exp = 5/3;
ps_theo = (4*const)./((1 + f_x).^exp);

end

Then call lsqcurvefit with something like this:

x0Start = [330.2,10]; % vector of initial parameters
x0_opt = lsqcurvefit(@ps_th,x0Start,xdata,ydata);
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top