ld.c

这是数字信号处理方面的一些源码

/***************************************************************************
LEVINSON-DURBIN 算法
****************************************************************************/
#include <stdio.h>
#include <math.h>

#define P 14
#define LENGTH 200

void     durbin( int n, int *buff, float *a  );

main()
{
int             i,buffer[LENGTH];         
float   a[P+1];

        for(i=0;i<LENGTH;i++)
        {
                buffer[i]=(int)(1000*sin(2*3.1415926*i/LENGTH));
        }                                            
        durbin(LENGTH,buffer,a);             
        for(i=0;i<=P;i++)
        {
                printf("a[%2d]=%f\n",i,a[i]);
        }
}

/***************************************************************************
*                                                                          *
*   int   n     : 输入数据缓冲区 buff 中的数据个数                         *
*   int   *buff : 输入数据缓冲区的地址，用于存放被分析的整数数据           *
*   float *a    : LPC 系数的缓冲区首地址，对于P 阶线性预测，缓冲区的长度   *
*                 为 P+1 个浮点数长度。                                    *
*                                                                          *
****************************************************************************/


void     durbin( int n, int *buff, float *a  )
{
int     i , j , m;
long    r[P+1], sum;
float   energy , sum_f;
float   aa1[P+1], aa2[P+1], k[P+1];
float   *p, *p_buf, *p_tmp;

        p=aa1;
        p_buf=aa2;
        for( m=0; m<=P; m++ )       /* caculate the relation R(m) */
        {
                sum=0;
                for( i=0,j=m; j<n; i++,j++)
                {
                        sum=sum+(long)buff[i]*(long)buff[j];
                }
                r[m]=sum;
        }                           /* end of caculate the relation R(m) */

        energy=(float)r[0];         /* initialize the coefficient */
        for( i=1; i<=P; i++)
        {
                p[i]=0; p_buf[i]=0;
        }
        p[0]=1; p_buf[0]=1;

        for( m=1; m<=P; m++ )
        {
                sum_f=0;                /* caculate PARCOR coefficient k[m] */
                for( i=0; i<m; i++ )
                {
                        j=m-i;
                        if( j<0 ) j=-j;
                        sum_f=sum_f+p[i]*r[j];
                }
                k[m]=-sum_f/energy;     /* end of caculate PARCOR coefficient k[m] */

                for( i=1; i<=m; i++ )   /* caculate LPC coefficient */
                {
                        p_buf[i]=p[i]+k[m]*p[m-i];
                }
                energy=( 1-k[m]*k[m] )*energy;
                p_tmp=p;
                p=p_buf;
                p_buf=p_tmp;
        }
        for(m=0; m<=P; m++)             /* end of caculate LPC coefficient */
        {
                a[m]=p[m];
        }
}