michael@0: /* Copyright (c) 2008-2011 Octasic Inc.
michael@0:    Written by Jean-Marc Valin */
michael@0: /*
michael@0:    Redistribution and use in source and binary forms, with or without
michael@0:    modification, are permitted provided that the following conditions
michael@0:    are met:
michael@0: 
michael@0:    - Redistributions of source code must retain the above copyright
michael@0:    notice, this list of conditions and the following disclaimer.
michael@0: 
michael@0:    - Redistributions in binary form must reproduce the above copyright
michael@0:    notice, this list of conditions and the following disclaimer in the
michael@0:    documentation and/or other materials provided with the distribution.
michael@0: 
michael@0:    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
michael@0:    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
michael@0:    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
michael@0:    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
michael@0:    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
michael@0:    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
michael@0:    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
michael@0:    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
michael@0:    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
michael@0:    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
michael@0:    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
michael@0: */
michael@0: 
michael@0: #ifdef HAVE_CONFIG_H
michael@0: #include "config.h"
michael@0: #endif
michael@0: 
michael@0: #include "opus_types.h"
michael@0: #include "opus_defines.h"
michael@0: 
michael@0: #include <math.h>
michael@0: #include "mlp.h"
michael@0: #include "arch.h"
michael@0: #include "tansig_table.h"
michael@0: #define MAX_NEURONS 100
michael@0: 
michael@0: #if 0
michael@0: static OPUS_INLINE opus_val16 tansig_approx(opus_val32 _x) /* Q19 */
michael@0: {
michael@0: 	int i;
michael@0: 	opus_val16 xx; /* Q11 */
michael@0: 	/*double x, y;*/
michael@0: 	opus_val16 dy, yy; /* Q14 */
michael@0: 	/*x = 1.9073e-06*_x;*/
michael@0: 	if (_x>=QCONST32(8,19))
michael@0: 		return QCONST32(1.,14);
michael@0: 	if (_x<=-QCONST32(8,19))
michael@0: 		return -QCONST32(1.,14);
michael@0: 	xx = EXTRACT16(SHR32(_x, 8));
michael@0: 	/*i = lrint(25*x);*/
michael@0: 	i = SHR32(ADD32(1024,MULT16_16(25, xx)),11);
michael@0: 	/*x -= .04*i;*/
michael@0: 	xx -= EXTRACT16(SHR32(MULT16_16(20972,i),8));
michael@0: 	/*x = xx*(1./2048);*/
michael@0: 	/*y = tansig_table[250+i];*/
michael@0: 	yy = tansig_table[250+i];
michael@0: 	/*y = yy*(1./16384);*/
michael@0: 	dy = 16384-MULT16_16_Q14(yy,yy);
michael@0: 	yy = yy + MULT16_16_Q14(MULT16_16_Q11(xx,dy),(16384 - MULT16_16_Q11(yy,xx)));
michael@0: 	return yy;
michael@0: }
michael@0: #else
michael@0: /*extern const float tansig_table[501];*/
michael@0: static OPUS_INLINE float tansig_approx(float x)
michael@0: {
michael@0: 	int i;
michael@0: 	float y, dy;
michael@0: 	float sign=1;
michael@0: 	/* Tests are reversed to catch NaNs */
michael@0:     if (!(x<8))
michael@0:         return 1;
michael@0:     if (!(x>-8))
michael@0:         return -1;
michael@0: 	if (x<0)
michael@0: 	{
michael@0: 	   x=-x;
michael@0: 	   sign=-1;
michael@0: 	}
michael@0: 	i = (int)floor(.5f+25*x);
michael@0: 	x -= .04f*i;
michael@0: 	y = tansig_table[i];
michael@0: 	dy = 1-y*y;
michael@0: 	y = y + x*dy*(1 - y*x);
michael@0: 	return sign*y;
michael@0: }
michael@0: #endif
michael@0: 
michael@0: #if 0
michael@0: void mlp_process(const MLP *m, const opus_val16 *in, opus_val16 *out)
michael@0: {
michael@0: 	int j;
michael@0: 	opus_val16 hidden[MAX_NEURONS];
michael@0: 	const opus_val16 *W = m->weights;
michael@0: 	/* Copy to tmp_in */
michael@0: 	for (j=0;j<m->topo[1];j++)
michael@0: 	{
michael@0: 		int k;
michael@0: 		opus_val32 sum = SHL32(EXTEND32(*W++),8);
michael@0: 		for (k=0;k<m->topo[0];k++)
michael@0: 			sum = MAC16_16(sum, in[k],*W++);
michael@0: 		hidden[j] = tansig_approx(sum);
michael@0: 	}
michael@0: 	for (j=0;j<m->topo[2];j++)
michael@0: 	{
michael@0: 		int k;
michael@0: 		opus_val32 sum = SHL32(EXTEND32(*W++),14);
michael@0: 		for (k=0;k<m->topo[1];k++)
michael@0: 			sum = MAC16_16(sum, hidden[k], *W++);
michael@0: 		out[j] = tansig_approx(EXTRACT16(PSHR32(sum,17)));
michael@0: 	}
michael@0: }
michael@0: #else
michael@0: void mlp_process(const MLP *m, const float *in, float *out)
michael@0: {
michael@0:     int j;
michael@0:     float hidden[MAX_NEURONS];
michael@0:     const float *W = m->weights;
michael@0:     /* Copy to tmp_in */
michael@0:     for (j=0;j<m->topo[1];j++)
michael@0:     {
michael@0:         int k;
michael@0:         float sum = *W++;
michael@0:         for (k=0;k<m->topo[0];k++)
michael@0:             sum = sum + in[k]**W++;
michael@0:         hidden[j] = tansig_approx(sum);
michael@0:     }
michael@0:     for (j=0;j<m->topo[2];j++)
michael@0:     {
michael@0:         int k;
michael@0:         float sum = *W++;
michael@0:         for (k=0;k<m->topo[1];k++)
michael@0:             sum = sum + hidden[k]**W++;
michael@0:         out[j] = tansig_approx(sum);
michael@0:     }
michael@0: }
michael@0: #endif