How to estimate spectral index of a time series??? Please help me
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
How to estimate spectral index of a time series??? Please help me
3 Kommentare
Akzeptierte Antwort
Wayne King
am 10 Jan. 2013
Bearbeitet: Wayne King
am 10 Jan. 2013
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
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 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Spectral Estimation finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!