HOW TO make a least-squares fit in a power spectrum.

The situation:
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) (wi)2

Then a power spectrum estimate has the following property

ν P(ω) / Pe is χ2-distributed with ν degrees of freedom,

where Pe is the expected value of power and  ν=N/(T I).
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, Pe is the power that we have to model..

Step 1:
Minimise

χ2=∑ (i=plock) [ √(P(ωi)) - √(Pe(p1, p2,...,ωi)) /σi ]2   ................  (1)

where σi = √(P(ωi)/ν), i.e. find the parameters p1, p2,... of your power spectrum model that yield this
minimum. Our b-model (power-law model of spectral power) is a two-parameter model,

Pe(P0, b, ω) = P0 ω-b  ........................................  (2)

This gives the first solution for b, still biased because of crude statistics.


Step 2:
Form the modified spectrum Q(ωi)
xi = ν P(ωi) / Pe(P0, b, ωi)   ∀ i=plock  ..........................  (3)

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

Q(ωi) =  √(Pe(P0, b, ωi)) g(xi) /μ   ...............................  (5)

and minimise

χ2=∑ (i=plock) [ Q(ωi) - √(Pe(P0, b, ωi)) ]2 ..........................  (6)

Now you have a better guess for b and P0. With those you can redo (3)...(6).  Hopefully the process converges.
The 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 = [x1 .. xn ]
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