# How to Prevent overflow/saturation in fixed point implementation of Fast Data Projection Algorithm

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

• I guess your processor uses higher internal accuracy (e.g. 32 bits) for MAC (multiply-accumulate). Are you sure that this is simulated properly in Matlab? In any case, what you need to find out is where exactly overflow occurs, so you can apply proper scaling to avoid it. That's of course easier said than done ... – Matt L. Mar 25 '14 at 18:05
• Yes of course, the processor has Two 16bit ALUs, and two 40bit accumulators with 8 guard bits. This is simulated with the fixed point toolbox. For example after an inner product of two vectors, when all the intermediate sums are done agains the 40bit Acc, the final resulted has to casted down to 16bit and when it was over 1 or -1, the processor saturates by hardware. I know "where" it occurs, what i don't know is how to stop it. – bone Mar 25 '14 at 18:08
• By "where" I meant where in the algorithm, of course. Divisions are dangerous; are you sure that the numerator never gets bigger than the denominator? Anyway, the only solution is proper scaling. – Matt L. Mar 25 '14 at 18:12
• Haha yes yes, that is what i meant. For exmaple the vector $\mathbf{y}_{k}$ starts to slowly drift appart from its double precisiòn counterpart and eventually all the components saturate, and with them the rest of the algorithm. I am still not sure which is the original cause of the saturation. Anyway, thanks for the answer, at least i know that i have to try millions of different ways of scaling until it works. – bone Mar 25 '14 at 18:16
• Since you have a floating-point version of it: keep track of the intermediate values and see when the values go out of range (or consistently are small). That would provide a decent starting point to figure out which numbers are more likely to need scaling. – Oscar Mar 27 '14 at 16:33