filtertest.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
 
#include <sp/spLib.h>
 
#define SIG_LENGTH 16
#define FIR1_ORDER 5
#define FIR1_LP_CUTOFF 0.3
#define FIR1_HP_CUTOFF 0.6
 
int main(int argc, char **argv)
{
    long order;
    double cutoff;
    double beta;
    double gain;
    spDVector x, y;
    spDVector h;
    spDVector poles;
    spDVectors coefs;
    spDVector t;
    spDVector theta1, theta2;
    spDVector x1, x2;
    spDVector y1;
    unsigned long options;
 
    x = xdvzeros(SIG_LENGTH);
    x->data[0] = 1.0;
 
    h = xdvfir1(FIR1_ORDER, FIR1_LP_CUTOFF);
    dvdump(h);
    xdvfree(h);
 
    h = xdvfir1bp(FIR1_ORDER, FIR1_LP_CUTOFF, FIR1_HP_CUTOFF);
    dvdump(h);
    xdvfree(h);
    
    h = xdvfir1bs(FIR1_ORDER, FIR1_LP_CUTOFF, FIR1_HP_CUTOFF);
    dvdump(h);
    xdvfree(h);
    
    h = xdvfir1(FIR1_ORDER + 1, FIR1_HP_CUTOFF);
    dvdump(h);
    xdvfree(h);
    
    h = xdvfir1hp(FIR1_ORDER, FIR1_HP_CUTOFF);
    dvdump(h);
    xdvfree(h);
    
    /* check for butterworth filter */
    poles = xdvbuttap(3);
    dvdump(poles);
    xdvfree(poles);
 
    coefs = xdvsbutter(3, 0.25, SP_FALSE);
    dvdump(coefs->vector[0]);
    dvdump(coefs->vector[1]);
    xdvsfree(coefs);
    
    coefs = xdvsbutter(3, 0.25, SP_TRUE);
    dvdump(coefs->vector[0]);
    dvdump(coefs->vector[1]);
    xdvsfree(coefs);
 
    /* check for Chebyshev type I filter */
    poles = xdvcheb1ap(5, 3.0, &gain);
    dvdump(poles);
    fprintf(stderr, "gain = %f\n\n", gain);
    xdvfree(poles);
 
    coefs = xdvscheby1(5, 3.0, 0.5, SP_FALSE);
    dvdump(coefs->vector[0]);
    dvdump(coefs->vector[1]);
    y = xdvfilter(coefs->vector[0], coefs->vector[1], x);
    dvdump(y);
    xdvfree(y);
    xdvsfree(coefs);
    
    coefs = xdvscheby1(5, 3.0, 0.5, SP_TRUE);
    dvdump(coefs->vector[0]);
    dvdump(coefs->vector[1]);
    y = xdvfilter(coefs->vector[0], coefs->vector[1], x);
    dvdump(y);
    xdvfree(y);
    xdvsfree(coefs);
    
    kaiserord(1000.0 / 4000.0, 1500.0 / 4000.0, 0.05, 0.01, SP_FALSE, &order, &cutoff, &beta);
    printf("kaiserord result: order = %ld, cutoff = %f, beta = %f\n", order, cutoff, beta);
 
    h = xdvkaiser(order + 1, beta);
    dvdump(h);
    xdvfree(h);
 
    xdvfree(x);
 
    t = xdvinit(0.0, 0.00025, 1.0);
    theta1 = xdvscoper(t, "*", 2.0 * PI * 30.0);
    theta2 = xdvscoper(t, "*", 2.0 * PI * 60.0);
    x1 = xdvsin(theta1);
    x2 = xdvsin(theta2);
    x = xdvoper(x1, "+", x2);
 
#if 1
    options = SP_DECIMATE_OPTION_KEEP_EDGE;
#else
    options = SP_DECIMATE_OPTION_USE_FIR/* | SP_DECIMATE_OPTION_KEEP_EDGE*/;
#endif
    y = xdvdecimateex(x, 4, 0, options);
    printf("-------- decimate result --------\n");
    y1 = NODATA;
    {
        long m;
        long framel;
        long max_output_buf_length;
        long delay;
        long nprocess;
        long ypos;
        spDecimateRec decimate;
 
        framel = 128;
        
        if ((decimate = spDecimateOpen(4, 0, framel, options, &max_output_buf_length, &delay)) != NULL) {
            spDebug(10, "main", "max_output_buf_length = %ld, delay = %ld\n", max_output_buf_length, delay);
            y1 = xdvzeros(y->length);
            ypos = 0;
            for (m = 0;; m++) {
                nprocess = spDecimateProcess(decimate, &x->data[m * framel], framel, &y1->data[ypos],
                                             SP_FALSE, &delay);
                spDebug(100, "main", "nprocess = %ld, delay = %ld\n", nprocess, delay);
                ypos += nprocess;
                
                if ((m + 1) * framel + framel >= x->length) {
                    spDebug(10, "main", "last frame: ypos = %ld, y1->length = %ld\n", ypos, y1->length);
                    nprocess = spDecimateProcess(decimate, &x->data[(m + 1) * framel], x->length - (m + 1) * framel, &y1->data[ypos],
                                                 SP_TRUE, &delay);
                    ypos += nprocess;
                    spDebug(100, "main", "last frame: nprocess = %ld, delay = %ld, ypos = %ld / %ld\n",
                            nprocess, delay, ypos, y1->length);
                    break;
                }
            }
 
            spDecimateClose(decimate);
        }
    }
    if (y1 != NODATA) {
        long k;
        spDVector diffv;
        double dist;
        diffv = xdvoper(y, "-", y1);
        dvsquare(diffv);
        for (k = 0; k < diffv->length; k++) {
            if (diffv->data[k] > 0.0001) {
                printf("%ld: big distortion: %f\n", k, diffv->data[k]);
            }
        }
        dist = dvmean(diffv);
        printf("distortion between y (length: %ld) and y1 (length: %ld) = %f\n", y->length, y1->length, dist);
        
        /*dvnfdump(stdout, y, y1, NULL);*/
        xdvfree(y1);
    } else {
        dvdump(y);
    }
    xdvfree(y);
    
    xdvfree(x);
    xdvfree(x1);
    xdvfree(x2);
    xdvfree(theta1);
    xdvfree(theta2);
    xdvfree(t);
    
    return 0;
}