Thursday, November 6, 2014

Principal Components Analysis (PCA)

To understand Principal Components Analysis, take the age and height of trees in a forest as an example, the data is distributed as follow:

It's clear that you don't need 2 dimensions (height and age) to be a representative of a particular tree. If you take an another representative $ha=\frac1{\sqrt{2}}(height+age)$, you then need only 1D to represent trees.
Or you can understand that you want to find the most important characteristics of the trees in this statistic. It may be $ha=\frac1{\sqrt{2}}(height+age)$
Or in mathematic, just answer the simple question: Whether we can reduce the number of dimension of the statistic to keep the most accurate properties of the statistic.
In practical, we have much more dimensions than the above example, in that case, the dimensional reduction is very necessary.
In essential, they are the CHANGE OF BASIS
The old (observed) basis is:
$$B=\begin{bmatrix}b_1\\b_2\\...\\b_m\end{bmatrix}=\begin{bmatrix}1&0&...&0\\0&1&...&0\\...\\0&0&...&1\end{bmatrix}=I$$
We would like to change the basis to another basis by a linear combination:
$$B_P=PB=P$$
The rows of P are a set of new basic vectors.
MATRIX P IS CALLED AS THE PRINCIPAL COMPONENTS OF X
The data is re-expressed as
$$Y=PX$$
That lead questions:
- What is the best way to re-express X?
- What is a good choice of basis P
It's obvious that A and B variable is uncorrelated if and only if AB=0. In practice, A and B is expressed as a set of observed data. Therefore  A and B are uncorrelated if and only if  $A^TB=0$
ANSWER THE QUESTIONS:
The BEST WAY to re-express X is that
$$CY=\begin{bmatrix}y_{11}&0&...&0\\0&y_{22}&...&0\\...\\0&0&...&y_{mm}\end{bmatrix}$$
It means that the re-expressed data is uncorrelated with the data along the new basic vectors
$$C_Y=\frac1mYY^T=\frac1mPXX^TP^T=PC_XP^T$$
$C_X=\frac1mXX^T$ is a symmetric matrix. It can be analyzed as $C_X=EDE^T$ where D is diagonal matrix and $\begin{bmatrix}e_1&e_2&...&e_m\end{bmatrix}$ where the column $e_i$ is eigenvector $e_i$.
$C_X$ is now become $C_Y=(PE)D(PE)^T$ We would like to find the matrix P according to the diagonal matrix C_Y. If we choose $P=E^T$, $C_Y$ will become $C_Y=D$ because E is eigenvector ($EE^T=E^TE=I$).
HOW TO SOLVE PCA
1. Eigenvector decomposition
- The principle components of X are the eigenvectors of $C_X=\frac1nXX^T$
- The $i^{th}$ diagonal value of $C_Y$ is the variance of X along $p_i$
2. Using Single Value Decomposition SVD
If X is an arbitrary nxm matrix and rank of  $X^TX$ is r, I will prove that the principle components of X will be the set of othorgonal eigenvector ${v_1,v_2, ..., v_r}$
Prove
Because $v_i$ is the eigenvector of X associated eigenvalue $\lambda_i$
$$(X^TX)v_i=\lambda_iv_i$$
If we denote $u_i=\frac1\sigma_iXv_i$ where $\sigma_i=\sqrt{\lambda_i}$
It derive to the column $i^{th}$: $Xv_i=\sigma_iu_i$, and in matrix view scope:
$$XV=U\Sigma$$
Where:
$$\Sigma=\begin{bmatrix}\sigma_1&0&...&0&...&0\\0&\sigma_2&...&0&...&0\\...\\0&0&...&\sigma_r&...&0\\0&0&...&0&...0\end{bmatrix}$$
$$V=\begin{bmatrix}v_1&v_2&...&v_m\end{bmatrix}$$
$$U=\begin{bmatrix}u_1&u_2&...&u_n\end{bmatrix}$$
The matrix V and U is added with (m-r) and (n-r) orthonormal vectors respectively. The properties of U and V is:
- The columns of V are orthonormal together (definition) and U are also orthonormal ($u_i^Tu_j=0$
- The $||v_i||=1$ (definition) and $||u_i||=1$
In conclusion, an arbitrary matrix n x m can be analyses as:
$$X=U{\Sigma}V^T$$
It's obvious that
$$X^TX=V{\Sigma}^TU^TU{\Sigma}V^T=VDV^T$$
Therefore, we can compute SVD of $X^TX$ to find the PC of X

Source code to test PCA on Matlab
The input data (X)

Output data (after projection) (Eigen Decomposition method)

Output data (SVD method)


 %Create the samples  
 X=[1:1000;1:1000];    % X is 2 dimensional with 1000 samples  
 X=X+100*randn(2,1000);  % Create some 'noise' around the true value  
 plot(X(1,:),X(2,:),'or'); %Visualize the samples  
 grid on;  
 axis([-2000 2000 -2000 2000])  
 %OBJECTIVE:  
 %FINDING THE PC OF X  
 %PROJECT X ONTO A NEW BASIC (FINDING Y) WHERE IT IS 'BETTER'  
 %1st method (EIGEN VECTOR DECOMPOSITION)  
 N=size(X,2); %Determine the number of samples  
 covX=1/(N-1)*X*X';  
 [PC V]=eig(covX);  
 %Extract diagnonal of matrix as a vector  
 %Sort the variances in decreasing order  
 diagV=diag(V);  
 [tmp,srt]=sort(-1*diagV);  
 V=diagV(srt);  
 PC=PC(:,srt);  
 Y=PC'*X;  
 figure;  
 plot(Y(1,:),Y(2,:),'o'); grid on;  
 axis([-2000 2000 -2000 2000])  
 %2nd method (SINGLE VALUE DECOMPOSITION)  
 mn=mean(X,2);  
 X=X-repmat(mn,1,N); %Remove the mean for each dimension  
 [u,S,PC]=svd(X');  
 S=diag(S);  
 V2=S.*S;  
 Y2=PC'*X;  
 figure;  
 plot(Y2(1,:),Y2(2,:),'o'); grid on;  
 axis([-2000 2000 -2000 2000])  

Reference 

[1] Shlens, Jonathon. "A tutorial on principal component analysis." arXiv preprint arXiv:1404.1100 (2014).

Thursday, October 23, 2014

Derivatives of Matrix to Vector/Scalar


$$y=Hx+n$$
$$y=\begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}=\begin{bmatrix}a_{11}&a_{12}&...&a_{1n}\\a_{21}&a_{22}&...& a_{2n}\\...\\a_{n1}&a_{n2}&...& a_{nn}\end{bmatrix}\begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}+\begin{bmatrix}n_1\\n_2\\...\\n_n \end{bmatrix}$$
$$\min_x((y-Hx)^2)$$

Wednesday, October 22, 2014

WHITE NOISE


White noise is defined to be a stationary random process having a constant spectral density function:
$$S_{wn}(j\omega)=\int_{-\infty}^{\infty}R_X(\tau)e^{-j\omega\tau}d(\tau)=A=const$$
where $$R_X(\tau)=\int_{-\infty}^{\infty}x(t)x(t+\tau)$$

Sunday, October 19, 2014

CW Interference

This post shows the theory of CW interference and Matlab source codes for validating the theory.
- The received satellite signal can be represented as
$$s(t)=A_sC(t)D(t)cos(2\pi(f_c+f_d)t+\Phi_0)$$

Tuesday, October 14, 2014

How to determine C/N0 and SNR


In this tutorial we will show a technique to determine CN0 and the behind story
-Suppose that I(m) and Q(m) is the results of tracking phase at mth millisecond
$$I(m)=\sum_{n=0}^{N-1}(s(m,n)+w(m,n))cos(2\pi f_{IF}n)PRN(n)$$
$$Q(m)=\sum_{n=0}^{N-1}(s(m,n)+w(m,n))sin(2\pi f_{IF}n)PRN(n)$$
We have
$$P_s=(\frac{1}{M}\sum_{m=0}^{M-1}I(m))^2$$
$$P_{total}=\frac{1}{M}\sum_{m=0}^{M-1}|I(m)+jQ(m)|^2$$
$$P_n=P_{total}-Ps$$
$$CN_0=10log_{10}(\frac{2P_sF_s}{P_nN})$$

Friday, September 26, 2014

Realtime GNSS SDR Receiver

In this post, i will describe the most struggles of implementing GNSS-SDR RX and how to overcome the limitations which based on our experiments in this field
1. PC-BASED RX
1.1 Analysis of the limitations
1.2 Solutions
2. FPGA-BASED RX

Wednesday, September 24, 2014

SDRs

1. PC Based Receiver
- Most popular and cited: AKOS Matlab source code, GNSS-SDR (Offline and Realtime), GNSS-SDRlib (very understandable source code -- need to more reviews)
2. FPGA Based Receiver
HelloGNSS
Namuru
SwiffNav
"Build your own GPS Receiver"

Frontends

1. SDRNAV00 (Michele)
http://michelebavaro.blogspot.it/
2. SIGE GN3S
https://www.sparkfun.com/products/10981
3. HG (HelloGNSS)
hellognss.com
4. USRP
5. RTL32U (TCXO added Version)
http://www.rtl-sdr.com/tag/rtl2832/
6. Updating ...

Beamforming Principle

1. What is beamforming
Beamforming