Your power spectrum is based on a windowed autocovariance (it always is. If you're unaware: implicitly). Assume

you have sampled an N-point time series, the window has nonzero coeffcients between -T and +T with an RMS-value of I

I
= ∑_{(-T ≤
i ≤ T)} (w_{i})^{2}

Then a power spectrum estimate has the following property

ν P(ω)
/ P_{e} is
χ^{2}-distributed with ν degrees of freedom,

where P

Thus a full-length rectangular window yields ν=1 (since autocovariance goes from -N and +N and T=N then).

That's the case in a standard periodogram.

Now, P

Step 1:

Minimise

χ^{2}=∑ _{(i=plock)} [
√(P(ω_{i})) - √(P_{e}(p_{1}, p_{2},...,ω_{i})) /σ_{i} ]^{2} ................
(1)

where σminimum. Our b-model (power-law model of spectral power) is a two-parameter model,

P_{e}(P_{0}, b, ω) = P_{0 }ω^{-b}
........................................ (2)

x_{i}
= ν P(ω_{i})
/
P_{e}(P_{0},
b, ω_{i}) ∀ i=plock ..........................
(3)

μ = √( chisqi(ν,
½) )
......................................... (4)

Q(ω_{i}) = √(P_{e}(P_{0},
b, ω_{i})) g(x_{i})
/μ
............................... (5)

and minimise

χ^{2}=∑ _{(i=plock)} [ Q(ω_{i}) - √(P_{e}(P_{0},
b, ω_{i})) ]^{2} ..........................
(6)

Now you have a better guess for b
and PThe function g that provides the mapping into a normally distributed random variable is

g(x, ν,
μ) = pnormali(
chisqu(
ν, x ), μ )
........................... (7)

and

function
y=chisq(nu,x);

y=gammainc(x/2,nu/2);

end

function y=chisqi(nu,p);

f=@(x)chisq(nu,x)-p;

y=fzero(f,[0 4*nu]);

end

function y=pnormal(x,mu)

y=0.5+((x>mu)-0.5).*gammainc((x-mu).^2/2,0.5);

end

function y=pnormali(p,mu);

f=@(x)pnormal(x,mu)-p;

y=fzero(f,[0 20]);

end

Because of the @(x)
handle construct, you cannot call chisqi and pnormi for a
vector x = [xy=gammainc(x/2,nu/2);

end

function y=chisqi(nu,p);

f=@(x)chisq(nu,x)-p;

y=fzero(f,[0 4*nu]);

end

function y=pnormal(x,mu)

y=0.5+((x>mu)-0.5).*gammainc((x-mu).^2/2,0.5);

end

function y=pnormali(p,mu);

f=@(x)pnormal(x,mu)-p;

y=fzero(f,[0 20]);

end

Since I am a nobody in matlab I'd appreciate help to code this up such that you can avoid an outer for-loop.

Finally, the plock scheme is to be designed such that successive spectral samples are separated by a full mainlobe of the window spectrum, in order to avoid correlation. And remember that the spectrum sampling starts a zero frequency, and the power there is not well defined. So plock should start picking beyond the fourth or so spectral bin.

hgs didit again