/* Midpoint Fraying and Splicing */

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define N	5		/* Half-signal length */
#define E	3		/* Must be less than N */

typedef float real;

#ifndef PI
# define PI	3.14159265358979323846264338327950288
# define PI2	1.57079632679489661923132169163975144
# define PI4	0.78539816339744830961566084581987572
#endif


/* Print a vector of reals as one row. */
static void
print_vector( 
	     const char *title,
	     real *x,
	     int n)
{
  int i;
  printf("%s (dim=%d): ", title, n);
  for(i=0; i<n; i++ ) printf("%7.4f ", x[i]);
  putchar('\n');
  return;
}

/* Evaluate the rising cutoff function r[1](t) at a point t. */
static real
r1( real t )
{
  real r=0;
  if( t >  -1.0 )
    if ( t >= 1.0 ) r = 1;
    else r = sin(0.25*PI*(sin(0.5*PI*t)+1.0));
  return r;
}


/* Fraying a sampled signal at the point between uneg[-1] and upos[0]: */
void
fray(
     real *uneg,		/* Valid indices: -e,...,-2,-1 */
     real *upos,		/* Valid indices:  0,1,...,e-1 */
     int e)			/* Reach of the indices */
{
  real rp, rn, up, un, t;  int k;
  assert(e>=0);  assert(uneg);  assert(upos);
  for(k=0; k<e; k++) {
    t = (k+0.5)/(real)e;    rp = r1(t);    rn = r1(-t);
    up = upos[k];    un = uneg[-1-k];
    upos[k] = rp*up + rn*un;    uneg[-1-k] = rp*un - rn*up;
  }
  return;
}

/* Splicing a sampled signal at the point between uneg[-1] and upos[0]: */
void
splice(
       real *uneg,		/* Valid indices: -e,...,-2,-1 */
       real *upos,		/* Valid indices:  0,1,...,e-1 */
       int e)			/* Reach of the indices */
{
  real rp, rn, up, un, t;  int k;
  assert(e>=0);  assert(uneg);  assert(upos);
  for(k=0; k<e; k++) {
    t = (k+0.5)/(real)e;    rp = r1(t);    rn = r1(-t);
    up = upos[k];    un = uneg[-1-k];
    upos[k] = rp*up - rn*un;    uneg[-1-k] = rp*un + rn*up;
  }
  return;
}

int
main(void)
{
  real u[2*N], v[2*N], rise[2*E];
  int k, k1;

  /* Show the rising cutoff function: */
  for(k=0, k1=(-E); k1<E; k++, k1++)
    rise[k] = r1((0.5+k1)/(real)E);
  print_vector("Rising cutoff function", rise, 2*E);

  /* Fill u[],u1[],u2[] with a smooth function: */
  for(k=0; k<2*N; k++)
    u[k] = v[k] = sin(0.3*PI*(k+0.5)/N);
  print_vector("Original signal", u, 2*N);
    
  /* Fray `u[]' at index N-1/2, then splice it back together: */
  fray( &u[N],  &u[N], E );
  print_vector("  Frayed signal", u, 2*N);
  splice( &u[N],  &u[N], E );
  print_vector("Restored signal", u, 2*N);

  /* Splice `v[]' at index N-1/2, then fray it back together: */
  splice( &v[N],  &v[N], E );
  print_vector(" Spliced signal", v, 2*N);
  fray( &v[N],  &v[N], E );
  print_vector("Restored signal", v, 2*N);
  


  exit(EXIT_SUCCESS);
}