I am trying to implement an algorithm in real-time on a Fixed point DSP (The Blackfin from Analog Devices). The algorithm does a lot of stuff, but in the middle it performs an algorithm called "Fast Data Projection Method" (FDPM), which goes something like this:
Let $\mathbf{x}_{k} = [x_{1}, x_{2},\dotsb,x_{M}]^{T}$ be a random vector which contains $M$ samples from a discrete signal. In the Process we take sequentially many vectors from the sampled signal with some degree of overlapping, but that is not relevant. We can assume that the vectors $\mathbf{x}_{k}$ come one after another.
The FDPM aims at obtaining a matrix $W \in \mathbb{R}^{M\times N}$ whose columns are the eigenvectors of the correlation matrix $R_{x} = E[\mathbf{x}\mathbf{x}^{T}]$ of the vector $\mathbf{x}$. So we initiallize the algorithm with a random matrix $W_{0}$ (Can also be the identity matrix), and perform the following steps:
1 - $\mathbf{y}_{k} = W_{k}^{T}\mathbf{x}_{k}$
2 - $\mathbf{a}_{k} = \mathbf{y}_{k} - \|\mathbf{y}_{k}\|\mathbf{e}_{1}$, (Where $\mathbf{e}_{1} = [1,0,0,\dotsb,0]^{T}$)
3 - $G_{k+1} = I-\frac{2}{\|\mathbf{a}_{k}\|^{2}}\mathbf{a}_{k}\mathbf{a}_{k}^{T}$
4 - $W_{k+1} = Normalize\{\left[W_{k}+\mu_{k}\mathbf{x}_{k}\mathbf{x}_{k}^{T}W_{k}\right]G_{k+1}\}$
Where $Normalize\{·\}$ stand for normalizing each of the columns of the matriz individually.
So the problem is the following. I have tested this algoritm in MATLAB with double precisiòn and it works fine, i am able to obtain the eigenvectors very close to the real signal eigenvectors and the algorithms does it's job, no problem.
The thing is, when i use the fixed point toolbox to perform this operations with 16bit numbers in Q15 format, then everything goes wrong. The first one or two iterations have small error, but quickly all the vectors start to saturate (The components go to either 1 or -1) and while the matrices do not saturate, they converge to a matrix which is orthogonal to the same matrix calculated with double precision.
I am guessing that the loss of precision occurred by casting to 16bit affects too much, so i think i might have to do some signal scaling or some "re normalization" every once in a while, but i have no idea how to analize the problem and how to know where i have to modify the algorithm to make it work in 16bit.
Has anyone ever worked with this type of problem? anyone has any idea what i could do?
Thanks!!