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

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

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

$$\rho_k = \sigma^2 \; \frac{ \vert B_k\vert^2}{\vert A_k\vert^2}.$$ (10.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. 10.22 with $2L$ replaced by $N$).

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

Thus,

$$\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}$$ (10.25)

where

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

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

$$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

$$\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

$$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).$$ (10.27)

Similarly,

$$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,$$ (10.28)

Using (10.24), (17.9),

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

and

$$\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}$$ (10.30)

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

$$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.$$ (10.31)

Note that the frequency-domain function

$$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

$$\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

$$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},$$ (10.32)

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

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

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

Similarly,

$$\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}$$ (10.33)

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

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

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

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

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

Also,

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

We also arrive at

$$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,

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

(See [31] 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;