math_lib.c

适用于单片机等mcu的FFT程序

/*

2.4 kbps MELP Proposed Federal Standard speech coder

Fixed-point C code, version 1.0

Copyright (c) 1998, Texas Instruments, Inc.  

Texas Instruments has intellectual property rights on the MELP
algorithm.  The Texas Instruments contact for licensing issues for
commercial and non-government use is William Gordon, Director,
Government Contracts, Texas Instruments Incorporated, Semiconductor
Group (phone [phone]).

The fixed-point version of the voice codec Mixed Excitation Linear
Prediction (MELP) is based on specifications on the C-language software
simulation contained in GSM 06.06 which is protected by copyright and
is the property of the European Telecommunications Standards Institute
(ETSI). This standard is available from the ETSI publication office
tel. [phone]. ETSI has granted a license to United States
Department of Defense to use the C-language software simulation contained
in GSM 06.06 for the purposes of the development of a fixed-point
version of the voice codec Mixed Excitation Linear Prediction (MELP).
Requests for authorization to make other use of the GSM 06.06 or
otherwise distribute or modify them need to be addressed to the ETSI
Secretariat fax: [phone].

*/
#include <stdio.h>
#include <math.h>
#include "spbstd.h"
#include "mathhalf.h"
#include "mathdp31.h"
#include "math_lib.h"
#include "constant.h"


Shortword DEBUG;

/***************************************************************************
 *
 *   FUNCTION NAME: sin_fxp
 *
 *   PURPOSE:
 *
 *     Compute the sine of x whose value is expressed in radians/PI.
 *
 *
 *   INPUTS:
 *
 *     x
 *                     16 bit short signed integer (Shortword) in Q15.
 *
 *   OUTPUTS:
 *
 *     none
 *
 *   RETURN VALUE:
 *
 *     ty
 *                     16 bit short signed integer (Shortword) in Q15.
 *
 *************************************************************************/

#define PI_Q13     24576    /* M_PI*(1<<13) */

Shortword sin_fxp(Shortword x)
{
    static Shortword table[129] =
    { 
            0,   402,   804,  1206,  1608,  2009,  2411,  2811,  3212,
	 3612,  4011,  4410,  4808,  5205,  5602,  5998,  6393,  6787,
	 7180,  7571,  7962,  8351,  8740,  9127,  9512,  9896, 10279,
	10660, 11039, 11417, 11793, 12167, 12540, 12910, 13279, 13646,
	14010, 14373, 14733, 15091, 15447, 15800, 16151, 16500, 16846,
	17190, 17531, 17869, 18205, 18538, 18868, 19195, 19520, 19841,
	20160, 20475, 20788, 21097, 21403, 21706, 22006, 22302, 22595,
	22884, 23170, 23453, 23732, 24008, 24279, 24548, 24812, 25073,
	25330, 25583, 25833, 26078, 26320, 26557, 26791, 27020, 27246,
	27467, 27684, 27897, 28106, 28311, 28511, 28707, 28899, 29086,
	29269, 29448, 29622, 29792, 29957, 30118, 30274, 30425, 30572,
	30715, 30853, 30986, 31114, 31238, 31357, 31471, 31581, 31686,
	31786, 31881, 31972, 32058, 32138, 32214, 32286, 32352, 32413,
	32470, 32522, 32568, 32610, 32647, 32679, 32706, 32729, 32746,
	32758, 32766, 32767
    };
    
    Shortword tx, ty;
    Shortword sign;
    Shortword index1,index2;
    Shortword m;
    Shortword temp;

    sign = 0;
    if (x < 0) 
	{
        tx = -x;     
		sign = -1;     
    }
    else 
	{
		tx = x;     
    }

    if (tx > X05_Q15)
    {
        tx = sub(ONE_Q15,tx);     
    }
    //将输入转换成0~128范围内
    index1 = shr(tx,7);     
    index2 = add(index1,1);     

    if (index1 == 128) 
	{
		if (sign != 0)
			return(-(table[index1]));
		else
			return(table[index1]);
    }
    
	
    m = sub(tx,shl(index1,7));
    m = shl(m,8);     

    temp = sub(table[index2],table[index1]);
    temp = mult(m,temp);
    ty = add(table[index1],temp);     

    if (sign != 0)
        return(-ty);
    else
        return(ty);
} /* sin_fxp */


/***************************************************************************
 *
 *   FUNCTION NAME: cos_fxp
 *
 *   PURPOSE:
 *
 *     Compute the cosine of x whose value is expressed in radians/PI.
 *
 *
 *   INPUTS:
 *
 *     x
 *                     16 bit short signed integer (Shortword) in Q15.
 *
 *   OUTPUTS:
 *
 *     none
 *
 *   RETURN VALUE:
 *
 *     ty
 *                     16 bit short signed integer (Shortword) in Q15.
 *
 *************************************************************************/
Shortword cos_fxp(Shortword x)
{
    static Shortword table[129] =
    { 
        32767, 32766, 32758, 32746, 32729, 32706, 32679, 32647, 32610,
        32568, 32522, 32470, 32413, 32352, 32286, 32214, 32138, 32058,
        31972, 31881, 31786, 31686, 31581, 31471, 31357, 31238, 31114,
        30986, 30853, 30715, 30572, 30425, 30274, 30118, 29957, 29792,
        29622, 29448, 29269, 29086, 28899, 28707, 28511, 28311, 28106,
        27897, 27684, 27467, 27246, 27020, 26791, 26557, 26320, 26078,
        25833, 25583, 25330, 25073, 24812, 24548, 24279, 24008, 23732,
        23453, 23170, 22884, 22595, 22302, 22006, 21706, 21403, 21097,
        20788, 20475, 20160, 19841, 19520, 19195, 18868, 18538, 18205,
        17869, 17531, 17190, 16846, 16500, 16151, 15800, 15447, 15091,
        14733, 14373, 14010, 13646, 13279, 12910, 12540, 12167, 11793,
        11417, 11039, 10660, 10279,  9896,  9512,  9127,  8740,  8351,
        7962,   7571,  7180,  6787,  6393,  5998,  5602,  5205,  4808,
        4410,   4011,  3612,  3212,  2811,  2411,  2009,  1608,  1206,
        804,     402,     0
    };
    
    Shortword tx, ty;
    Shortword sign;
    Shortword index1,index2;
    Shortword m;
    Shortword temp;

    sign = 0;
 
    if (x < 0) 
	{
		tx = -x;     
    }
    else 
	{
        tx = x;     
    }

    if (tx > X05_Q15)
    {
        tx = sub(ONE_Q15,tx);     
        sign = -1;     
    }
	//将输入转换成0~128范围内
    index1 = shr(tx,7);     
    index2 = add(index1,1);     
    
    if (index1 == 128)
        return((Shortword)0);
    m = sub(tx,shl(index1,7));
    m = shl(m,8);     

    temp = sub(table[index2],table[index1]);
    temp = mult(m,temp);
    ty = add(table[index1],temp);     

    if (sign != 0)
        return(-ty);
    else
        return(ty);
}