CR bound analysis of circular ARMA parameters pdf_arma_circ.m.

The circular ARMA model lends itself better to analysis. The CR bound for the circular ARMA model parameters

   $\displaystyle \mbox{\boldmath$\theta$}$$\displaystyle = [\sigma^2, a_1, \ldots a_P, b_1, \ldots b_Q]^\prime
$

are derived now. We analyze the PDF (9.18) for the CR bound. We apply the results of Section 16.2.2 with

$\displaystyle \rho_k = \sigma^2 \; \frac{ \vert B_k\vert^2}{\vert A_k\vert^2}.$ (9.23)

noting that the ARMA parameters are spectral parameter $ \theta$. Following these results, we need the derivatives of the spectral values $ \rho_k$ with respect to the parameters. The first derivatives are (from eq. 9.22 with $ 2L$ replaced by $ N$).

\begin{displaymath}\begin{array}{rcl} \left( {\displaystyle \frac{\partial \rho_...
...t) & =& \frac{\vert B_k\vert^2}{\vert A_k\vert^2 }. \end{array}\end{displaymath} (9.24)

Thus,

\begin{displaymath}\begin{array}{rcl} D(b_m) &=& - \sum_{k=0}^{N-1} \; T_1(k) {\...
... T_1(k) \frac{\vert B_k\vert^2}{\vert A_k\vert^2 }. \end{array}\end{displaymath} (9.25)

where

$\displaystyle T_1(k) = \frac{1}{\rho_k} - \frac{\vert X_k\vert^2}{N\rho_k^2}$ (9.26)

We can simplify (9.25). For example, we can rewrite $ D(a_m)$ as

$\displaystyle D(a_m) = \sum_{k=0}^{N-1} \;
{\rm Re}\left( \frac{e^{j 2\pi k m...
...\bar{A}_k} \right)
(1 - \frac{\vert X_k \vert^2\vert A_k\vert^2}{N\sigma^2}).
$

The $ 1$ in the last term can be ignored since

$\displaystyle \sum_{k=0}^{N-1} {\rm Re}\left( \frac{e^{j 2\pi k m / N}}{\bar{A}_k} \right) = 0
$

because $ 1/\bar{A}_k$ is the DFT of an anti-causal sequence. This leaves

$\displaystyle D(a_m) = - \sum_{k=0}^{N-1} \; {\rm Re}\left( \frac{e^{j 2\pi k m / N} \bar{X}_k \; X_k \; A_k }{N \; \vert B_k\vert^2 \; \sigma^2 } \right).$ (9.27)

Similarly,

$\displaystyle D(b_m) = \sum_{k=0}^{N-1} \; \frac{\vert X_k\vert^2}{N\rho_k} \; ...
...ft( {1\over \bar{B}_k} \; e^{j 2\pi k m / N} \right), \;\;\;\; 1 \leq m \leq Q,$ (9.28)

Using (9.24), (16.9),

$\displaystyle I(\sigma^2, \sigma^2) = \frac{N}{2\sigma^4}$ (9.29)

and

\begin{displaymath}\begin{array}{rcl} I(b_m,b_l) &=& 2 \sum_{k=0}^{N-1} \; \frac...
...m Re}\left( e^{j 2\pi k l / N} / \bar{B}_k \right). \end{array}\end{displaymath} (9.30)

A simple form for $ I(b_m,b_l)$ can be found. We re-write (9.30) as

$\displaystyle I(b_m,b_l) = \sum_{k=0}^{N-1} \; {\rm Re}\left( e^{j 2\pi (m-l) k...
...} / \bar{B}_k^2 \right) \;\;\;\;\; 0 \leq m \leq Q, \;\;\;\;\; 0 \leq l \leq Q.$ (9.31)

Note that the frequency-domain function

$\displaystyle H_k = 1 / \bar{B}_k^2 $

is the DFT of an anti-causal function since both $ B_k$ and $ 1/B_k $ are causal, so

$\displaystyle \sum_{k=0}^{N-1} e^{j 2\pi (m+l) k / N} / \bar{B}_k^2 = 0$

for $ 0 < (l+m) < N-Q$ giving

$\displaystyle I(b_m,b_l) = {\rm Re}\left( \sum_{k=0}^{N-1} \; {e^{j 2\pi (m-l) k / N} \over \vert B_k\vert^2} \right) = N \; r^b_{\vert m-l\vert},$ (9.32)

where $ r^b_\tau$ is the ACF of the $ Q$-th order AR process with parameters

$\displaystyle \sigma^2=1, \;\; [b_1 \ldots b_Q].$

   rb = rlevinson(b,1);
   Ibb = toeplitz(rb(1:Q)) * N;
Similarly,

\begin{displaymath}\begin{array}{rcl} I(a_m,a_l) &=& 2 \sum_{k=0}^{N-1} \; \frac...
...\left( e^{j 2\pi k l / N} / \bar{A}_k \right)   \end{array}\end{displaymath} (9.33)

which is the same form as (9.30), so we can copy the form of (9.32) to arrive at

\begin{displaymath}\begin{array}{rcl} I(a_m,a_l) &=& {\rm Re}\left( \sum_{k=0}^{...
...vert^2} \right)   &=& N\; r^a_{\vert m-l\vert}, \end{array}\end{displaymath} (9.34)

where $ r^a_\tau$ is the ACF of the $ P$-th order AR process with parameters

$\displaystyle \sigma^2=1, \;\; [a_1 \ldots a_P].$

   ra = rlevinson(a,1);
   Iaa = toeplitz(ra(1:P)) * N;
Also,

\begin{displaymath}\begin{array}{rcl} I(a_m,b_l) &=& -2 \sum_{k=0}^{N-1} \; \fra...
...j 2\pi k (m-l) / N} / (B_k \bar{A}_k) \right)   \end{array}\end{displaymath} (9.35)

We also arrive at

$\displaystyle I(\sigma^2, b_m) = \frac{1}{2} \sum_{k=0}^{N-1} \; \frac{\vert B_...
...} \; {\rm Re}\left( \sum_{k=0}^{N-1} \bar{B}_k e^{j 2\pi k m / N} \right) = 0,
$

because $ b_m$ is zero for negative $ m$. Similarly,

$\displaystyle I(\sigma^2, a_m) = 0.
$

(See [25] p. 277 for additional information). A MATLAB implementation for $ I(a_m,a_l)$, $ I(a_m,b_l)$, $ I(b_m,b_l)$ is provided below (See software/pdf_arma_circ.m).
>>   A=fft([a(:); zeros(N-P-1,1)]);
>>   B=fft([b(:); zeros(N-Q-1,1)]);
>>   E=exp(-1i*2*pi/N*[0:N-1]' * [0:max(P,Q)]);
>>   Hb = real( E(:,2:Q+1) .* repmat(1./B,1,Q) );
>>   Ha = real( E(:,2:P+1) .* repmat(-1./A,1,P) );
>>   I.a_a = 2 * Ha' * Ha;
>>   I.a_b = 2 * Ha' * Hb;  
>>   I.b_b = 2 * Hb' * Hb;

Baggenstoss 2017-05-19