Saturday 25 June 2016

Exact Formula for Asymptotic Convergence of Fourier Transform of Uniform Random Variables

It is well known that the distribution of sums of random variables converges to Gaussian. In this post, we look look at the convergence of Fourier transform of uniformly distributed random variables.

Exact and analytical formulas for convergence are interesting as they help develop an intuition about the asymptotic convergence phenomenon and understand the convergence error for finite number  of samples. Under the scanner today is the convergence behaviour of Fourier transform of uniformly distributed random variables (URVs). 

Let $u_1, u_2, \cdots, u_N$ denote complex-valued zero mean i.i.d URVs with finite support $[-q_n, q_n]\times [-q_n, q_n]$ for $n = 1\cdots N$. Their Fourier transform is given by
\begin{align}U_k = \frac{1}{\sqrt{N}} \sum_{n = 1}^{N} u_n W^{kn} = \frac{1}{\sqrt{N}} \sum_{n = 1}^{N} v_n,\end{align}
where $W^{kn}$ are the roots of unity, and $v_n = u_n W^{kn}$.

Let $U^r_k = \mathfrak{Re}(U_k)$, and $v^r_n = \mathfrak{Re}(v_n)$. Variable $v^r_n$ preserves the distribution of $\mathfrak{Re}(u_n)$. Therefore, we have the characteristic function of the real part of the Fourier transformed variable as
\begin{align}\text{E}(e^{itU_k^r}) &= \prod_{n = 1}^{N} \text{E}(e^{\frac{1}{\sqrt{N}}it v^r_n})\nonumber \\ &= \prod_{n = 1}^{N} \int_{-q_n}^{q_n} \frac{1}{2q_n} e^{\frac{1}{\sqrt{N}}it v^r_n} dv^r_n \nonumber\\ &= \prod_{n = 1}^{N} \frac{e^{\frac{1}{\sqrt{N}}itq_n} - e^{-\frac{1}{\sqrt{N}}itq_n}}{\frac{2}{\sqrt{N}}itq_n} \nonumber\\ &= \exp\left( \sum_{n = 1}^{N} \log \frac{\sin\left(\frac{1}{\sqrt{N}}tq_n\right)}{\frac{1}{\sqrt{N}}tq_n}\right).\end{align}

Interestingly, for $f(t) = \log \frac{\sin(t\alpha)}{t\alpha}$, the limits of the derivatives $\lim_{t\to 0} \frac{\text{d}^n f(t)}{\text{d}t^n}$, for $n = 1, 2, \cdots$, exist! Therefore, the Taylor expansion of $f(t)$ around $t=0$ is given by
\begin{align}f(t) = -\frac{t^2\alpha^2}{6} - \frac{t^4\alpha^4}{180} - \frac{t^6\alpha^6}{2835} + \cdots.\end{align}
The Taylor series expansion given above may be easily verified using symbolic solvers like Wolfram Alpha Online.

Substituting the Taylor series expansion in the characteristic function above, we have (only the first two terms of the Taylor series are shown)
\begin{align}\text{E}(e^{itU_k^r}) = \exp\left(-\frac{1}{6} \frac{t^2 \sum_{n = 1}^{N} q^2_n}{N} - \frac{1}{180}\frac{t^4 \sum_{n = 1}^{N} q^4_n}{N^2} + \cdots\right). \label{eqn:cg}\end{align}

Therefore, as $N\to \infty$, 
\begin{align}\lim_{N\to \infty} \text{E}(e^{itU_k^r}) = \exp\left(-\frac{t^2}{2} \frac{\bar{q}^2}{3}\right), \label{eqn:acg} \end{align}
where
\begin{align}\bar{q}^2 = \lim_{N\to \infty} \frac{1}{N}\sum_{n = 1}^{N} q_n^2.\end{align}

The right hand side expression in (\ref{eqn:acg}) showing the asymptotic convergence of the characteristic function of $U^r_k$ is readily identified as the characteristic function of Gaussian distribution with zero mean and variance $\bar{q}^2/3$. A similar expression may also be found for $\lim_{N\to \infty} \text{E}(e^{itU_k^i})$ in a straightforward manner.

Therefore, the exact convergence formula for the Fourier transform of URVs is as shown in (\ref{eqn:cg})!

ARK


Saturday 30 January 2016

Complex-valued Durbin Algorithm with MATLAB code

This post is a rehash of my earlier post [1] where the complex-valued Durbin's algorithm was used indirectly for matrix inversion. Here, the algorithm is written down separately as it may be useful for some to have it "out of the box."

The complex-valued versions of Durbin's and Levinson's algorithms may be arrived upon by straightforward application of description in [2] sec. 4.7.2 "Solving the Yule-Walker Equations", and 4.7.3 "The General Right Hand Side Problem".

Let $r_0 = 1, r_1, r_2, \cdots, r_n \in \mathbb{C}$ and $T \in \mathbb{C}^{\text{N}\times\text{N}}$ be a positive definite Hermitian symmetric Toeplitz matrix constructed with the scalars $r_0, r_1, \cdots, r_{n-1}$ [see example below], Durbin's algorithm solves the system $Ty=-(r_1, r_2, \cdots, r_n)$. The algorithm is as follows.

$z(1) = -r_1^*; \beta=1; \alpha=-r_1^*$
$\text{for}~ k = 1:n-1$
    $\beta = (1-|\alpha|^2)\beta$
    $\alpha = - (r^*(k+1) - r(k:-1:1)^H z(1:k))/\beta$
    $z(1:k) = z(1:k) + \alpha z^*(k:-1:1)$
    $z(k+1) = \alpha$
$\text{end}$
$y = z^*$

MATLAB code:

>> r = complex(rand(10,1),rand(10,1)) ;
>> r(1) = 1 ;
>> T = toeplitz(r(1:9),conj(r(1:9))) ;
>> k = r(2:end) ;
>> y1 = inv(T)*-k ;
>> norm(y1-durbin(r))
ans = 1.0709e-10

function Y = durbin(R)
% Y = durbin(R)
% where R is a vector with r0, r1, r2, ..., rn.
% Solves RY = -K, where R is the NxN symmetric Toeplitz matrix 
% with elements r0, ..., r(n-1) and K is a vector with elements
% r1, ..., rn.
%

R = R/R(1,1) ;
N = length(R)-1 ;
R = R(2:end) ;
Z = zeros(N,1) ;

Z(1) = -conj(R(1)) ;
A = -conj(R(1)) ;
B = 1 ;

for k=1:N-1
    B = (1-abs(A)^2)*B ;
    A = -(conj(R(k+1)) + R(k:-1:1)'*Z(1:k))/B ;
    Z(1:k) = Z(1:k) + A*conj(Z(k:-1:1)) ;
    Z(k+1) = A ;
end
Y = conj(Z) ;

References:
[1] Inverse Cholesky Factor of a positive definite Hermitian symmetric Toeplitz matrix using Durbin recursions, Applied Mathematics and Engineering Blog, September 2015.
[2] Gene H. Golub, Charles F. Van Loan, Matrix Computations, Third Edition.