Spectral Estimation
Arash
Komaee, Afshin Sepehri
Department
of Electrical and Computer Engineering
{akomaee, afshin}@eng.umd.edu
1.
Introduction
Auto-correlation function is an important statistics describing a stochastic process. This function describes the time-correlation between samples of a stochastic process. The corresponding frequency-domain representation of auto-correlation is power spectral density. By definition, power spectral density of a stochastic process is the Fourier transform of its auto-correlation function.
Spectral estimation in general is the procedure of estimating the power spectral density of a random signal based on finite number of observations. The accuracy of the estimated spectrum is strictly related to the number of observations, the more observation we have, we can estimate the spectrum more accurate. Another important consideration regarding to the spectral estimation is that most often the observations of the signal are corrupted by noise.
There are two main categories of spectral estimation methods, non-parametric or classical methods, and parametric methods. Historically, non-parametric method has been invented before parametric methods. This method and its modifications are based on direct computation of Fourier transform of observation vector. Non-parametric methods specially, for small number of observation does not present satisfactory performance. In contrast to non-parametric methods, parametric methods avoid of direct computation of Fourier transform, instead they try to fit a suitable model to the provided data. These methods generally have better characteristics.
In this project, we examine the estimation of the power spectral density of four stochastic processes, based on observation vectors of length 64. During the project we try both classical and parametric methods, with emphasis on parametric methods.
2.
Spectral Estimation – An Overview
As it was mentioned, there are two different approaches for spectral estimation. The first one, which includes classical power spectral estimation methods, is based on the periodogram. Second group of methods are model-based or parametric methods.
In both these the approaches, we try to find power spectral density of a signal from finite-duration observations. This finite length of data has negative impact on the quality of the power spectrum estimation. In continue, we will quickly discuss the method.
2.1. Non-Parametric
Methods- Periodogram
Provided N samples 
of a random
process, we can obtain the periodogram of the signal using
(2.1)
As it can be seen, it is directly related to the Fourier transform of
the signal. It can be shown that the variance of 
does not decay to zero
as 
. Hence we can conclude that periodogram is not a consistent
estimate of the true power density spectrum. Another problem is the effect of
limited number of observation that results a low-resolution and distorted estimate
for the spectrum of the signal. To overcome the mentioned problems, several
modifications have been proposed. We will discuss them in this section.
2.1.1.
To solve the problem of convergence for periodogram,
Therefore, the average periodogram can be obtained as
(2.2)
where
(2.3)
If we calculate
the variance of 
, it decays to zero as . That is, the variance of
2.1.2.
Welch Method – Averaging Modified Priodograms
In Welch method, we make two modifications to the
where D
is the size of the segments and 
. According to this segmentation, we make the second
modification, which is to window data segments prior to computing the
periodogram. As a result
(2.4)
where
(2.5)
and the Welch spectrum estimate is the average of these modified periodograms,
(2.6)
2.1.3.
Windowing
Using a limited number of observations in spectrum analysis and estimation is equivalent to multiplying the signal in time-domain with a rectangular window. Limited number of samples is then processed to resolve closely spaced sinusoids present in the signal. In fact the Fourier transform of the window is convolved with the signal in frequency domain. One approach to improve the results is to use a window other that rectangular, because this window does not satisfy some frequency requirements well. Instead, according to the application and based on the spectrum, some other specific window shapes can give better answers.
There are some popular windows such as Hamming, Hann,
Figure 2.1 – Four windows in time
domain – Left to right and top to the bottom: Rectangular, Hamming, Hanning,
Figure 2.2 – Four windows in frequency
domain – Left to right and top to the bottom: Rectangular, Hamming, Hanning,
There are two important parameters in window design: Window Length and Window Shape. These two parameters determine the performance measures of the main lobe width and the side lobes’ power.
The window length determines the main lobe width. In general, as a window’s length is increased, the main lope width decreases. A narrow main lobe width provides good resolution. That is, the ability of the filter to separate very close frequencies in the signal.
The window shape is in fact the type of the window chosen from the above list. The effect of the window shape is on the leakage from the other frequencies to the main lobe. This causes problem in detecting the amplitude of the sinusoids.
As it can be concluded from the impacts of the two parameters in window, we should select an appropriate window based on our application. For the sake of this project, we apply Rectangular and Hamming windows to the signals and estimate spectrums as non-parametric results and compare them with parametric modeling result.
2.2. Parametric
Methods
The non-parametric methods discussed about in previous subsection are simple to understand and to apply. However they suffer some disadvantages.
First, to achieve reasonable resolution for many application, we need a long set of data. Second, due to windowing, these methods suffer leakage. This leakage may restrict us in finding weak signals, which exist in the data.
The basic limitation of these methods is that they assume auto-correlation coefficients outside the window zero which is an unrealistic assumption. As a solution to these weaknesses, we try to fit a model to our signal so it can be extrapolated easily to the future times. This approach decreases the number of the parameters and also increases estimation accuracy.
In this section we look at the four common models of random processes, Auto-Regressive (AR), Moving Average (MA), Auto-Regressive Moving-Average (ARMA), and sinusoidal models. In fact the first two models are particular cases of the general case ARMA.
2.2.1.
Auto-Regressive Models
An AR process model is defined as
(2.7)
where 
is a realization of
the process and v(n) is a white-noise process with variance s2.
Given the AR parameters 
, the power spectrum density (PSD) of 
is found as
(2.8)
Auto-regressive model is the most widely used model of the three linear models introduced earlier, because first it is suitable for representing spectra with narrow peaks and second a number of linear equations need to be solved for finding the AR parameters.
· Yule-Walker Method: In this method, we simply estimate the autocorrelation coefficients from the data and use
(2.9)
The auto-correlation is estimated using
(2.10)
· Burg Method
The Burg method is a least-squares optimization problem with the constraint that the reflection coefficients obtained from the problem must satisfy Levinson-Durbin recursion, i.e.
(2.11)
subject to constrain
(2.12)
which results
(2.13)
The Burg method has three advantages.
(1) It has high frequency resolution.
(2) It results stable AR model.
(3) Because it uses Levinson-Durbin recursion.
2.2.2.
Moving Average Models
A MA process can be defined as
(2.14)
Given the MA
parameters 
, the power spectrum density (PSD) of is found as
(2.15)
To find , we are dealing with some non-linear equations, so we prefer
to convert it to AR form and then solve it. Since most of the times, the MA
model for spectral estimation, is a high-order process, we prefer not to use
it.
2.2.3.
ARMA Models
An ARMA process is a combination of AR and MA models. In this model we have both zeros and poles.
(2.16)
Given the MA
parameters , the power spectrum density (PSD) of is found as
(2.17)
As it was shown in AR section, the Burg algorithm provides a reliable, high-resolution spectrum estimates based on the AR model. An ARMA model can decrease the order of the AR solution, especially when the data has been corrupted by noise. We use this property in estimating the models in next section.
2.2.4.
Sinusoidal Model
Another generally used model for stochastic processes is a linear
combination of sinusoidal signals with constant amplitudes and frequencies and
a phase randomly distributed in 
. This model can be written as
(2.18)
It can be shown that the auto-correlation function of such random process is
(2.19)
where 
is defined as
(2.20)
The power spectral density of the process as Fourier transform of the above function is obtained as
(2.21)
For this process we can see that
(2.22)
We will use the above equations for spectral estimation of signals in this project.
2.2.5.
Selection of Model Order
Selection of the model order of the AR process is an important issue. If we do not choose the order high enough, we may come with a too smooth spectrum. On the other side, too low order may cause low-level peaks in the spectrum.
There are many different researches and criteria on choosing an appropriate order for the model. We mention here just three of them:
· The first criterion called the final prediction error (FPE) criterion tries to minimize the performance index
(2.23)
where 
is the estimated variance of the linear prediction error.
· The second criterion called Akaike information criterion (AIC), is selecting the order by minimizing
(2.24)
· The third criterion is to minimize the description length (MDL) where MDL is defined as
(2.25)
During this project we will use AIC to obtain the AR order, but to be familiar with two other criteria, we make the plot of them.
3.
Spectral Estimation of Provided Signals
We will use same approach to estimate the spectral density of all four records of observations. Although or main focus is not using non-parametric methods, for each signal we first obtain the periodogram of the signal based on rectangular window, Hamming window, and Welch method. These periodograms will serve as a basis of comparison in next steps. They also provide an initial guess about the shape and order of parametric model.
In the next step we found the order of an AR model, which best fits the data. To perform this task we find the AIC of the signal based on minimizing (2.24). After finding the required order, we can find an appropriate AR model based on Yule-Walker or Burg methods. Although we obtain the AR model based on both methods, we observed that Burg method could perform the task in more reliable way, so for next steps we have used the results given by Burg method.
To represent the model in more compact form, for each spectrum obtained by Burg method, we tried to fit a lower order MA or ARMA to the AR model. To perform this task, one can use the MATLAB function “invfreqz”. In all cases we have obtained the order of ARMA model such that it fits the AR model with small enough error.
At last step we fit the results of Burg estimation to the sinusoidal
model. Note that each sinusoid represents a pick in spectrum of signal, so in
equation (2.19), M is equal to the number of local maximums. Also the
frequencies 
in (2.19) are
equivalent to the frequencies that the spectrum has local maxima, thus a
computer program searching in the spectrum can find the value of . In the other hand Lacoss (1971) showed that spectral
picks obtained from AR models are proportional to square of corresponding
sinusoidal power. This fact is used to obtain . Also from (2.20) one can easily obtain 
of (2.18). Finally to
make sure that the obtained frequencies of sinusoidal model are valid, we use
the MUSIC algorithm as a basis of comparison. In using MUSIC algorithm we set
the order of signal space equal to 2M.
The results of the above procedure (with the mentioned order) are presented in the reminder of this section. For each signal numbered by K=1, 2, 3, 4, Figure 3.K.1 shows the periodogram of given observations for rectangular window, Hamming window, and Welch method. Figure 3.K.2 illustrate the functions log(FPE(p)), AIC(p), and MDL(p)/N. Note that to obtain the order of AR model we only use the minimum of AIC. In Figure 3.K.3 we have illustrated the spectrum of signal estimated by AR model using Yule-Walker and Burg methods. Finally as a comparison between results obtained by parametric and non-parametric methods, the spectrum estimated based on Burg method and Hamming windowed periodogram is given in Figure 3.K.4.
For each signal, Table 3.K.1, Table 3.K.2, and Table 3.K.3 show the AR parameters, ARMA parameters and sinusoidal model parameters respectively.
3.1. First
Signal
3.1.1.
Non-Parametric Method
Figure 3.1.1 – Periodogram spectral
estimations of first signal
Blue: Rectangular window,
Green: Hamming window, Red: Welch Method
3.1.2.
Parametric Method
Figure 3.1.2 – Selecting the optimum model
order
Blue: log(FPE), Green: AIC,
Red: MDL, Optimal Order based on AIC is 21
Figure 3.1.3 – Spectral estimation on Blue:
Burg Method, Green: Yule-Walker Method
3.1.3.
Comparison of parametric and non-Parametric
Methods
Figure 3.1.4 – Comparison between
parametric (Blue: Burg Method) and non-parametric (Green: Periodogram with
hamming window)
3.1.4.
AR, ARMA, and Sinusoidal Models Parameters
K
|
Coefficients
|
K |
Coefficients |
K |
Coefficients |
K |
Coefficients
|
|
0 |
1.0000 |
7 |
0.2226 |
14 |
-0.1679 |
21 |
0.1924 |
|
1 |
-0.2582 |
8 |
0.3267 |
15 |
-0.2140 |
22 |
|
|
2 |
-0.2284 |
9 |
0.2029 |
16 |
0.2845 |
23 |
|
|
3 |
-0.3674 |
10 |
-0.1957 |
17 |
-0.0455 |
24 |
|
|
4 |
0.2318 |
11 |
-0.2125 |
18 |
-0.1315 |
25 |
|
|
5 |
-0.0011 |
12 |
0.4062 |
19 |
0.0050 |
26 |
|
|
6 |
-0.4841 |
13 |
-0.1764 |
20 |
0.2332 |
27 |
|
|
Variance: 0.6472 |
|||||||
Table 3.1.1 - AR mode parameters based on
Burg method
K |
Denominator |
Numerator |
|
0 |
1 |
1.0157 |
|
1 |
-1.3541 |
-1.1003 |
|
2 |
0.592 |
0.5364 |
|
3 |
-0.2217 |
0.0259 |
|
4 |
-0.3224 |
-0.8026 |
|
5 |
0.3594 |
0.6073 |
|
6 |
-0.4031 |
-0.086 |
|
7 |
0.7663 |
-0.1777 |
|
8 |
-0.1072 |
0.4066 |
|
9 |
-0.0846 |
-0.1132 |
|
10 |
0.0688 |
-0.1038 |
|
11 |
-0.5202 |
|
|
12 |
0.5474 |
|
|
13 |
-0.6312 |
|
|
14 |
0.2564 |
|
|
15 |
0.1275 |
|
|
16 |
0.0549 |
|
Variance |
0.6472 |
|
|
p |
16 |
|
|
q |
10 |
|
Table 3.1.2 - ARMA model parameters based
on AR modeling of original spectrum and model reduction procedure
|
k |
Frequency |
Power |
a Parameter |
|
1 |
0.1296 |
0.7856 |
1.2534 |
|
2 |
0.3858 |
0.8819 |
1.3280 |
|
3 |
0.7647 |
0.0619 |
0.3518 |
|
4 |
1.0515 |
0.1364 |
0.5223 |
|
5 |
1.2901 |
0.1032 |
0.4543 |
|
6 |
1.8408 |
0.1202 |
0.4903 |
|
7 |
2.2442 |
0.1971 |
0.6279 |
|
8 |
2.9268 |
0.1389 |
0.5270 |
Table 3.1.3 - Sinusoidal model parameters
3.2. Second
Signal
3.2.1.
Non-Parametric Method
Figure 3.2.1 – Periodogram spectral
estimations of second signal
Blue: Rectangular window,
Green: Hamming window, Red: Welch Method
3.2.2.
Parametric Method
Figure 3.2.2 – Selecting the optimum model
order
Blue: log(FPE), Green: AIC,
Red: MDL, Optimal Order based on AIC is 9
Figure 3.2.3 – Spectral estimation on Blue:
Burg Method, Green: Yule-Walker Method
3.2.3.
Comparison of parametric and non-Parametric
Methods
Figure 3.2.4 – Comparison between
parametric (Blue: Burg Method) and non-parametric (Green: Periodogram with
hamming window)
3.2.4. AR, ARMA, and Sinusoidal Models Parameters
K
|
Coefficients
|
K |
Coefficients |
K |
Coefficients |
K |
Coefficients
|
|
0 |
1.0000 |
3 |
1.3081 |
6 |
0.1665 |
9 |
-0.2744 |
|
1 |
1.0310 |
4 |
2.0480 |
7 |
1.0392 |
10 |
|
|
2 |
0.1583 |
5 |
0.4443 |
8 |
0.3141 |
11 |
|
|
Variance: 1.3764 |
|||||||
Table 3.2.1 - AR mode parameters based on
Burg method
K |
Denominator |
Numerator |
|
0 |
1.0000 |
0.9990 |
|
1 |
1.3879 |
0.3530 |
|
2 |
0.6539 |
0.1235 |
|
3 |
1.5416 |
0.0432 |
|
4 |
2.5983 |
0.0123 |
|
5 |
1.3717 |
|
|
6 |
0.6563 |
|
|
7 |
1.2736 |
|
|
8 |
0.7687 |
|
Variance |
1.3764 |
|
|
p |
8 |
|
|
q |
4 |
|
Table 3.2.2 - ARMA model parameters based
on AR modeling of original spectrum and model reduction procedure
|
K |
Frequency |
Power |
a Parameter |
|
1 |
0.9296 |
13.1029 |
5.1192 |
|
2 |
1.1022 |
6.4994 |
3.6054 |
|
3 |
2.5364 |
6.7342 |
3.6699 |
|
4 |
2.9729 |
6.3818 |
3.5726 |
Table 3.2.3 - Sinusoidal model parameters
3.3. Third
Signal
3.3.1.
Non-Parametric Method
Figure 3.3.1 – Periodogram spectral
estimations of third signal
Blue: Rectangular window,
Green: Hamming window, Red: Welch Method
3.3.2.
Parametric Method
Figure 3.3.2 – Selecting the optimum model
order
Blue: log(FPE), Green: AIC,
Red: MDL, Optimal Order based on AIC is 26
Figure 3.3.3 – Spectral estimation on Blue:
Burg Method, Green: Yule-Walker Method
3.3.3.
Comparison of parametric and non-Parametric
Methods
Figure 3.3.4 – Comparison between
parametric (Blue: Burg Method) and non-parametric (Green: Periodogram with
hamming window)
3.3.4. AR, ARMA, and Sinusoidal Models Parameters
K
|
Coefficients
|
K |
Coefficients |
K |
Coefficients |
K |
Coefficients
|
|
0 |
1.0000 |
7 |
-0.8543 |
14 |
0.7216 |
21 |
-0.0625 |
|
1 |
0.3216 |
8 |
-0.0313 |
15 |
0.1555 |
22 |
-0.2151 |
|
2 |
-0.0222 |
9 |
0.0247 |
16 |
-0.6540 |
23 |
0.2375 |
|
3 |
-0.0020 |
10 |
0.7417 |
17 |
-0.5113 |
24 |
0.0819 |
|
4 |
0.9818 |
11 |
-0.2257 |
18 |
-0.2112 |
25 |
0.1936 |
|
5 |
-0.0143 |
12 |
-0.4981 |
19 |
0.5479 |
26 |
-0.3345 |
|
6 |
0.0447 |
13 |
-0.0661 |
20 |
-0.1764 |
27 |
|
|
Variance: 0.1796 |
|||||||
Table 3.3.1 - AR mode parameters based on
Burg method
K |
Denominator |
Numerator |
|
0 |
1.0000 |
1.0211 |
|
1 |
0.8713 |
0.5345 |
|
2 |
0.1713 |
0.0517 |
|
3 |
0.5163 |
0.6250 |
|
4 |
1.9282 |
0.8666 |
|
5 |
0.3918 |
-0.4204 |
|
6 |
-0.2891 |
-0.2270 |
|
7 |
-0.1479 |
0.2568 |
|
8 |
0.1119 |
-0.2038 |
|
9 |
-1.0405 |
-0.7181 |
|
10 |
0.0326 |
-0.0125 |
|
11 |
-0.2600 |
0.0128 |
|
12 |
-0.6624 |
-0.0701 |
|
13 |
-0.2065 |
0.2456 |
|
14 |
1.4986 |
0.5398 |
|
15 |
0.3208 |
0.2345 |
|
16 |
-0.1620 |
0.3111 |
|
17 |
0.4721 |
0.4012 |
|
18 |
0.8014 |
0.2301 |
Variance |
0.1796 |
|
|
p |
18 |
|
|
q |
18 |
|
Table 3.3.2 - ARMA model parameters based
on AR modeling of original spectrum and model reduction procedure
|
K |
Frequency |
Power |
a Parameter |
|
1 |
0.4042 |
0.0171 |
0.1851 |
|
2 |
0.9135 |
0.5320 |
1.0315 |
|
3 |
1.0078 |
0.7197 |
1.1998 |
|
4 |
1.5723 |
3.9877 |
2.8241 |
|
5 |
1.9788 |
0.0169 |
0.1838 |
|
6 |
2.1913 |
3.6671 |
2.7082 |
|
7 |
2.5019 |
1.1344 |
1.5062 |
|
8 |
2.6139 |
0.8730 |
1.3213 |
|
9 |
2.7313 |
1.6899 |
1.8384 |
Table 3.3.3 - Sinusoidal model parameters
3.4. Forth
Signal
3.4.1.
Non-Parametric Method
Figure 3.4.1 – Periodogram spectral
estimations of forth signal
Blue: Rectangular window,
Green: Hamming window, Red: Welch Method
3.4.2.
Parametric Method
Figure 3.4.2 – Selecting the optimum model
order
Blue: log(FPE), Green: AIC,
Red: MDL, Optimal Order based on AIC is 11
Figure 3.4.3 – Spectral estimation on Blue:
Burg Method, Green: Yule-Walker Method
3.4.3.
Comparison of parametric and non-Parametric
Methods
Figure 3.4.4– Comparison between parametric
(Blue: Burg Method) and non-parametric (Green: Periodogram with hamming window)
3.4.4. AR, ARMA, and Sinusoidal Models Parameters
|
K |
Coefficients
|
K |
Coefficients |
K |
Coefficients |
K |
Coefficients
|
|
0 |
1.0000 |
3 |
1.0708 |
6 |
1.3192 |
9 |
0.5051 |
|
1 |
0.4382 |
4 |
1.3401 |
7 |
0.7995 |
10 |
0.2796 |
|
2 |
1.3885 |
5 |
1.0528 |
8 |
0.8694 |
11 |
0.2020 |
|
Variance: 0.8304 |
|||||||
Table 3.4.1 - AR mode parameters based on
Burg method
K |
Denominator |
Numerator |
|
0 |
1.0000 |
0.9985 |
|
1 |
-0.2354 |
-0.6746 |
|
2 |
1.5310 |
0.4319 |
|
3 |
0.0494 |
-0.2802 |
|
4 |
1.2853 |
0.1710 |
|
5 |
0.1911 |
-0.1251 |
|
6 |
1.1880 |
0.0785 |
|
7 |
0.0033 |
-0.0629 |
|
8 |
0.8442 |
|
|
9 |
-0.0551 |
|
|
10 |
0.2991 |
|
Variance |
0.8304 |
|
|
p |
10 |
|
|
q |
7 |
|
Table 3.4.2 - ARMA model parameters based
on AR modeling of original spectrum and model reduction procedure
|
K |
Frequency |
Power |
a Parameter |
|
1 |
0.7386 |
0.6699 |
1.1575 |
|
2 |
1.4872 |
2.0776 |
2.0384 |
|
3 |
2.5357 |
0.6858 |
1.1711 |
Table 3.4.3 - Sinusoidal model parameters
4.
Summery and Conclusion
In this project we performed spectral estimation on four different records of observations. We used the Burg method as the center of our work. For each signal after obtaining AR model, we extract a more efficient model representing the spectrum. In this way, for each signal we built ARMA and sinusoidal model based on AR model. A summery of results is gathered in Table 4.
|
Model |
AR |
ARMA |
Sinusoidal |
|
|
Order |
Numerator Order |
Denominator Order |
Number of Sinusoids |
|
|
21 |
10 |
16 |
8 |
|
|
9 |
4 |
8 |
4 |
|
|
26 |
18 |
18 |
9 |
|
|
11 |
7 |
10 |
3 |
|
Table 4.1 – Summary of the results
Also the frequency components obtained by sinusoidal model is illustrated in Figure 4.1. In this figure, the lines are located at sinusoid frequencies and their magnitude is proportional to corresponding sinusoid’s power.
Figure 4.1 - Frequency components of
signals.
(up left signal 1, up right
signal 2, bottom left signal 3, and bottom right signal4)
In this project we observed that
- In non-parametric methods using a rectangular window results a distorted spectrum comparing to Hamming window.
- The frequency resolution of rectangular window is better than Hamming window.
- The Welch method results a very smooth spectrum, ignoring the very near adjacent frequency components.
- In comparison between Yule-Walker and Burg methods, the second one is superior in frequency resolution. As it is clear from Figure 3.K.3, Yule-Walker method cannot distinguish the nearby frequencies.
- Using ARMA model instead of AR can reduce the complexity of model in some cases. In contrast for some other cases it cannot reduce the complexity.
- The signals provided in this project can be well approximated using sinusoidal model.
- MUSIC method validates the results obtained by Burg method.
References
[1] Proakis John G., Rader Charles M., Ling
Fuyun, Nikias Chrysostomos L., “Advanced Digital Signal Processing”, Macmillan,
1992
[2] Talha Ahsan S., “Optimal Window Design, comparing it with the Kaiser Windw”, IEEE potentials, December 2002, pp. 40-44
[3] MATLAB Signal Processing Toolbox: http://www.mathworks.com/access/helpdesk/help/toolbox/signal/signal.shtml