polyphas.c

FIR2LIFT is a program to factor a wavelet FIR filter pair into lifting steps.

/*
FILE : POLYPHAS.C

Routines to handle polyphase matrices.

(C) C. Valens

Created     : 14/09/1999
Last update : 24/09/1999
*/


#include "redir.h"
#include "zpoly.h"
#include "polyphas.h"
#include <stdlib.h>
#include <stdio.h>


/*
 * Print a polyphase pair.
 */
void polyphase_write(zpoly_pair evenodd)
{
  redir_printf("{ ");
  zpoly_write(evenodd.p1);
  redir_printf("} + z^(-1) { ");
  zpoly_write(evenodd.p2);
  redir_printf("}");
  redir_printf("\n");
}


/*
 * Split a polynomial in z in its polyphase components:
 *
 *   x(z) = { x_even(z^2) } + z^(-1) { x_odd(z^2) }
 */
zpoly_pair polyphase_split(zpoly *poly)
{
  zpoly_pair evenodd;
  zpoly *pt1;
  zpoly *pt2;

  evenodd.p1 = zpoly_copy(poly);
  evenodd.p2 = zpoly_copy(poly);

  /*
   * Remove even powers from odd part.
   * Also add 1 to odd powers.
   */
  pt1 = evenodd.p2;
  pt2 = NULL;
  while (pt1) {
    if ((pt1->power&0x1)==0) {
      if (pt2==NULL) {
        evenodd.p2 = pt1->next;
        free(pt1);
        pt1 = evenodd.p2;
      }
      else {
        pt2->next = pt1->next;
        free(pt1);
        pt1 = pt2->next;
      }
    }
    else {
      pt1->power++;
      pt2 = pt1;
      pt1 = pt1->next;
    }
  }

  /*
   * Remove odd powers from even part.
   */
  pt1 = evenodd.p1;
  pt2 = NULL;
  while (pt1) {
    if ((pt1->power&0x1)==1) {
      if (pt2==NULL) {
        evenodd.p1 = pt1->next;
        free(pt1);
        pt1 = evenodd.p1;
      }
      else {
        pt2->next = pt1->next;
        free(pt1);
        pt1 = pt2->next;
      }
    }
    else {
      pt2 = pt1;
      pt1 = pt1->next;
    }
  }

  return evenodd;
}


/*
 * Merge two polyphase components into a polynomial in z:
 *
 *   x(z) = { x_even(z^2) } + z^(-1) { x_odd(z^2) }
 */
zpoly *polyphase_merge(zpoly_pair evenodd)
{
  zpoly *merged;
  zpoly *p;

  merged = zpoly_create(evenodd.p2->power,evenodd.p2->coeff);
  merged = zpoly_add_to(merged,evenodd.p2->next);
  p = merged;
  while (p) {
    p->power--;
    p = p->next;
  }
  merged = zpoly_add_to(merged,evenodd.p1);

  return merged;
}


/*
 * Internal function.
 *                             2
 * Transforms a polynomial in z  to one in z.
 */
zpoly *polyphase_halve_degree(zpoly *poly)
{
  zpoly *p = poly;

  while (p) {
    p->power >>= 1;
    p = p->next;
  }

  return poly;
}


/*
 * Internal function.
 *                                         2
 * Transforms a polynomial in z to one in z .
 */
zpoly *polyphase_double_degree(zpoly *poly)
{
  zpoly *p = poly;

  while (p) {
    p->power <<= 1;
    p = p->next;
  }

  return poly;
}


/*
 * Builds a polyphase matrix:
 *
 *  /                                           \
 *  | row1.p1 = poly1_even  row1.p2 = poly1_odd |
 *  | row2.p1 = poly2_even  row2.p2 = poly2_odd |
 *  \                                           /
 *                                         2
 * Note that the input polynomials are in z  while the output polynomials
 * are in z.
 */
polyphase_matrix polyphase_matrix_build(zpoly_pair evenodd)
{
  polyphase_matrix m;
  m.row1 = polyphase_split(evenodd.p1);
  m.row2 = polyphase_split(evenodd.p2);
  polyphase_halve_degree(m.row1.p1);
  polyphase_halve_degree(m.row1.p2);
  polyphase_halve_degree(m.row2.p1);
  polyphase_halve_degree(m.row2.p2);
  return m;
}


/*
 * Decomposes a polyphase matrix in two polynomials in z.
 * Make sure that the polyphase components are in the right place.
 */
zpoly_pair polyphase_matrix_decompose(polyphase_matrix m)
{
  zpoly_pair evenodd;
  zpoly_pair temp;

  temp.p1 = zpoly_copy(m.row1.p1);
  temp.p2 = zpoly_copy(m.row1.p2);
  temp.p1 = polyphase_double_degree(temp.p1);
  temp.p2 = polyphase_double_degree(temp.p2);
  evenodd.p1 = polyphase_merge(temp);
  zpoly_destroy(temp.p1);
  zpoly_destroy(temp.p2);

  temp.p1 = zpoly_copy(m.row2.p1);
  temp.p2 = zpoly_copy(m.row2.p2);
  temp.p1 = polyphase_double_degree(temp.p1);
  temp.p2 = polyphase_double_degree(temp.p2);
  evenodd.p2 = polyphase_merge(temp);
  zpoly_destroy(temp.p1);
  zpoly_destroy(temp.p2);

  return evenodd;
}


/*
 * Destroy a polyphase matrix;
 */
void polyphase_matrix_destroy(polyphase_matrix m)
{
  zpoly_destroy(m.row1.p1);
  zpoly_destroy(m.row1.p2);
  zpoly_destroy(m.row2.p1);
  zpoly_destroy(m.row2.p2);
}


/*
 * Returns 1 if two polyphase matrices are equal.
 */
int polyphase_matrix_equal(polyphase_matrix m1, polyphase_matrix m2)
{
  if ((zpoly_equal(m1.row1.p1,m2.row1.p1)) &&
      (zpoly_equal(m1.row1.p2,m2.row1.p2)) &&
      (zpoly_equal(m1.row2.p1,m2.row2.p1)) &&
      (zpoly_equal(m1.row2.p2,m2.row2.p2))) {
    return 1;
  }
  else {
    return 0;
  }
}


/*
 * Create a copy of a polyphase matrix;
 */
polyphase_matrix polyphase_matrix_copy(polyphase_matrix m)
{
  polyphase_matrix m_copy;
  m_copy.row1.p1 = zpoly_copy(m.row1.p1);
  m_copy.row1.p2 = zpoly_copy(m.row1.p2);
  m_copy.row2.p1 = zpoly_copy(m.row2.p1);
  m_copy.row2.p2 = zpoly_copy(m.row2.p2);
  return m_copy;
}


/*
 * Transposes a polyphase matrix:
 *
 *  /                                      \
 *  | row1.p1 = row1.p1  row1.p2 = row2.p1 |
 *  | row2.p1 = row1.p2  row2.p2 = row2.p2 |
 *  \                                      /
 */
polyphase_matrix polyphase_matrix_transpose(polyphase_matrix m)
{
  zpoly *temp = m.row1.p2;
  m.row1.p2 = m.row2.p1;
  m.row2.p1 = temp;
  return m;
}


/*
 * Calculate the determinant of a polyphase matrix.
 *
 *     /      \
 * M = | a  b |, determinant = |M| = (ad-bc).
 *     | c  d |
 *     \      /
 *
 * Note that the determinant of a polyphase matrix is always monomial.
 */
zpoly *polyphase_matrix_determinant(polyphase_matrix m)
{
  zpoly *d1 = zpoly_mul(m.row1.p1,m.row2.p2);
  zpoly *d2 = zpoly_mul(m.row2.p1,m.row1.p2);

  d1 = zpoly_sub_from(d1,d2);
  free(d2);

  return d1;
}


/*
 * Calculate the inverse of a polyphase matrix (Cramer's rule).
 *
 *     /      \   -1      1    /     \
 * M = | a  b |, M   = ------- | d -c|
 *     | c  d |        (ad-bc) |-b  a|
 *     \      /                \     /
 */
polyphase_matrix polyphase_matrix_inverse(polyphase_matrix m)
{
  polyphase_matrix minv = {NULL,NULL,NULL,NULL};
  zpoly *det = polyphase_matrix_determinant(m);
  zpoly *temp;

  /*
   * Check if matrix can be inverted (determinant has to be a monomial).
   */
  if (zpoly_degree(det)==0) {

    minv = polyphase_matrix_copy(m);

    /*
     * Swap a & d.
     */
    temp = minv.row1.p1;
    minv.row1.p1 = minv.row2.p2;
    minv.row2.p2 = temp;

    /*
     * Swap b & c.
     */
    temp = minv.row1.p2;
    minv.row1.p2 = minv.row2.p1;
    minv.row2.p1 = temp;

    /*
     * Multiply b & c by -1.
     */
    zpoly_scale(minv.row1.p2,-1.0);
    zpoly_scale(minv.row2.p1,-1.0);

    /*
     * Divide by determinant.
     */
    det->power *= -1;
    det->coeff = 1.0/det->coeff;
    minv = polyphase_matrix_scale(minv,det);

  }

  zpoly_destroy(det);

  return minv;
}


/*
 *  /      \ /      \    /              \
 *  | a  b | | e  f | =  | ae+bg  af+bh |
 *  | c  d | | g  h |    | ce+dg  cf+dh |
 *  \      / \      /    \              /
 */
polyphase_matrix polyphase_matrix_mul(polyphase_matrix m1, polyphase_matrix m2)
{
  polyphase_matrix mnew;
  zpoly *temp;

  mnew.row1.p1 = zpoly_mul(m1.row1.p1,m2.row1.p1);
  temp = zpoly_mul(m1.row1.p2,m2.row2.p1);
  mnew.row1.p1 = zpoly_add_to(mnew.row1.p1,temp);
  zpoly_destroy(temp);

  mnew.row1.p2 = zpoly_mul(m1.row1.p1,m2.row1.p2);
  temp = zpoly_mul(m1.row1.p2,m2.row2.p2);
  mnew.row1.p2 = zpoly_add_to(mnew.row1.p2,temp);
  zpoly_destroy(temp);

  mnew.row2.p1 = zpoly_mul(m1.row2.p1,m2.row1.p1);
  temp = zpoly_mul(m1.row2.p2,m2.row2.p1);
  mnew.row2.p1 = zpoly_add_to(mnew.row2.p1,temp);
  zpoly_destroy(temp);

  mnew.row2.p2 = zpoly_mul(m1.row2.p1,m2.row1.p2);
  temp = zpoly_mul(m1.row2.p2,m2.row2.p2);
  mnew.row2.p2 = zpoly_add_to(mnew.row2.p2,temp);
  zpoly_destroy(temp);

  return mnew;
}


/*
 *   /      \   /          \
 * K*| a  b | = | K*a  K*b |
 *   | c  d |   | K*c  K*d |
 *   \      /   \          /
 */
polyphase_matrix polyphase_matrix_scale(polyphase_matrix m, zpoly *poly)
{
  m.row1.p1 = zpoly_mul_by(m.row1.p1,poly);
  m.row1.p2 = zpoly_mul_by(m.row1.p2,poly);
  m.row2.p1 = zpoly_mul_by(m.row2.p1,poly);
  m.row2.p2 = zpoly_mul_by(m.row2.p2,poly);
  return m;
}


/*
 * Reverse the order of all elements.
 * If reverse_time==1 all powers are multiplied by -1.
 */
polyphase_matrix polyphase_matrix_reverse(polyphase_matrix m, int reverse_time)
{
  polyphase_matrix mnew;

  mnew.row1.p1 = zpoly_reverse(m.row1.p1,reverse_time);
  mnew.row1.p2 = zpoly_reverse(m.row1.p2,reverse_time);
  mnew.row2.p1 = zpoly_reverse(m.row2.p1,reverse_time);
  mnew.row2.p2 = zpoly_reverse(m.row2.p2,reverse_time);

  return mnew;
}


void polyphase_matrix_write(polyphase_matrix m)
{
  redir_printf("(0,0)=[");
  if (zpoly_degree(m.row1.p1)>=0) {
    zpoly_write(m.row1.p1);
  }
  else {
    redir_printf("0");
  }
  redir_printf("], (0,1)=[");
  if (zpoly_degree(m.row1.p2)>=0) {
    zpoly_write(m.row1.p2);
  }
  else {
    redir_printf("0");
  }
  redir_printf("], (1,0)=[");
  if (zpoly_degree(m.row2.p1)>=0) {
    zpoly_write(m.row2.p1);
  }
  else {
    redir_printf("0");
  }
  redir_printf("], (1,1)=[");
  if (zpoly_degree(m.row2.p2)>=0) {
    zpoly_write(m.row2.p2);
  }
  else {
    redir_printf("0");
  }
  redir_printf("]\n");
}