# Finding the linear prediction coefficients for a sampled sinusoid

I am unable to work my way through this problem -

A sinusoid with random phase is given as 𝑦(𝑛) = 2 sin(0.25π𝑛 + 𝜃). If this signal is sampled at 10 KHz, estimate the coefficients of a second order LP model. Find the magnitude and phase values of the second order LP error filter at 𝜔 = 0.5𝜋. Estimate the bias (error) in the peak frequency (frequency at the maximum magnitude response) of 𝐻(𝜔)

I am thinking of proceeding with Yule-Walker equations but unable to solve this.

Can you please provide a way to solve this.

• If you rewrite it as a sin + cos, you'll have it in a form where you can use a linear regression model to estimate the coefficients. Bloomfeld explains the details of how to do this at the beginning of his text. In fact, his setting is more complex but reduces to yours for the special case where there is only one frequency. – mark leeds May 14 '20 at 13:52
• Is the phase θ really random, as in a random variable, or is it just unspecified? If it is a random variable, then its probability distribution must be given in order to compute the coefficients R ( i ) = E [ x ( n ) x ( n + i )]. – Joe Mack May 15 '20 at 15:37
• No theta does not hold any special relevance here – ann May 21 '20 at 17:30

In the absence of more details, I assume that

1. the phase $$\Theta$$ is a random variable and
2. $$\Theta$$ is uniformly distributed on the interval $$[0,2\pi]$$, so that its probability density function is $$\frac{1}{2\pi}$$ on that interval.

This is a common choice of model for such problems, and I will show that its guarantees that the stochastic process $$x$$ is wide-sense stationary (WSS).

$$$$\begin{split} \mathsf{E}[x[n]] &=~ \mathsf{E}[2\sin(0.25\pi n + \Theta)]\\ &=~ 2\mathsf{E}[\sin(0.25\pi n + \Theta)]\\ &=~ 2\int_{0}^{2\pi}\sin(0.25\pi n + \theta)\frac{1}{2\pi}d\theta\\ &=~ \frac{1}{\pi}\int_{0}^{2\pi}\sin(0.25\pi n + \theta)d\theta \end{split}$$$$ The integral is the integral of sine over a $$2\pi$$-wide interval, so it is equal to 0. Hence $$\mathsf{E}[x[n]] = 0$$, independent of $$n$$.

$$$$\begin{split} \mathsf{Var}[x[n]] &=~ \mathsf{E}\left[\left(x[n] - \mathsf{E}[x[n]]\right)^2\right]\\ &=~ \mathsf{E}\left[\left(x[n] - 0\right)^2\right]\\ &=~ \mathsf{E}\left[x[n]^2\right]\\ &=~ \mathsf{E}\left[4\sin^2(0.25\pi n + \Theta)\right]\\ &=~ 4\int_{0}^{2\pi}\sin^2(0.25\pi n + \theta)\frac{1}{2\pi}d\theta\\ &=~ \frac{2}{\pi}\int_{0}^{2\pi}\sin^2(0.25\pi n + \theta)d\theta\\ &=~ \frac{2}{\pi}\times\pi ~=~ 2, \end{split}$$$$ which is independent of $$n$$. We see that $$x$$ is a WSS process.

This is important because we rely on expected values of the form $$$$\begin{split} \mathsf{E}[x[n]x[m]] &=~ \mathsf{E}[4\sin(0.25\pi n + \Theta)\sin(0.25\pi m + \Theta)]\\ &=~ \frac{4}{2\pi}\int_{0}^{2\pi}\sin(0.25\pi n + \theta)\sin(0.25\pi m + \theta)d\theta \end{split}$$$$ I leave it as an exercise to show that the value of this integral depends only on the difference $$n-m$$ and not on $$n$$ and $$m$$ separately.

In particular, we can define $$R[i] = \mathsf{E}[x[n]x[n-i]]$$ and know that $$R[i]$$ is really a function of $$i$$ alone and has no dependence on $$n$$.

For second-order linear predictive coding (LPC), we estimate $$x[n]$$ from the previous two values: $$$$\widehat{x}[n] = a_1x[n-1] + a_2x[n-2].$$$$ The error is then $$$$\begin{split} e[n] &=~ x[n] - \widehat{x}[n]\\ &=~ x[n] - a_1x[n-1] - a_2x[n-2]. \end{split}$$$$ The coefficients $$a_1$$ and $$a_2$$ must be chosen to minimize mean square error, which is $$\mathsf{E}[e[n]^2]$$.

To minimize mean square error, the error must be stochastically orthogonal to the random variables that make up the estimate:

$$\begin{eqnarray} \mathsf{E}[e[n]x[n-1]] &=& 0,\\ \mathsf{E}[e[n]x[n-2]] &=& 0. \end{eqnarray}$$ These will give us linear equations whose solution reveals the proper choice of $$a_1$$ and $$a_2$$.

$$$$\begin{split} \mathsf{E}[e[n]x[n-1]] &=~ \mathsf{E}\left[\left(x[n] - a_1x[n-1] - a_2x[n-2]\right)x[n-1]\right)]\\ &=~ \mathsf{E}[x[n]x[n-1]] - a_1\mathsf{E}[x[n-1]^2] - a_2\mathsf{E}[x[n-2]x[n-1]]\\ &=~ R[1] - a_1R[0] - a_2R[1] \end{split}$$$$ Since we want this quantity to be equal to zero, we have the linear equation $$$$R[0]a_1 + R[1]a_2 = R[1].$$$$
$$$$\begin{split} \mathsf{E}[e[n]x[n-2]] &=~ \mathsf{E}\left[\left(x[n] - a_1x[n-1] - a_2x[n-2]\right)x[n-2]\right)]\\ &=~ \mathsf{E}[x[n]x[n-2]] - a_1\mathsf{E}[x[n-1]x[n-2]] - a_2\mathsf{E}[x[n-2]^2]\\ &=~ R[2] - a_1R[1] - a_2R[0] \end{split}$$$$ Since we also want this quantity to be zero, we have the linear equation $$$$R[1]a_1 + R[0]a_2 = R[2].$$$$
We now have the system of equations $$$$\left(\begin{array}{cc} R[0] & R[1]\\ R[1] & R[0] \end{array}\right) \left(\begin{array}{c}a_1\\a_2\end{array}\right) = \left(\begin{array}{c}R[1]\\R[2]\end{array}\right),$$$$ whose solution is $$$$\left(\begin{array}{c}a_1\\a_2\end{array}\right) = \frac{1}{R[0]^2 - R[1]^2} \left(\begin{array}{c} R[0] & -R[1]\\-R[1] & R[0] \end{array}\right) \left(\begin{array}{c}R[1]\\R[2]\end{array}\right)$$$$
I leave it as an exercise to compute the correlations $$R[0]$$, $$R[1]$$, and $$R[2]$$.

An aside: In higher-order LPC we take advantange of the very nice structure of the matrix in the linear equations (it's a Toeplitz matrix). That structure allows one to use the so-called Levinson-Durbin algorithm to solve it numerically. While the matrix is real symmetric and thus is numerically amenable to numerous algorithms, Levinson-Durbin is the one used in practice because of its speed and because it yields other useful quantities as it runs. These reflection coefficients (RCs) that it yields along the way, are proxies for the linear prediction coefficients (LPCs).
I should note here that there are sometimes differences in sign. See, for example, the help page for MATLAB's lpc function. But once the signs are chosen, there is consistency afterward. Check the definition of LPC in your example before starting.
The filter to compute the error is $$x[n] - a_1x[n-1] - a_2x[n-2]$$, so the impulse response is $$\mathbf{h} = (h_0, h_1, h_2) = (1, -a_1, -a_2)$$. The corresponding frequency response is $$$$H(e^{i\omega}) = 1 - a_1e^{-i\omega} - a_2e^{-2i\omega}.$$$$
In some applications (such as vocoders), the error ($$x[n] - \widehat{x}[n]$$) is assumed to be a noise process. If the LPCs have been received, then a realization of the noise process is generated, and an approximation of the original signal is created by using the inverse filter, which is an IIR filter with $$z$$-transform $$$$\frac{1}{H(z)} = \frac{1}{1 - a_1z^{-1} - \cdots - a_kz^{-k}}.$$$$