_NUMERICAL C AND DSP_
by Kevin Leary

Listing One

#include <stdio.h>
#include <math.h>
#include <dspc.h>

static inline float sigmoid (float x){ return 1.f/(1.f+expf(-x));}
/* Simulating a backpropogation network: a feed-forward network. All signals 
move foward from the input to the outputs. xp ip op the 'p' is used to denote
the pth exemplar. Multiply 'xp' our input vector by wts on the hidden layer.
      * * * * * * *  input layer.
      \ | | | | | /
       + + + + + +   hidden layer.
      / | | | | | \
     o  o o o o o  o output layer 
     ip output of the hidden layer.  op output of the output layer.
The matrix wts_h is used to map the input layer to the hidden layer, each node
in the input layer has a weight from that node to every node in hidden layer.
An zero in the weight matrix implies there is no connection. So the ip vector 
or the output of the input layer.
                     h
   ip  = SUM xp * wts
     i    j    i     ji
                     o
   op  = SUM ip * wts
     k    j    j     kj       the sigmoid function is used as a thresholding 
function on the output signal of the layer.*/
void
bpn_simulate(int nni, int nho, int nno, float xp[nni], float ip[nho], 
             float op[nno], float pm wts_h[nho][nni], float pm wts_o[nno][nho])
{
  iter i = nni, j = nho, k = nno;
  ip[j] = sigmoid (sum (xp[i]*wts_h[j][i]));
  op[k] = sigmoid (sum (ip[j]*wts_o[k][j]));
}
/* Training is the crux of the algorithim and where network got its name--
Backpropogation. Compute the error at the output layer and adjust the weights
on the output layer. Take some weighted multiple of the error and propogate 
error back to hidden layer. Given xp the input vector, yp the expected output 
vector, and a set of weights wts_h, wts_o, the control variable eta is used to
control how fast algorithim learns. erro = yp-op is the difference of what we 
expected and what we got. Now we find the gradiant of the error surface and 
build delta_o and delta_h from delta_o. Modify the weights with delta_h and 
delta_o. Return half of the gradiant squared to be used by a control algorithm
to determine how much more we need to train. */
float
bpn_train(int nni, int nho, int nno, float xp[nni], float yp[nno],
                float pm wts_h[nho][nni], float pm wts_o[nno][nho], float eta)
{
  iter i = nni, j = nho, k = nno;
  float s;
  float ip[nho], op[nno], delta_h[nho], delta_o[nno];
  /* simulate xp to get ip, op */
  ip[j] = sigmoid (sum (xp[i]*wts_h[j][i]));
  op[k] = sigmoid (sum (ip[j]*wts_o[k][j]));
  /* error * derivitive of 1/(1+exp(op)) {simplifed} */
  delta_o[k] = (yp[k] - op[k]) * op[k] * (1 - op[k]);
  delta_h[j] = ip[j] * (1 - ip[j]) * sum (delta_o[k] * wts_o[k][j]);
  /* adjust the weights on the output layer. */
  wts_o[k][j] += eta * delta_o[k] * ip[j];
  /* adjust the weights on the hidden layer */
  wts_h[j][i] += eta * delta_h[j] * xp[i];
  /* return half of the gradiant squared  d2 */
  return .5f * sum (delta_o[k] * delta_o[k]);
}
/* Initialize the weights to random values between 0 1. not too important */
void
bpn_wts_init(int nni, int nho, int nno,
                           float pm wts_h[nho][nni], float pm wts_o[nno][nho])
{
  float op[nno],ip[nho];
  for (p) { 
    bpn_simulate (nni,nho,nno,exemplars[p],ip,op,wts_h,wts_o);
#if 0
    printf("{"); printf (" %f,",exemplars[p][i]); printf("}\n");  
    printf("{"); printf (" %f,",desired[p][j]); printf("}\n");  
    printf("{"); printf (" %f,",op[j]); printf("}\n");  
    printf ("------------\n"); 
#endif
  }
}
void
bpn_training_loop(int nni, int nho, int nno,int nex, float exemplars[nex][nni],
   float outputs[nex][nno],float pm wts_h[nho][nni], float pm wts_o[nno][nho],
   float eta,float epsilon,int printp)
{
  float error = 100;
  int    iteration=0;
  iter   p = nex;
  while (error > epsilon)
    {
      error = 0;
      for (p) 
       error += bpn_train (nni, nho, nno,exemplars[p], outputs[p], wts_h, 
                                                                   wts_o, eta);
      error /= nex;
      if (printp)
         printf ("%5d: err=%10lf\n",iteration,error);
      iteration++;
      if (printp && (iteration % printp) == 0)
       bpn_simulate_set (nni, nho, nno, nex, exemplars, outputs, wts_h, wts_o);
    }
}