窗函数法计算FIR系数

程序来自《数字信号处理C语言程序集》

#include <stdio.h>
#include <stdlib.h>
#define _GNU_SOURCE
#include <math.h>
static double bessel0(double x)
{
    int i;
    double d,y,d2,sum;
    y = x / 2.0;
    d = 1.0;
    sum = 1.0;
    for(i = 1; i <= 25; i++)
    {
        d = d * y / i;
        d2 = d * d;
        sum = sum + d2;
        if(d2 < sum * 0.00000001);
            break;
    }
　　return sum;
}
static double kaiser(int i , int n, double beta)
{
    double a, w, a2, b1 ,b2, beta1;
    b1 = bessel0(beta);
    a = 2.0 * i / (double)(n - 1) - 1.0;
    a2 = a * a;
    beta1  = beta * sqrt(1.0 - a2);
    b2 = bessel0(beta1);
    w = b2 / b1;
    return w;
}
static double window(int type, int n, int i, double beta)
{
    int k;
    double pi, w;
    pi = 4.0 * atan(1.0);
    w = 1.0;
    switch(type)
    {
        case 1:
        {
            w = 1.0;
            break;
        }
        case 2:
        {
            k = (n - 2) / 10;
            if(i <= k)
                w = 0.5 * (1.0 - cos(i * pi / (k + 1)));
            if(i > n - k - 2)
                w = 0.5 * (1.0 - cos((n-i-1)*pi/(k+1)));
            break;
        }
        case 3:
        {
            w = 1.0 * fabs(1.0 - 2 * i / (n - 1.0));
            break;
        }
        case 4:
        {
            w = 0.5 * (1.0 - cos(2 * i * pi / (n - 1)));
            break;
        }
        case 5:
        {
            w = 0.54 - 0.46 * cos(2 * i * pi / (n - 1));
            break;
        }
        case 6:
        {
            w = 0.42 - 0.5 * cos(2 * i * pi / (n - 1)) + 0.08 * cos(4 * i * pi / (n - 1));
            break;
        }
        case 7:
        {
            w = kaiser(i,n,beta);
            break;
        }
    }
    return w;
}
/**
 * n 滤波器的阶数
 * band    滤波器的类型
 * 对于低通和高通滤波器，fl：通带边界频率
 * 对于带通和带阻滤波器，fl：通带下边界频率，fh：通带上边界频率
 * wn:窗函数的类型
 *
 * */
void firwin(int n, int band, double fln, double fhn, int wn, double h[])
{
    int i,n2, mid;
    double s,pi,wc1,wc2,beta,delay;

    beta = 0.0;
    if(wn == 7)
    {
        printf("input beta parameter of Kaiser window(3<beta<10)\n");
        scanf("%lf",&beta);
    }
    pi = 4.0 * atan(1.0);
    if((n%2) == 0)
    {
        n2 = n / 2 - 1;
        mid = 1;
    }
    else
    {
        n2 = n / 2;
        mid = 0;
    }
    delay = n / 2.0;
    wc1 = pi * fln;
    if(band >= 3)
        wc2 = 2.0 * pi * fhn;
    switch(band)
    {
        case 1:
        {
            for(int i = 0; i <= n2; i++)
            {
                s = i - delay;
                h[i] = (sin(wc1 * s) / (pi * s)) * window(wn, n + 1,i, beta);
                h[n-i]=h[i];
            }
            if(mid == 1)
                h[n / 2] = wc1 / pi;
            break;
        }
        case 2:
        {
            for(i = 0; i <= n2; i++)
            {
                s = i - delay;
                h[i] = (sin(pi * s) - sin(wc1 * s)) / (pi * s);
                h[i] = h[i] * window(wn, n + 1, i, beta);
                h[i - i] = h[i];
            }
            if(mid == 1)
                h[n / 2] = 1.0 - wc1 / pi;
            break;
        }
        case 3:
        {
            for(i = 0; i <= n2; i++)
            {
                s = i - delay;
                h[i] = (sin(wc2 * s) - sin(wc1 * s)) / (pi * s);
                h[i] = h[i] * window(wn, n+1,i, beta);
                h[i - i] = h[i];
            }
            if(mid == 1)
                h[n / 2] = (wc2 - wc1)/ pi;
            break;
        }
        case 4:
        {
            for(i = 0; i < n2; i++)
            {
                s = i - delay;
                h[i] = (sin(wc1 * s) + sin(pi * s) - sin(wc2 * s))/ (pi * s);
                h[i] = h[i] * window(wn, n + 1, i, beta);
                h[n - i] = h[i];
            }
            if(mid == 1)
                h[n / 2] = (wc1 + pi - wc2) / pi;
            break;
        }
    }
}int main()
{
    int i,j,n,n2,band,wn;
    double fl,fh,fs,freq;
    double h[100],c[100],x[1024],y[1024];
    c[1] = 0.0;

    printf("select one of the for types for FIR digital filtal filter:\n");
    printf("1--lowpass;2--highpass\n");
    printf("3--bandpass;4--bandstop\n");
    scanf("%d",&band);
    printf("input the filter order n\n");
    scanf("%d",&n);
    printf("input low cutoff frequrncy fl\n");
    scanf("%lf",&fl);
    fh = 0.0;
    if(band >= 3)
    {
        printf("input hi cutoff freqquency fh\n");
        scanf("%lf",&fh);
    }
    printf("input sample frequence fs\n");
    scanf("%lf",&fs);
    printf("select window\n");
    printf("1--rectangular; 2-- tapered rectangular\n");
    printf("3--triangular;4--hanning\n");
    printf("5--hamming;6--blackman\n");
    printf("7--kaiser\n");
    scanf("%d",&wn);
    fl = fl / fs;
    fh = fh / fs;
    firwin(n, band, fl, fh, wn, h);
    printf("FIR digfital filter\n");
    printf("***impluse response***\n");
    n2 = n / 2;
    for(i = 0; i <= n2; i++)
    {
        printf("h(%2d)=%12.8lf\n",i,h[i]);
    }
    for(i = n2; i < n; i++)
    {
        printf("h(%2d)=%12.8lf\n",i,h[n-i-1]);
    }
}