main.cpp

我自己改写的Remez,FIR滤波器参数计算程序

//    Remez Exchange FIR Filter Design Program
//    Copyright (C) 1999  Jim A. Ledin ([e-mail])
//
//    This program is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation; either version 2 of the License, or
//    (at your option) any later version.
//
//    This program is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with this program; if not, write to the Free Software
//    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include "remez.h"
#include "remez.tab.h"

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

#include "types.h"

// Global variables filled in by yyparse()
BandEdgeVec BandEdge;
AmplitudeVec Amplitude;
WeightVec Weight;
int FilterType = BANDPASS;
int Order = -1;

int main(int argc, char *argv[])
{
    FILE* flptr;
    extern int ErrorCount, LineNumber;
    extern int yyparse();

    if (argc < 2)
    {
        fprintf(stderr,
            "usage: %s filename.rem [freqresp.csv] [filename.m]\n"
            "       filename.rem  = Input filter specifications\n"
            "       freqresp.csv  = Spreadsheet format frequency response data\n"
            "       filename.m    = Filter parameters in Matlab format\n"
            "       C source output file = filter.dat\n",
            argv[0]);
        return 0;
    }

    if ((flptr = freopen(argv[1], "r", stdin)) == NULL) {
        char buf[200];
        sprintf(buf, "ERROR: Can't open %s for input\n", argv[1]);
        fprintf(stderr, buf);
        return 0;
    }

    yyparse();

    if (ErrorCount > 0) {
        char buf[200];
        sprintf(buf, "%d error%s detected\n", ErrorCount,
            ErrorCount == 1? "": "s");
        fprintf(stderr, buf);
        return 0;
    }

    // Check input:
    //  1) Length of BandEdges and Amplitude must be equal
    int bes = BandEdge.size();
    int ws = Weight.size(), as = Amplitude.size();

    if (bes != 2*as)
    {
        fprintf(stderr, 
            "Number of amplitudes must match number of band edge pairs\n");
        return 0;
    }

    //  2) Length of Weight must be zero or match BandEdge
    if (ws == 0)
    {
        for (int i=0; i<bes/2; i++)
            Weight.push_back(1.0);  // Fill in default weights
    }
    else
    {
        if (bes != ws*2)
        {
            fprintf(stderr, 
                "Number of weights must match number of band edge pairs\n");
            return 0;
        }
    }

    //  3) Each band must have the lower frequency listed first
    for (int i=0; i<bes-1; i++)
        if (BandEdge[i] >= BandEdge[i+1])
        {
            fprintf(stderr, 
                "Band edge frequencies must appear in ascending order\n");
            return 0;
        }

    //  4) Lowest band must start at 0
    if (BandEdge[0] != 0)
    {
        fprintf(stderr, 
            "First band edge must be at frequency = 0\n");
        return 0;
    }

    //  5) Highest band must end at 1
    if (fabs(BandEdge[bes-1] - 1) > 1e-6)
    {
        fprintf(stderr, 
            "Last band edge must be at frequency = 1\n");
        return 0;
    }

    //  6) Weights must be > 0
    for (i=0; i<bes/2; i++)
        if (Weight[i] <= 0)
        {
            fprintf(stderr, 
                "Weights must all be > 0\n");
            return 0;
        }

    //  7) Order must be specified
    if (Order < 0)
    {
        fprintf(stderr, 
            "Order must be specified\n");
        return 0;
    }

    int numtaps = Order + 1;
    int numband = bes/2;

    double* h       = new double[numtaps];
    double* bands   = new double[bes];
    double* des     = new double[bes];
    double* weight  = new double[bes/2];

    // Set up input arrays for remez function
    for (i=0; i<bes; i++)
    {
        bands[i] = BandEdge[i] * 0.5;
        des[i] = Amplitude[i];

        if (i < bes/2)
            weight[i] = Weight[i];
    }

    // Prepare the output file
    FILE* c_iov = fopen("filter.dat", "w");
    if (!c_iov)
    {
        fprintf(stderr, 
            "Unable to open filter.dat for writing\n");
        return 0;
    }

    // If specified, set up for the frequency response output file
    char* freqresp_file = 0;
    if (argc > 2)
        freqresp_file = argv[2];

    // If specified, open the Matlab output file
    FILE* matlab_iov = 0;
    if (argc > 3)
    {
        matlab_iov = fopen(argv[3], "w");

        if (!matlab_iov)
        {
            fprintf(stderr, 
                "Unable to open %s for writing\n", argv[4]);
            return 0;
        }
    }

    // Design the filter
    int converge_failed = remez(h, numtaps,
           numband, bands, des, weight, FilterType, freqresp_file);

    if (converge_failed)
    {
        fprintf(stderr, 
            "Design failed to converge\n");
        fclose(c_iov);
        if (matlab_iov)
            fclose(matlab_iov);

        return 0;
    }

    // Write the C output file
    fprintf(c_iov, "#define N %d\n\n", numtaps);
    fprintf(c_iov, "static double H[N] =\n{\n");
    for (i=0; i<numtaps; i++)
        fprintf(c_iov, "%30.18le%s\n", h[i], i<numtaps-1 ? "," : "");

    fprintf(c_iov, "};\n");
    fclose(c_iov);

    // If specified, write the Matlab file
    if (matlab_iov)
    {
        fprintf(matlab_iov, "N = %d\n", Order);

        fprintf(matlab_iov, "F = [ ...\n");
        for (i=0; i<bes/2; i++)
            fprintf(matlab_iov, "    %.20le %.20le ...\n",
            BandEdge[2*i], BandEdge[2*i+1]);
        fprintf(matlab_iov, "    ];\n\n");

        fprintf(matlab_iov, "A = [ ...\n");
        for (i=0; i<bes/2; i++)
            fprintf(matlab_iov, "    %.20le %.20le ...\n", 
            Amplitude[i], Amplitude[i]);
        fprintf(matlab_iov, "    ];\n\n");

        fprintf(matlab_iov, "W = [ ...\n");
        for (i=0; i<numband; i++)
            fprintf(matlab_iov, "    %.20le ...\n", Weight[i]);
        fprintf(matlab_iov, "    ];\n\n");

        fprintf(matlab_iov, "H = [ ...\n");
        for (i=0; i<numtaps; i++)
            fprintf(matlab_iov, "    %.20le ...\n", h[i]);
        fprintf(matlab_iov, "    ];\n\n");
        fclose(matlab_iov);
    }

    // Clean up
    delete[] h;
    delete[] bands;
    delete[] des;
    delete[] weight;

    return 0;
}