/* Filename: decimator8to1.c
 * Description: Efficient implementation of the multiplierless 8-to-1
 *              decimator.
 * Author: Juha Yli-Kaakinen
 * References: 
 * [1] J. Yli-Kaakinen and T. Saramäki, "A Systematic Algorithm for
 * the Synthesis of Multiplierless Lattice Wave Digital Filters,"
 * invited book chapter in Digital filters, Fausto Pedro García
 * Márquez (Ed.), ISBN: 978-953-307-190-9, Intech, pp. 257-289,
 * Apr. 2011.
 * [2] Gazsi, L, "Explicit formulas for lattice wave digital filters,"
 * IEEE Trans. Circuits Syst., Vol. CAS-32, No. 1, pp. 68--88, 1985.
 * 
 * (c) Juha Yli-Kaakinen,  6.2.2012
 */

#include <stdio.h>

#define BUFFER_SIZE 64  
#define SAMPLE short  

void decimator(short *buffer) {
  unsigned int n;
  // Internal variables
  short inp1, inp2;
  short D;
  // Delay elements
  static short d11 = 0, d21 = 0, d22 = 0;
  static short d31 = 0, d32 = 0, d33 = 0;

  /* First stage (1st-order filter) */
  for (n = 0; n < BUFFER_SIZE>>1; ++n) {
    // ========== 1st brach ==========    
    inp1 = buffer[(n<<1)+1]>>1;
    // 1st adaptor : Type 3 [a = -2^(-1)+2^(-3)+2^(-5)+2^(-8)]
    D = inp1 - d11; 
    D = (D>>1)-(D>>3)-(D>>5)-(D>>8);
    buffer[n] = D - d11 + (buffer[(n<<1)]>>1);     
    d11 = D - inp1; 
  }
 
  /* Second stage (3rd-order filter) */
  for (n = 0; n < BUFFER_SIZE>>2; ++n) {
    // ========== 2nd brach ==========
    // 1st adaptor : Type 4 [a = -2^(-1)-2^(-4)]
    D = d22 - (buffer[(n<<1)]>>1);
    d22 = (D>>1)-(D>>4) - d22; 
    buffer[n] = d22 - D;     
 
    // ========== 1st brach ==========
    inp1 = buffer[(n<<1)+1]>>1;
    
    // 1st adaptor : Type 3 [a = -2^(-3)]
    D = inp1 - d21; 
    D = D>>3;
    buffer[n] += D - d21;
    d21 = D - inp1;
  }

  /* Third stage (5th-order filter) */
  for (n = 0; n < BUFFER_SIZE>>3; ++n) {
    // ========== 2nd brach ==========
    inp1 = buffer[(n<<1)]>>1;
  
    // 1st adaptor : Type 3 [a = -2^(-2)-2^(-4)]
    D = inp1 - d33;
    D = (D>>2)+(D>>4);
    buffer[n] = D - d33;
    d33 = D - inp1; 
  
    // ========== 1st brach ==========
    inp1 = buffer[(n<<1)+1]>>1;
  
    // 1st adaptor : Type 3 [a = -2^(-4)-2^(-6)]
    D = inp1 - d31;
    D = (D>>4)+(D>>6);
    inp2 = D - d31;
    d31 = D - inp1; 

    // 2nd adaptor : Type 4 [a = 1+2^(-2)+2^(-5)+2^(-7)]
    D = d32 - inp2; 
    d32 = (D>>2)+(D>>5)+(D>>7) - d32; 
    buffer[n] += d32 - D;
  }
}
 
int main(){
  unsigned int counter = 0;
  char *inputFileName = "data/sines.bin";
  char *outputFileName = "data/output.bin";
  FILE *inputFile = fopen(inputFileName, "rb");
  FILE *outputFile = fopen(outputFileName, "wb");
  short buffer[BUFFER_SIZE];

  if (!inputFile) {
    fprintf(stderr, "Unable to open file %s\n", inputFileName);
    return 1;
  }

  if (!outputFile) {
    fprintf(stderr, "Unable to open file %s\n", outputFileName);
    return 1;
  }

  while (fread(buffer, sizeof(short), BUFFER_SIZE, inputFile) > 0) {
    decimator(buffer);
    fwrite(buffer, sizeof(short), BUFFER_SIZE>>3, outputFile);
    counter++; 
  }
  printf("Processed %i samples\n", counter*BUFFER_SIZE);

  fclose(inputFile);
  fclose(outputFile);
  
  return 0;
}