数字滤波器的时域和频域响应
直接型：
1、直接型结构简单直观，用到的延迟器最少，为（N/M中较大者的个数）
2、系数对滤波器的控制关系不直接，因此调整不方便。
3、具体实现滤波器时，系数的量化误差将使滤波器的频响产生很大的改变，甚至影响到系统的稳定性。
4、直接型结构一般用以实现低阶系统，对高阶系统，级联型和并联型结构更有优越性。
级联型：
1、级联型结构便于准确的实现系统的零、极点，也便于进行滤波器性能的调整，从整体看，容易控制零点也难于调整极点
2、级联型结构的零极点陪对方式和基本节级联次序具有很大的灵活性，但由于有限字长的影响对于不同的排列，运算误差各有不同。
3、每个基本节都有相同的结构，硬件实现时可以用一个二阶节进行时分复用，故只需要很少的顺序存储单元和运算部件。
并联型：
1、运算速度快，各基本节的误差无不影响，总误差低于级联型结构
2、从整体上看，并联结构容易调整极点位置，不容易直接控制零点。
4.1 数字滤波器的频率响应
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
#include <math.h>
/**
* @brief 计算数字滤波器的频率响应、幅频响应和相频响应
*
* @param b ：长度为（m+1），存放滤波器分子的系数b
* @param a ：长度为（n+1），存放滤波器分母的系数a
* @param m ：滤波器分子多项式的阶数
* @param n ：滤波器分母多项式的阶数
* @param x ：长度为len，sign=0时存放滤波器频率响应的实部，
* sign=1时存放滤波器幅频响应H(w),sign=2时存放分贝表示的滤波器幅频响应H(w)
* @param y ：长度为len，sign=0时存放滤波器频率响应的虚部，
* sign=1,2的时候存放滤波器的相频响应
* @param len ：滤波器响应的长度
* @param sign ：sign=0时为滤波器频率响应的实部和虚部，sign=1时输出幅频响应和相频响应
* sign=2时输出dB形式的响应
*/
void gain(double *b, double *a, int m, int n, double *x, double *y, int len, int sign)
{
int i, k;
double ar, ai, br, bi, zr, zi, im, re, den, numi, numr, freq, temp;
for (k = 0; k < len; k++)
{
freq = k * 0.5 / (len - 1);
zr = cos(-8.0 * atan(1.0) * freq);
zi = sin(-8.0 * atan(1.0) * freq);
br = 0.0;
bi = 0.0;
for (i = m; i > 0; i--)
{
re = br;
im = bi;
br = (re + b[i]) * zr - im * zi;
bi = (re + b[i]) * zi + im * zr;
}
ar = 0.0;
ai = 0.0;
for (i = n; i > 0; i--)
{
re = ar;
im = ai;
ar = (re + a[i]) * zr - im * zi;
ai = (re + a[i]) * zi + im * zr;
}
br = br + b[0];
ar = ar + 1.0;
numr = ar * br + ai * bi;
numi = ar * bi - ai * br;
den = ar * ar + ai * ai;
x[k] = numr / den;
y[k] = numi / den;
switch (sign)
{
case 1:
temp = sqrt(x[k] * x[k] + y[k] * y[k]);
y[k] = atan2(y[k], x[k]);
x[k] = temp;
break;
case 2:
temp = x[k] * x[k] + y[k] * y[k];
y[k] = atan2(y[k], x[k]);
x[k] = 10.0 * log10(temp); // 20.0*log10(temp)?
break;
default:
break;
}
}
}
测试函数：
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
#include <stdio.h>
#include "filter_process.h"
int main(int argc, char **argv)
{
int i;
double a[]={1,-0.184213803077536,0.177578118927698};
double b[]={0.248341078962541,0.496682157925081,0.248341078962541};
double f,x[300],y[300];
FILE *fp=NULL;
gain(b,a,2,2,x,y,300,2);
if((fp=fopen("gainam.dat","w"))==NULL)
{
printf("Error opening gainam file\n");
return 1;
}
for(i=0;i<300;i++)
{
f=i*0.5/299;
fprintf(fp,"%lf %lf\n",f,x[i]);
}
fclose(fp);
if((fp=fopen("gainph.dat","w"))==NULL)
{
printf("Error opening gainph file\n");
return 1;
}
for(i=0;i<300;i++)
{
f=i*0.5/299;
fprintf(fp,"%lf %lf\n",f,y[i]);
}
fclose(fp);
printf("done\n");
}
4.2 级联型数字滤波器的频率响应
级联型数字滤波器可以使用卷积，将多级滤波器生成一个高阶滤波器，再使用上面的gain函数就可以计算了。
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
/**
* @brief 计算级联型数字滤波器的频率响应、幅频响应和相频响应
*
* @param b ：长度为ns*(n+1)，存放滤波器分子的系数b[j][i]，j为阶数i为每阶系数
* @param a ：长度为ns*(n+1)，存放滤波器分母的系数a[j][i]
* @param n ：级联型滤波器每节的阶数
* @param ns ：级联型滤波器的n阶节数L
* @param x ：长度为len，sign=0时存放滤波器频率响应的实部，
* sign=1时存放滤波器幅频响应H(w),sign=2时存放分贝表示的滤波器幅频响应H(w)
* @param y ：长度为len，sign=0时存放滤波器频率响应的虚部，
* sign=1,2的时候存放滤波器的相频响应
* @param len ：滤波器响应的长度
* @param sign ：sign=0时为滤波器频率响应的实部和虚部，sign=1时输出幅频响应和相频响应
* sign=2时输出dB形式的响应
*/
void gainc(double *b, double *a, int n, int ns, double *x, double *y, int len, int sign)
{
int i, j, k, n1;
double ar, ai, br, bi, zr, zi, im, re, den, numi, numr, freq, temp;
double hr, hi, tr, ti;
n1 = n + 1;
for (k = 0; k < len; k++)
{
freq = k * 0.5 / (len - 1);
zr = cos(-8.0 * atan(1.0) * freq);
zi = sin(-8.0 * atan(1.0) * freq);
x[k] = 1.0;
y[k] = 0.0;
for (j = 0; j < ns; j++)
{
br = 0.0;
bi = 0.0;
for (i = n; i > 0; i--)
{
re = br;
im = bi;
br = (re + b[j * n1 + i]) * zr - im * zi;
bi = (re + b[j * n1 + i]) * zi + im * zr;
}
ar = 0.0;
ai = 0.0;
for (i = n; i > 0; i--)
{
re = ar;
im = ai;
ar = (re + a[j * n1 + i]) * zr - im * zi;
ai = (re + a[j * n1 + i]) * zi + im * zr;
}
br = br + b[j * n1 + 0];
ar = ar + 1.0;
numr = ar * br + ai * bi;
numi = ar * bi - ai * br;
den = ar * ar + ai * ai;
hr = numr / den;
hi = numi / den;
tr = x[k] * hr - y[k] * hi;
ti = x[k] * hi + y[k] * hr;
x[k] = tr;
y[k] = ti;
}
switch (sign)
{
case 1:
temp = sqrt(x[k] * x[k] + y[k] * y[k]);
if (temp != 0.0)
{
y[k] = atan2(y[k], x[k]);
}
else
{
y[k] = 0.0;
}
x[k] = temp;
break;
case 2:
temp = x[k] * x[k] + y[k] * y[k];
if (temp != 0.0)
{
y[k] = atan2(y[k], x[k]);
}
else
{
temp = 1e-40;
y[k] = 0.0;
}
x[k] = 10.0 * log10(temp); // 20.0*log10(temp)?
break;
}
}
}
测试函数
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
int main(int argc, char **argv)
{
int i,n,ns;
// a[0][3]={1.0f,-0.8f,0.0f};
// a[1][3]={1.0f,-0.95f,0.9f};
// b[0][3]={0.2f,0.0f,0.0f};
// b[1][3]={0.05f,-0.1f,0.0f};
double a[2][3]={0};
double b[2][3]={0};
double f,x[300],y[300];
FILE *fp=NULL;
n=2,ns=2;
gainc((double *)b,(double *)a,n,ns,x,y,300,1);
if((fp=fopen("gainam.dat","w"))==NULL)
{
printf("Error opening gainam file\n");
return 1;
}
for(i=0;i<300;i++)
{
f=i*0.5/299;
fprintf(fp,"%lf %lf\n",f,x[i]);
}
fclose(fp);
if((fp=fopen("gainph.dat","w"))==NULL)
{
printf("Error opening gainph file\n");
return 1;
}
for(i=0;i<300;i++)
{
f=i*0.5/299;
fprintf(fp,"%lf %lf\n",f,y[i]);
}
fclose(fp);
printf("done\n");
}
本文由作者按照 CC BY 4.0 进行授权