filter.c

大师写的二代小波经典之作

/*
 *  -*- Mode ANSI C -*-
 *  $Id: filter.c,v 1.4 1996/08/16 17:25:47 fernande Exp $
 *  $Source: /sgi.acct/sweldens/cvs/liftpack/Lifting/filter.c,v $
 *  Author: Gabriel Fernandez
 *
 *  This module finds the filter coefficients used in the wavelet
 *  transform to generate the wavelet coefficients (gamma) using
 *  the contribution of the lambda values. The process uses the
 *  Neville's algorithm (neville.c) for the polynomial interpolation
 *  and a given value, N, to know how many filter coefficients are
 *  needed. The results are returned in a Matrix of N/2+1 rows by
 *  N columns in a linear representation.
 */
/* do not edit anyhting above this line */

/* System header files */
#include <assert.h>
/* FLWT header files */
#include <flwtdef.h>   /* Matrix type and error codes definitions */
#include <filter.h>    /* external declaration of GetFilter() */
#include <mem.h>
#include <neville.h>   /* Neville() function */
#include <util.h>      /* ODD() macro */

/* code */

/*
 * GetFilter function: finds the filter coefficients used in the
 *                     Gamma coefficients prediction routine. The
 *                     Neville polynomial interpolation algorithm
 *                     is used to find all the possible values for
 *                     the filter coefficients. Thanks to this
 *                     process, the boundaries are correctly treated
 *                     without including artifacts.
 *                     Results are ordered in matrix as follows:
 *                         0 in the left and N   in the right
 *                         1 in the left and N-1 in the right
 *                         2 in the left and N-2 in the right
 *                                        .
 *                                        .
 *                                        .
 *                         N/2 in the left and N/2 in the right
 *                     For symmetry, the cases from
 *                         N/2+1 in the left and N/2-1
 *                     to
 *                         N in the left and 0 in the right
 *                     are the same as the ones shown above, but with
 *                     switched sign.
 */
Vector
GetFilter ( const int N )
{
    static Flt xa[40], ya[40], x;
    static int row, col, cc, Nrows;
    Vector filter, ptr;

    assert ( N < sizeof xa / sizeof *xa );   /* N is greater than dim(xa) */

    /* Check for validity of N */
    if ( ODD (N) || (N == 0) )
        return (Vector)NULL;   /* N is not even and greater than zero */

    /* Number of cases for filter calculations */
    Nrows = (N>>1) + 1;    /* N/2 + 1 */

    /* Allocate memory for filter matrix */
    filter = vector( 0, (long)(Nrows*N - 1) );
    ptr = filter;

    /* Generate values of xa */
    xa[0] = (Flt)0.5*(Flt)(1 - N);    /* -N/2 + 0.5 */
    for (col = 1 ; col < N ; col++)
        xa[col] = xa[col-1] + (Flt)1;

    /* Find filter coefficient values */
    filter += ( (Nrows*N) - 1 );   /* go to last position in filter matrix */
    for (row = 0 ; row < Nrows ; row++) {
        x = (Flt)row;
        for (col = 0 ; col < N ; col++) {
            for (cc = 0 ; cc < N ; cc++)
                ya[cc] = (Flt)0;
            ya[col] = (Flt)1;
            *(filter--) = Neville (xa, ya, N, x);
        }
    }

    return ptr;
}