文章

FIR数字滤波器的设计

发表于
作者 lixinghui 9 分钟阅读

6.1 窗函数方法

设N-1阶FIR数字滤波器的单位冲激响应为h(n)，则传递函数H(z)为：

\[H(z)=\sum_{n=0}^{N-1}h(n)z^{-n}\]

窗函数法的设计步骤如下：

1、根据给定的理想频率响应\(H_d(e^{jw})\)，利用傅里叶反变换，求出单位冲激响应\(h_d(n)\)

\[h_d(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}H_d(e^{jw})e^{jwn}dw\]

2、将\(h_d(n)\)乘以窗函数w(n)，得到所要求的FIR滤波器系数h(n)

\[h(n)=w(n)h_d(n),0\leq n \leq N-1\]

常用的窗函数有：

矩形窗：

\[w(n)=1, 0 \leq n \leq N-1\]

三角窗：

\[w(n)=1-|1-\frac{2n}{N-1}|,0 \leq n \leq N-1\]

汉宁(Hanning)窗：

\[w(n)=0.5-0.5cos(\frac{2\pi n}{N-1}),0 \leq n \leq N-1\]

海明(Hamming)窗：

\[w(n)=0.54-0.46cos(\frac{2\pi n}{N-1}),0 \leq n \leq N-1\]

布拉克曼(Blackman)窗：

\[w(n)=0.42-0.5cos(\frac{2\pi n}{N-1})+0.08cos(\frac{4\pi n}{N-1}),0 \leq n \leq N-1\]

凯塞(Kaiser)窗：

\[w(n)=\frac{I_0(\beta \sqrt{1-(1-2n/(N-1))^2})}{I_0(\beta)},0 \leq n \leq N-1\]

其中\(I_0(\beta)\)是第一类零阶修正贝塞耳函数。β是控制窗函数形状的参数，β越大，\(w(n)\)窗越窄，频谱的旁瓣越小，但主孵也相应加宽。β的典型值为4≤β≤9。β=0时，凯塞窗变为矩形窗；β=5.44时，凯塞窗与海明窗接近；β=8.5时，凯塞窗与布拉克曼窗接近。

fir_design.h

1
2
3
4
5
6

#ifndef FIR_DESIGN_H
#define FIR_DESIGN_H

void firwin(int n, int band, double fln, double fhn, int wn, double *h);

#endif

fir_design.c

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "fir_design.h"

static double window(int type, int n, int i, double beta);
static double kaiser(int i, int n, double beta);
static double bessel0(double x);

/**
 * @brief : 用窗函数方法设计线性相位FIR数字滤波器
 *
 * @param n : 整型变量。滤波器的阶数。
 * @param band : 整型变量。滤波器的类型。取值为1、2、3和4，分别对应低通、高通、带通和带阻滤波器。
 * @param fln 对于低通和高通滤波器， fln：通带边界频率；对于带通和带阻滤波器， fln：通带下边界频率,
 * @param fhn 通带上边界频率
 * @param wn : 对于低通和高通滤波器， fln：通带边界频率；对于带通和带阻滤波器， fln：通带下边界频率,fhn:通带上边界频率
 * @param h : 双精度实型一维数组，长度为（n+1）。存放FIR滤波器的系数。
 */
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 = 2.0 * pi * fln;

    if (band >= 3)
    {
        wc2 = 2.0 * pi * fhn;
    }
    switch (band)
    {
    case 1:
    {
        for (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[n - 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[n - 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;
    }
    }
}

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);
}

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 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 * (1.0e-8))
        {
            break;
        }
    }
    return (sum);
}

测试函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

#include "fir_design.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int main(int argc, char **argv)
{
    int i, j, n, n2, band, wn;
    double fl, fh, fs, freq;
    double h[100];

    printf("select one of the four types for FIR digital 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;
    if (band >= 3)
    {
        printf("input high cutoff frequency 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 digital filter \n");
    printf(" * * * * impulse response * * * * \n\n");

    n2 = n / 2;
    for (i = 0; i <= n2; i++)
    {
        j = n - i;
        printf(" h(%2d) = %12.8lf = h(%2d) \n", i, h[i], j);
    }

    return 0;
}

6.2 频域最小误差平方设计

DSP, C语言集
DSP_C语言集
本文由作者按照 CC BY 4.0 进行授权