diff options
Diffstat (limited to '2.3-1/src/c/signalProcessing/lev/sleva.c')
-rw-r--r-- | 2.3-1/src/c/signalProcessing/lev/sleva.c | 66 |
1 files changed, 66 insertions, 0 deletions
diff --git a/2.3-1/src/c/signalProcessing/lev/sleva.c b/2.3-1/src/c/signalProcessing/lev/sleva.c new file mode 100644 index 00000000..87c24660 --- /dev/null +++ b/2.3-1/src/c/signalProcessing/lev/sleva.c @@ -0,0 +1,66 @@ +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008 - INRIA - Arnaud TORSET + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + + +/*Resolve the Yule-Walker equations: + + |r(0) r(1) ... r(N-1)|| a(1) | |sigma2| + |r(1) r(0) ... r(n-1)|| a(2) | | 0 | + | : : ... : || : |=| 0 | + | : : ... : || : | | 0 | + |r(N-1) r(N-2) ... r(0) ||a(N-1)| | 0 | + + using Levinson's algorithm. + r :Correlation coefficients + ar :Auto-Regressive model parameters + sigma2 :Scale constant + rc :Reflection coefficients +*/ + +#include <stdlib.h> +#include "lev.h" + +float sleva(float* in,int size, float* ar, float* rc){ + int i=0, j=0; + float accu=0; + float* ak1; + float sigma2; + + /* FIXME : malloc here */ + ak1=(float*)malloc((unsigned int)size*sizeof(float)); + + /* initialize levinson's algorithm */ + ar[0]=-in[1]/in[0]; + rc[0]=ar[0]; + sigma2=(1-ar[0]*ar[0])*in[0]; + ak1[0]=0; + + /* iterative solution to yule-walker equations */ + for (i=1;i<size-1;i++){ + accu=0; + for (j=0;j<i;j++){ + accu+=ar[j]*in[i-j]; + } + ak1[i]=-(in[i+1]+accu)/sigma2; + rc[i]=ak1[i]; + for (j=0;j<i;j++){ + ak1[j] = ar[j]+ak1[i]*ar[i-1-j]; + } + sigma2=(1-ak1[i]*ak1[i])*sigma2; + for (j=0;j<=i;j++){ + ar[j]=ak1[j]; + } + } + free(ak1); + return sigma2; +} + |