/ Home
S-Archive Download Fortran - Download Script
pronyfreq Frequency Estimation
Estimates a sum of sinusoidal signals using an eigenanalysis method.
pronyfreq(y, nfreq=1, constant=T, tol=1e-5, maxit=40, trace=F, warnings=T)
y numeric vector.
nfreq number of frequencies to estimate.
constant logical variable, indicating whether to include a constant in the model.
tol tolerance indicating required accuracy.
maxit maximum number of iterations.
trace logical variable. If true then the top half of the Prony vector and the smallest and largest eigenvalues of the Prony matrix are output at each iteration.
warnings logical variable. If false, then run quietly without issuing warnings about non-convergence. Useful when using pronyfreq for a few iterations to get starting values for another estimation technique.
A list with the following components:
prony numeric vector of coefficients of the Prony difference equation.
freq numeric vector of frequencies.
coef numeric vector of coefficients. The first component is the constant term, then the coefficients of cos(f*t) for each frequency, then the coefficients of sin(f*t) for each frequency.
fitted numeric vector. The value of the sum of sinusoidal frequencies.
residuals numeric vector, equal to y - fitted
The signal is assumed to be of the form
yt = m + a1 cos(w1 t) + b1 sin(w1 t) + ... +  ak cos(wk t) + bk sin(wk t)
for t = 0 ... length(y)-1 and k = nfreq. If  constant=F  then  m  is not included. The frequences  wj   and the coefficients are estimated by an algorithm described in Smyth (2000). The estimators are not exactly least squares, but are equivalent or better than least squares in mean square error. On output, the component prony holds the coefficients of the polynomial with roots 1, and exp( i wj) for each frequency. The component freq holds (w1,...,wk). The component coef holds (m,a1,...,ak,b1,...,bk).

It is possible for the algorithm to return fewer frequencies than requested. In this case the data sequence is not well fitted by a sum of sinusoids. The algorithm works extremely well on data sets of moderate size. The algorithm used differences of y, and therefore will be numerically unstable if the length of y is large and the observations are sampled very densely relative to the frequencies.

The algorithm implemented by pronyfreq is the Osborne/Breslow/Macovski algorithm from Smyth (2000). Pisarenko's covariance method is used to start the iteration.

The function pronyfreq uses fortran subroutines for solving banded linear systems.
Smyth, G. K. (2000). Employing symmetry constaints for improved frequency estimation by eigenanalysis methods. Technometrics 42, 277-289.
Consider the variable star data set.
> out <- pronyfreq(star,2)
> out$freq
[1] 0.2166596 0.2617983

The estimated frequencies are 0.2167 and 0.2618.
> 2*pi/out$freq
[1] 29.00027 24.00010

These correspond to periods of almost exactly 24 and 29 days.
> out$coef
[1] 17.0857822411 7.6477245734 0.0009274355 6.4906186247 7.0845674915

The constant term is 17.086.
> plot(star)
> lines(0:599,out$fitted,col=2)
pronyfr1.gif (6851 bytes)


S-Archive Download Fortran - Download Script

Gordon Smyth. Copyright © 1996-2016. Last modified: 10 February 2004