Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

How to estimate spectral index of a time series??? Please help me

Asked by Chris Martin on 4 Jan 2013

How to estimate spectral index of a time series??? Please help me

3 Comments

Rick Rosson on 7 Jan 2013

Please explain in more detail.

Greg Heath on 7 Jan 2013

In particular, what is your definition of a spectral index?

Chris Martin on 7 Jan 2013

The power spectra, P, of many geophysical phenomena are well approximatebdy a power-lawd ependencoen frequencfy of the form [Agnew, 1992] P(f)=Po*f*exp(-a) where 'a' is the spectral index and Po is constant. My intend is to see the type of noise such as white noise,flicker noise present in my time series by estimating the spectral index as white noise has spectral index of 0, flicker noise has spectral index of 1 and random walk has a spectral index of 2. I will be gratefull if you help me out in this occasion

Chris Martin

Products

No products are associated with this question.

1 Answer

Answer by Wayne King on 10 Jan 2013
Edited by Wayne King on 10 Jan 2013
Accepted answer

I'm assuming your equation is incorrect or you are using nonstandard notation. The power law process should be the frequency raised to a power. So something like

 P(f) = Pof^{-a}

where the frequency is raised to the power -a. Or equivalently, 1/f is raised to the power a.

What you can do is to estimate the power spectrum, I would use pwelch() to get a smooth estimate of the power spectrum and then fit a least-squares model to the power spectrum based on a log-log model

Given the above equation, take the log of both sides

log(P(f)) = log(Po)-alog(f)

Now you see this is a linear model in the parameter a

I'll use the publicly available M-file

http://people.sc.fsu.edu/~jburkardt/m_src/cnoise/f_alpha_gaussian.m

for simulating a power law process with a = 2

   x = f_alpha_gaussian (1000, 1, 2 );
   % Now get the Welch estimate
   [pxx,f] = pwelch(x,300,250,length(x),1);
   % examine a loglog plot
   loglog(f,pxx)

You see above that in a loglog plot, the power spectrum is approximately a line with a negative slope, that slope will be your (-a) parameter. Obviously, this will be sensitive to the input parameters you specify for pwelch() so you may have to do some work there.

Now to estimate the slope you'll want to avoid DC (0 frequency) because the log of that will go to -infinity

    %design matrix for linear model
    X = ones(length(pxx(2:end)),2);
    fvar = log(f(2:end));
    X(:,2) = fvar;
    y = log(pxx(2:end));
    % get the least squares estimate
    X\y

For my particular realization of x, I get

   -2.3739
   -1.8764

-2.3739 is the estimate of log(Po) and -1.8764 is the estimate of -a, which in my case I entered to be -2, so -1.8764 is not bad.

2 Comments

Chris Martin on 11 Jan 2013

Sir could you suggest more details about the least square estimate. I am new to matlab

Chris Martin on 3 Mar 2013

How to get the standard error from this estimation?? please help me

Wayne King

Contact us