Libfilth(一个滤波器C库)使用

Libfilth使用说明

winshton

2009年2月

（*本文大部分翻译自libfilth,还有一部分是个人使用实践

*时间水平均有限，翻译的不完整，尤其第二章可以忽略）

版本历史修改记录

版本	作者	日期	备注
V1.0	winshton	2009-2-1	创建

目录

版本历史修改记录    1

1. 概述    5

2. 库文件分析    5

2.1. filth.h/filth.c    6

2.1.1. quantize()    6

2.1.2. filter_strerror()    6

2.2. filter.h/filter.c    6

2.2.1. fir()    6

2.2.2. pfir()    7

2.2.3. updatepq()    7

2.3. window.h/window.c    8

2.3.1. wd_boxcar()    8

2.3.2. wd_triang()    8

2.3.3. wd_hanning()    8

2.3.4. wd_hamming()    9

2.3.5. wd_blackman()    9

2.3.6. wd_flattop()    9

2.3.7. wd_kaiser()    9

2.4. filtanalysis.c    10

2.4.1. fa_ampfunc()    10

2.4.2. fa_response()    10

2.4.3. fa_cerror()    11

2.4.4. ia_response()    11

2.4.5. ia_sresp()    12

2.4.6. ia_response_allp()    12

2.5. firwindow.c    13

2.5.1. fd_window()    13

2.6. firls.c    13

2.6.1. fd_ls()    13

2.6.2. fd_cls()    14

2.6.3. fd_clsgd()    14

2.7. firremez.c    15

2.7.1. fd_remez()    15

2.8. firminphase.c    16

2.8.1. fd_minphase()    16

2.8.2. fd_minphaseminimax()    16

2.8.3. fd_minphasels()    17

2.8.4. fd_cepminphase()    18

2.9. firminimax.c    19

2.9.1. fd_minimax()    19

2.9.2. fd_cminimax()    20

2.9.3. fd_cminimaxgd()    20

2.10. firoptwin.c    21

2.10.1. fd_qpoptwin()    21

2.10.2. fd_minimaxoptwin()    22

2.10.3. fa_erroroptwin()    23

2.11. fireq.c    23

2.11.1. fd_cminimaxeq()    23

2.11.2. fd_cminimaxgdeq()    24

2.12. firfb.c    25

2.12.1. fd_halfband()    25

2.12.2. fd_dftfbadv()    26

2.12.3. fd_dftfb    26

2.12.4. ff_dftfb()    27

2.12.5. fd_pdftfb()    27

2.12.6. ff_pdftfb()    27

2.13. firtrans.c    28

2.13.1. ft_parallel()    28

2.13.2. ft_quantize()    28

2.14. iiranalog.c    29

2.14.1. id_bessel()    29

2.14.2. id_butterworth()    30

2.14.3. id_chebyshev()    30

2.14.4. id_ichebyshev()    31

2.14.5. id_elliptic()    32

2.15. iirtrans.c    33

2.15.1. it_lp2lp()    33

2.15.2. it_lp2hp()    34

2.15.3. it_lp2bp()    34

2.15.4. it_lp2bs()    35

2.15.5. it_bilinear()    35

2.15.6. it_bilinear_adv()    36

2.15.7. it_sos2poly()    37

2.15.8. it_sos2allpass()    38

2.15.9. it_quant_allpass()    38

2.15.10. it_allpass2poly()    38

2.16. matrix.h/matrix.c    39

2.16.1. mat_calloc()    39

2.16.2. mat_free()    39

2.16.3. mat_conv()    39

2.16.4. mat_cconv()    40

2.16.5. mat_flip()    40

2.16.6. mat_cflip()    40

2.16.7. mat_mult()    41

2.16.8. mat_cmult()    41

2.16.9. mat_vmult()    41

2.16.10. mat_cvmult()    42

2.16.11. mat_add()    42

2.16.12. mat_cadd()    42

2.16.13. mat_sub()    42

2.16.14. mat_csub()    43

2.16.15. mat_omult()    43

2.16.16. mat_comult()    43

2.16.17. mat_odiv()    44

2.16.18. mat_codiv()    44

2.16.19. mat_real()    44

2.16.20. mat_imag()    45

2.16.21. mat_conj()    45

2.16.22. mat_inv()    45

2.16.23. mat_cinv()    45

2.16.24. mat_eye()    46

2.16.25. mat_ceye()    46

2.16.26. mat_solve()    46

2.16.27. mat_csolve()    46

2.16.28. mat_linspace()    47

2.16.29. mat_clinspace()    47

2.16.30. mat_clinspace_pol()    47

3. Libfilth的windows移植    48

4. libfilth的测试和使用分析    49

4.1. 设计IIR滤波器    50

4.1.1. 处理IIR滤波器设计的频率响应    52

4.1.2. 处理IIR滤波器设计的阶数    54

4.1.3. IIR滤波器实现    56

4.1.4. 模拟滤波器变换    58

4.1.5. IIR设计实现总结    62

4.2. 使用Remes算法设计线性FIR滤波器    62

4.2.1. 使用Remes变换算法求解多项式    63

4.2.2. 利用Remes算法计算FIR滤波器    64

4.2.3. 利用Remes算法计算FIR低通滤波器    66

4.2.4. 利用Remes算法计算FIR高通滤波器    67

4.3. 数字全通自适应递归滤波器    68
概述

Libfilth 是一套设计、分析和应用数字和模拟滤波器的程序库。它包含一些基本的滤波器功能。Libfilth为滤波器的设计、分析和变换提供以下类型：

带有线性相位和最小二乘法的FIR滤波器。
大有复杂规范和最小二乘法的FIR滤波器。
使用线性编程提供线性相位和上下限设计的FIR滤波器。
使用Remes算法提供线性相位和上下限设计的FIR滤波器。
提供最小相位和最佳幅值的FIR滤波器。
提供复杂规范和上下限设计的FIR滤波器。
提供组延迟限制的FIR滤波器。
二次编程实现最优窗的FIR滤波器。
提供最小相位谱的FIR滤波器。
使用倒谱设计的FIR滤波器。
使用模拟滤波器设计的贝塞尔-汤姆逊，巴特沃斯，切比雪夫I型， II和椭圆函数滤波器。
模拟-模拟与模拟-数字滤波器的转换。
实现全通IIR滤波器。
FIR 半带滤波器
DFT的滤波器组和平行的DFT滤波器组的设计和实施。
FIR、IIR和模拟滤波器的性能计算分析。
窗函数功能。

带有fir和iir前缀的文件完成FIR和IIR滤波器的设计和转换功能。Filtannalysis.c完成所有滤波器的分析功能，window.c和window.h完成窗函数的计算功能。

AIP的函数前缀使用如下约定：

fd_ FIR Design
fe_ FIR Execute
fa_ FIR Analysis
ft_ FIR Transform
ff_ FIR Free (memory)
id_ IIR Design
ia_ IIR Analysis
it_ IIR Transform
wd_ Calculate window function

库文件分析
filth.h/filth.c

定义库中使用的所有全局数据类型和错误管理函数。

两个类型定义：

typedef double _ftype_t

设计函数中使用的双精度浮点指针类型。

typedef float _fshort_t

执行函数中使用的单精度浮点指针类型。

有两个函数：
quantize()

将输入值转换为小数表示法。

参数:

x     Floating point value

n     Number of bits [1,64]

nx     Numeartor

dx     Eponent for denominator [-128,128]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.
filter_strerror()

返回字符串描述的错误代码，拓展strerror()。

参数:

errnum 错误号

返回:

字符串描述的错误代码。
filter.h/filter.c

设计，实施，转化和分析不同类型的模拟和数字滤波器。这两个文件是滤波器的总纲，滤波器需要的所有类型定义和宏定义等基本资源都从这两个文件获得，而且滤波器的库接口头文件就是filter.h，所有相关的函数声明都在这里了。

其中包括若干数据结构定义、宏定义和类型定义。

此外在filter.c中包含三个基本的处理函数：
fir()

被pfir调用，内联函数。

C implementation of FIR filter , *表示卷积。

本函数似乎是卷积运算函数，输入数组进行一次求积并累加，但这种计算又不像是卷据运算。

参数:

n

过滤阀数，where mod

w

过滤阀

x

输入信号必须是一个顺序索引的环形缓冲区

返回:

卷积结果。

pfir()

C implementation of parallel FIR filter , *表示卷积。

本函数似乎是采集信号值的卷积运算，但看不太懂。

参数:

n	过滤阀数 where mod
d	滤波器数
xi	唤醒缓冲区当前索引
w	过滤阀 [k,n]
x	输入信号必须是一个顺序索引的环形缓冲区
y	输出缓冲区
s	输出缓冲步幅

返回:

输出缓冲y的指针。

updatepq()

内联函数。

添加数据到循环队列，用与并行FIR滤波器设计。

Add new data to circular queue designed to be used with a parallel FIR filter.

参数:

n

过滤阀数 where mod

d

滤波器数

xi

唤醒缓冲区当前索引

xq

环形缓冲n*2,k]

in

输入新数据*s]

s

输入冲步幅

返回:

Xq的新索引。
window.h/window.c

该文件提供计算窗函数。窗函数包括：Boxcar, Triang, Hanning, Hamming, Blackman, Flattop 和Kaiser.
wd_boxcar()

截断窗，又叫Rectangular窗，效果如同没有加窗一样，它的作用只是将信号截短。其谱泄露最大。该窗可以用来分析持续时间比窗短的信号。

Boxcar

参数:

n

窗长

w

窗参数缓冲区 [n]
wd_triang()

Triang a.k.a Bartlett

参数:

n

Window length

w

Buffer for window parameters [n]
wd_hanning()

Hanning

参数:

n

Window length

w

Buffer for window parameters [n]
wd_hamming()

Hamming

参数:

n

Window length

w

Buffer for window parameters [n]
wd_blackman()

Blackman

参数:

n

Window length

w

Buffer for window parameters [n]
wd_flattop()

Flattop

参数:

n

Window length

w

Buffer for window parameters [n]
wd_kaiser()

Kaiser

The beta parameter trades the rejection of the low pass filter against the transition width from passband to stop band. Larger Beta means a slower transition and greater stop band rejection. See Rabiner and Gold (Theory and Application of DSP) under Kaiser windows for more about Beta. The following table from Rabiner and Gold gives some feel for the effect of Beta:

All the passband ripple and the stop-band ripple is in dB, width of transition band .

参数:

n

Window length

w

Buffer for window parameters [n]

b

Beta parameter for Kaiser window, Beta >= 1
filtanalysis.c

该文件函数功能用于分析模拟和数值滤波器的特点。
fa_ampfunc()

为对称型和反对称型1-4FIR滤波器计算振幅，该函数的速度比fa_response()快。

参数:

n

过滤阀数

w

滤波器数

k

频率点数

f

频率点 [0,0.5] [k]

a

得出的振幅 [k]

flags

类型标志, see filter.h

返回:

成功返回0 ， -1返回的错误，and errno is set appropriately.
fa_response()

为任何实型滤波器计算响应(量值，功率，相位，群延迟)。

Calculate filter responses (magnitude, power, phase, group delay) for any type of real fir filter.

参数:

n

过滤阀数

w

滤波器数

k

频率点数

f

频率点矢量 [0,0.5] [k]

r

返回结果 [k]

flags

类型标志, see filter.h

返回:

成功返回0 ， -1返回的错误，and errno is set appropriately.

fa_cerror()

为所有实型FIR滤波器计算权重误差。

Calculate the weighted error function for any type of real FIR filter.

参数:

n	过滤阀数
w	滤波器数
k	频率点数
f	频率点矢量 [0,0.5] [k]
hd	理想的频率响应 used during the design of w [k]
v	权重 vector (if no weighting desired set to NULL) [k]
r	返回结果 [k]
flags	类型标志, see filter.h

返回:

成功返回0 ， -1返回的错误，and errno is set appropriately.

ia_response()

IIR滤波器分析函数，为由b/a给出的实型IIR滤波器计算频率响应、量值、功率、相位响应或群延迟。

Analysis function for IIR filters, calculates frequency response, magnitude, power, phase response or group delay for the real IIR filter given by b/a.

参数:

n	滤波器多项式分子数
b	滤波器多项式分子数组[n]
m	滤波器多项式分母数
a	滤波器多项式分母数组[m]
k	频率响应输出点数
f	频率采样点(范围[0,0.5])数组 [k]
r	返回频率响应结果 [k]
flags	分析功能类型标志, see filter.h

返回:

成功返回0 ， -1返回的错误，and errno is set appropriately.

ia_sresp()

为模拟滤波器滤波器分析函数，为由num/den给出的实型模拟滤波器计算频率响应、量值、功率、相位响应或群延迟。

参数:

n	Filter order
num	Numerator filter taps [3n/2] if n is even [3n/2-1] if n is odd
den	Denominator filter taps [3n/2] if n is even [3n/2-1] if n is odd
g	滤波器增益
k	频率点数
f	频率点矢量 [0,0.5] [k]
r	返回结果 [k]
flags	类型标志, see filter.h

返回:

成功返回0 ， -1返回的错误，and errno is set appropriately.

ia_response_allp()

作为两段全通带IIR滤波器的分析函数。为实型IIR滤波器计算频率响应、量值、功率、相位响应或群延迟。

Analysis function for IIR filters implemented as sum of two all-pass sections. The function calculates frequency response, magnitude, power or phase response for the real IIR filter.

参数:

a1	Linked list of all-pass sections number 1
a2	Linked list of all-pass sections number 2
k	频率点数
f	频率点矢量 [0,0.5] [k]
r	返回结果 [k]
flags	类型标志, see filter.h

返回:

成功返回0 ， -1返回的错误，and errno is set appropriately.

firwindow.c

该文件中的函数在为设计有限脉冲响应滤波器提供窗功能。

fd_window()

设计FIR滤波器加窗，该函数会使用window.c文件中定义的各种窗函数。

参数:

n	Number of filter taps, must be odd for high-pass and band-stop filters
w	Buffer for the filter taps [n]
fc	Cutoff frequencies [0,0.5] ([1] for low-pass and high-pass, [2] for band-pass and band-stop)
flags	Window and filter type as defined in filter.h
opt	Beta constant used only when designing using Kaiser windows

返回:

成功返回0 ， -1返回的错误，and errno is set appropriately.

firls.c

使用最小二乘法和复/实规格设计有限脉冲响应滤波器。

fd_ls()

使用最小二乘频率抽象法实现线性FIR滤波器。

该滤波器用于设计解决如下方程。

其中W是一个加权对角线矩阵的对角线V，从频率矢量 f中生成，ad是理想的幅度响应，a0是滤波器W的一半。

where W is a diagonal weighting matrix with the diagonal v, is generated from the frequency vector f, ad is the desired amplitude response and a0 is half the filter w.

参数:

n	Filter length must be odd for HP and BS filters
w	Buffer for the filter taps [n]
k	Number of frequency sample points
f	Frequency vector [k] [0,0.5]
ad	Desired amplitude function [k]
v	Weighting vector (if no weighting desired set to NULL) [k]
flags	Type flag (type 1-4), see filter.h

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

fd_cls()

使用最小二乘频率采样法设计非线性FIR滤波器。

Design non linear FIR filter using the least squares frequency sampling method.

The filter is designed by solving the equation

where W is a diagonal weighting matrix with the diagonal v, is generated from the frequency vector f and ad is the complex desired amplitude response. The complex frequency response is defined as

参数:

n	Filter length
w	Buffer for the filter taps [n]
k	Number of frequency sample points
f	Frequency vector [k] [0,0.5]
ad	Desired amplitude function [k]
gd	Desired group-delay [k]
v	Weighting vector (if no weighting desired set to NULL) [k]
ph	Phase offset [-0.5. 0.5]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

fd_clsgd()

使用最小二乘法频率采样法和群延迟抑制设计非线性FIR滤波器。

Design non linear FIR filter using the least squares frequency sampling method and group-delay constraint.

The filter is designed by solving the equation

where W is a diagonal weighting matrix with the diagonal v, phi and psi are generated from the frequency vector f and ad is the complex desired amplitude response. The group-delay constraints are given by the diagonal elements vg of the matrix U. The complex frequency response is defined as

Set for don‘t care values.

参数:

n	Filter length
w	Buffer for the filter taps [n]
k	Number of frequency sample points
f	Frequency vector [k] [0,0.5]
ad	Desired amplitude function [k]
v	Amplitude weighting vector (if no weighting desired set to NULL) [k]
gd	Group delay [k]
vg	Group-delay constraints [k]
ph	Phase offset [-0.5. 0.5]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

firremez.c

remes算法，用于设计有限脉冲响应滤波器。

fd_remez()

给出频率带边界、期望的响应和误差权重设计最优FIR滤波器。

Design the optimal (in the Chebyshev/minimax sense) FIR filter given a set of band edges, the desired reponse on those bands, and the weight given to the error in those bands.

参数:

n	Filter length
w	Filter taps [n]
k	Number of frequency bands in specification
f	Frequency band edges [0,0.5] [2 * k]
ad	Desired amplitude function [k]
v	Error weights [k]
flags	Type of filter, see filter.h

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

firminphase.c

设计带有实/复频特性和最小相位的FIR滤波器。

fd_minphase()

使用remes变换算法设计最优量值(在切比雪夫/极限判断里)的线性最小相位FIR滤波器。脉冲响应给出频率带边界、期望的响应和误差权重。滤波效果是较低的延迟但会有非线性相位。

Design optimum magnitude (in the Chebyshev/minimax sense) linear minimum-phase FIR filter using Remes exchange algorithm. The impulse response given a set of band edges, the desired response on those bands, and the weight given to the error in those bands. The resulting filter has very low delay but non linear phase.

警告:

Unpredictable result if v contains more than 2 different weight values.

参数:

n	Number of filter taps
w	Filter taps [n]
k	Number of frequency bands in specification
f	Frequency band edges [0,0.5] [2 * k]
ad	Desired amplitude function [k]
v	Error weights [k]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

fd_minphaseminimax()

利用谱分解设计带有改良极限准则的最小相位FIR滤波器。

Design minimum-phase FIR filter with modified minimax criterion using spectral factorization.

The filter is designed using the following specification:

where r is a symmetric intermediate filter of length , v is a weighting vector, is generated from the frequency vector f and ad is the amplitude function. The vector g represents frequencies for which the frequency response of the resulting filter will be zero. The index kl separates the passband and the stop-band. Frequency points with indexes below kl must be in the passband, and above kl in the stop-band. The solution to the above linear programming problem is found using simplex.

The minimum phase filter w is obtained from r using spectral factorisation.

警告:

This design method is extremely parameter sensitive exhaustive search over the number of filter taps may be required.

参数:

n	Filter length
w	Buffer for the filter taps [n]
k	Number of frequency sample points
kl	Border between passband and stop-band 0 <= kl <= k-1
f	Frequency vector [k] [0,0.5]
ad	Desired frequency response [k]
v	Weighting vector [k]
l	Number of null constraints
g	Null constraint vector [l] [0,0.5]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

fd_minphasels()

利用谱分解设计带有改良最小二乘准则的最小相位FIR滤波器。

Design minimum phase FIR filter with modified least squares criterion using spectral factorization.

The filter is designed using the following specification:

where r is a symmetric intermediate filter of length , is generated from the frequency vector f and au and al are the upper and lower limit for the amplitude function respectively. The vector g represents frequencies for which the frequency response of the resulting filter will be zero. The solution to the above linear programming problem is found using simplex.

The minimum phase filter w is obtained from r using spectral factorisation.

注解:

This method only designs low-pass filters.

警告:

This design method is extremely parameter sensitive exhaustive search over the number of filter taps may be required to find working solution.

参数:

n	Filter length
w	Buffer for the filter taps [n]
k	Number of frequency sample points
f	Frequency vector [k] [0,0.5]
au	Upper limit for desired frequency response [k]
al	Lower limit for desired frequency response [k]
fs	Stop band limit [0,0.5];
l	Number of null constraints
g	Null constraint vector [l] [0,0.5]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

fd_cepminphase()

利用倒频谱技术设计带有极限准则的最小相位FIR滤波器。该方法会使很长的滤波变快。

Design minimum phase FIR filter with minimax criterion using cepstrum technique. This method is good for designing very long filters fast.

参数:

n	Filter length
w	Buffer for the filter taps [n]
k	Number of frequency sample points limited to
au	Linearly spaced upper bound for frequency response [k/2+1]
al	Linearly spaced lower bound for frequency response [k/2+1]
t	Truncation tolerance (set to 0 for default value 0.5dB)

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

firminimax.c

使用线性编程设计带有实/复频极限特性的有限脉冲响应滤波器。

fd_minimax()

使用最小二乘频率采样法设计FIR滤波器。

The filter is designed using the following linear program specification:

where w is the resulting filter taps v is a weighting vector, is generated from the frequency vector f and ad is the amplitude function. The solution to the above linear programming problem is found using simplex.

参数:

n	Filter length must be odd for HP and BS filters
w	Buffer for the filter taps [n]
k	Number of frequency sample points
f	Frequency vector [k] [0,0.5]
ad	Desired amplitude function [k]
v	Weighting vector (if no weighting desired set to NULL) [k]
flags	Type flag (type 1-4), see filter.h

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

fd_cminimax()

使用复极限准则如Remes算法来设计非线性FIR滤波器。

Design non linear FIR filter with complex minimax specification see Remes filter design for filters with real specification.

The filter is designed using the following specification:

where w is the resulting filter taps v is a weighting vector, is generated from the frequency vector f and hd is the complex frequency response. The solution to the above linear programming problem is found using simplex.The complex frequency response is defined as

参数:

n	Filter length
w	Buffer for the filter taps [n]
k	Number of frequency sample points
f	Frequency vector [k] [0,0.5]
ad	Desired frequency response [k]
gd	Desired group delay [k]
v	Weighting vector (if no weighting desired set to NULL) [k]
ph	Phase offset [-0.5. 0.5]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

fd_cminimaxgd()

使用复极限准则和群延迟限制来设计非线性FIR滤波器。

Design non linear FIR filter with complex minimax specification and group-delay constraint.

The filter is designed using the following specification:

where w is the resulting filter taps v is a weighting vector, and psi are generated from the frequency vector f. The group-delay specification is given by gd[i]. The complex frequency response is defined as

Index values for which are excluded from the group-delay constraints. The solution to the above linear programming problem is found using the simplex algorithm.

参数:

n	Filter length
w	Buffer for the filter taps [n]
k	Number of frequency sample points
f	Frequency vector [k] [0,0.5]
ad	Desired amplitude function [k]
v	Weighting vector (if no weighting desired set to NULL) [k]
gd	Desired group-delay in samples [k]
vg	Group-delay weighting vector [k]
ph	Phase offset [-0.5. 0.5]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

firoptwin.c

使用最优窗来设计有限脉冲响应滤波器。
fd_qpoptwin()

使用最优窗的二次规划设计来设计FIR滤波器。

The filter is designed as

where

where and are a frequency matrices, a constant vector for the frequency band 0 and s is the maximum stop-band ripple. If s is too small the algorithm will not converge. The design problem is solved using quadratic programming.

注解:

This function only designs low-pass Type 1 and Type 2 linear FIR filters but could easily be extended to Type 3 and 4 filters as well.

参数:

n

Filter length

w

Buffer for the filter taps [n]

fs

Cutoff frequency [0,0.5]

s

Maximum stop-band ripple > 0

flags

Filter type according to filter.h

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.
fd_minimaxoptwin()

使用最优窗的极限设计来设计FIR滤波器。

where

where is a frequency matrices, a constant vector for the frequency band 0 and the objective function for the optimisation. The design problem is solved using simplex.

注解:

This function only designs low-pass Type 1 and Type 2 linear FIR filters but could easily be extended to Type 3 and 4 filters as well.

参数:

n

Filter length

w

Buffer for the filter taps [n]

fs

Cutoff frequency [0,0.5]

flags

Filter type according to filter.h

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.
fa_erroroptwin()

使用最优窗方法来计算FIR低通滤波器的停止带能量。

Calculate stop-band energy for a low-pass FIR filter designed using optimum window method.

参数:

n

Number of filter taps

w

Filter taps [n]

fs

Stop-band frequency

e

Error [1] (return value)

flags

Filter type according to filter.h

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.
fireq.c

fd_cminimaxeq()

Design FIR phase equaliser with complex specification using minimax optimisation criteria.

The filter is designed using the following specification:

for and

where w is the resulting filter taps, is generated from the frequency vector f, g is the complex channel to equalise, sigma is the stop-band limit and hd is the complex frequency response. The solution to the above linear programming problem is found using simplex.

参数:

n	Filter length must be odd for HP and BS filters
w	Buffer for the filter taps [n]
kp	Number of frequency sample points in passband
fp	Passband frequency vector [kp] [0,0.5]
hd	Desired frequency response [kp+ks]
g	Channel to equalise [kp]
ks	Number of frequency sample points in stop-band
fs	Stop-band frequency vector [ks] [0,0.5]
s	The stop-band limit sigma [ks]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

fd_cminimaxgdeq()

Design FIR phase equaliser with complex specification using minimax optimisation criteria and group delay constraints.

The filter is designed using the following specification:

for and

where w is the resulting filter taps, psi is generated from the frequency vector f, g is the complex channel to equalise, sigma is the stop-band limit and hd is the complex frequency response. The group-delay specification is given by gd[i] and is applied in the pass-band. The group delay constraining vector psi is calculated according to

The solution to the above linear programming problem is found using simplex.

参数:

n	filter length must be odd for HP and BS filters
w	buffer for the filter taps [n]
kp	number of frequency sample points in pass-band
fp	pass-band frequency vector [kp] [0,0.5]
hd	desired frequency response [kp]
g	channel to equalise [kp]
ggd	group delay for the channel to equalise [kp]
gth	phase response for the channel to equalise [kp]
gd	desired group delay in samples [kp]
vg	group delay weighting vector [kp]
ks	number of frequency sample points in stop-band
fs	stop-band frequency vector [ks] [0,0.5]
s	the stop-band limit sigma [ks]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

firfb.c

fd_halfband()

设计最小相位完美复现半袋滤波器。

Design minimum phase prefect reconstruction halfband filterbank

参数:

n	Filter length for each of the analysis and synthesis
h0	Analysis low-pass filter
h1	Analysis high-pass filter
f0	Synthesis low-pass filter
f1	Synthesis high-pass filter
fc	Cutoff-frequency for halfband filter

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

fd_dftfbadv()

Design of oversampled DFT filterbank, advanced interface. This function designs the prototype filter and creates workspaces. The number of subbands k must be 4x bigger than p.

参数:

k	Number of subbands
p	Oversampling factor 1 for critically down-sampled filterbank
l	Length of prototype filter l=n*k where n > 1
h0	Prototype filter [l] (if NULL this function will design one)
fb	Pointer to filterbank struct
stride	Time data stride
offset	Time data offset
flags	Analysis or synthesis filterbank, see filter.h

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

fd_dftfb

Design of oversampled DFT filterbank. This function designs the prototype filter, creates workspaces.

参数:

k	Number of subbands
p	Oversampling factor 1 for critically downsampled filterbank
l	Length of prototype filter where n > 1
h0	Prototype filter [l] (if NULL this function will design one)
fb	Pointer to filterbank struct
flags	Analysis or synthesis filterbank, see defines

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

ff_dftfb()

Free resources for oversampled DFT filterbank.

参数:

fb

Pointer to filterbank struct

fd_pdftfb()

Design of parallel oversampled DFT filterbank. This function designs the prototype filter and creates workspaces. The number of subbands k must be 4x bigger than p.

参数:

k	Number of subbands
p	Oversampling factor 1 for critically downsampled filterbank
l	Length of prototype filter where n > 1
i	Number of channels
h0	Prototype filter [l] (if NULL this function will design one)
fb	Pointer to filterbank struct
flags	Analysis or synthesis filterbank, see defines

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

ff_pdftfb()

Free resources for oversampled parallel DFT filterbank.

参数:

fb

Pointer to filterbank struct
firtrans.c

用于FIR变换。

Functions for performing transformations on finite impulse response filters.

ft_parallel()

把标准FIR滤波器变换为多相FIR滤波器。

Transform a prototype FIR filter into polyphase FIR filter.

参数:

n	Length of prototype filter
k	Number of polyphase components
w	Prototype filter taps
pw	Parallel FIR filter
g	Filter gain
flags	FIR_PARALLEL_FWD forward indexing FIR_PARALLEL_REW reverse indexing FIR_PARALLEL_ODD multiply every 2nd filter tap by -1, this results in a high-pass filter

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

ft_quantize()

Quantise FIR filter taps from _ftype_t to low precision format. The quantisation algorithm calculates the maximum gain for the input signal of the quantised filter it is returned in g, and may be < 1.

参数:

n

Number of filter taps

in

Input vector [n]

out

Quantised output vector [n]

g

Gain factor (optional set to null if not desired) [1]

type

Output type, see filter.h

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.
iiranalog.c

Classical analog filter design methods for Butterworth, Chebyshev I, Chebyshev II and Elliptic filters.

The length of the numerator and denominator polynomials (vectors) returned from the design functions is n3/2 if n is even and n3/2-1 if n is odd, where n is the order of the filter.

The below implementations are based on:

T. W. Parks and C. S. Burrus, Digital Filter Design, John Wiley & Sons, 1987, chapter 7.
id_bessel()

Calculate Bessel-Thomson polynomials for realisation in cascade or parallel form.

Numerator is always 1.0 calculation of denominator polynomial is performed using the following iterative method:

The first and second order sections are found using root solving.

参数:

n

Filter order

num

Polynomials for numerator

den

Polynomials for denominator

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.
id_butterworth()

Calculate Butterworth polynomials for realisation in cascade or parallel form.

For n even

for n odd

where F(s) is the Laplace transform of the filter.

参数:

n

Filter order

num

Polynomials for numerator

den

Polynomials for denominator

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.
id_chebyshev()

Calculate Chebyshev polynomials for realisation in cascade or parallel form.

For n even

for n odd

where F(s) is the Laplace transform of the filter and

where

e is part of the function call and can be calculated according to

where a is the desired passband ripple in dB.

参数:

n

Filter order

num

Polynomials for numerator

den

Polynomials for denominator

g

Over all filter gain (input and return value) [1]

e

Design parameter, see above

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.
id_ichebyshev()

Calculate inverse Chebyshev polynomials for realisation in cascade or parallel form.

For n even

for n odd

where F(s) is the Laplace transform of the filter and

where

e is part of the function call and can be calculated according to

where b is the desired passband ripple in dB.

参数:

n

Filter order

num

Polynomials for numerator

den

Polynomials for denominator

g

Over all filter gain (input and return value) [1]

e

Design parameter, see above

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

id_elliptic()

设计椭圆滤波器。参数e和k1的计算公式如下：

Elliptic function filter design. The design parameters e and k1 and be calculated according to

a是通带波纹，b是阻带衰减。

where a is passband ripple in dB and b is stop-band attenuation in dB.

参数:

n	滤波器阶数
num	生成滤波器多项式分子
den	生成的滤波器分母
g	增益
e	Design parameter see above
k1	Second modulus of k design parameter, see above

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

iirtrans.c

模拟和数字IIR滤波器变换。

Analog and digital infinite impulse response filter transformations.

The length of the numerator and denominator polynomials (vectors) used to represent the analog filter polynomials is n3/2 if n is even and n3/2-1 if n is odd, where n is the order of the filter, unless otherwise noted.

it_lp2lp()

低通到低通模拟滤波器变换（改变截断频率）。

Low-pass to low-pass transformation of analog filter (changes cutoff frequency).

参数:

n	Filter order
num	Polynomials for numerator
den	Polynomials for denominator
fc	Cutoff frequency [0,]
g	Over all filter gain (input and return value) [1]
numt	Polynomials for transformed numerator
dent	Polynomials for transformed denominator

注解:

dent and den can be the same pointer same goes for numt and num.

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

it_lp2hp()

模拟滤波器低通到高通变换。

Low-pass to high-pass transformation of analog filter.

参数:

n	Filter order
num	Polynomials for numerator
den	Polynomials for denominator
fc	Cutoff frequency [0,]
g	Over all filter gain (input and return value) [1]
numt	Polynomials for transformed numerator
dent	Polynomials for transformed denominator

注解:

dent and den can be the same pointer same goes for numt and num.

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

it_lp2bp()

模拟滤波器低通到带通变换。

Low-pass to band-pass transformation of analog filter.

参数:

n	Filter order
num	Polynomials for numerator
den	Polynomials for denominator
fl	Lower cutoff frequency [0,fh]
fh	Upper cutoff frequency [0,]
g	Over all filter gain (input and return value) [1]
numt	Polynomials for transformed numerator [n*3]
dent	Polynomials for transformed denominator [n*3]

注解:

This function doubles the order of the filter, numt and dent are therefore twice the length of num and den.

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

it_lp2bs()

模拟滤波器低通到带阻变换。

Low-pass to band-stop transformation of analog filter.

参数:

n	Filter order
num	Polynomials for numerator
den	Polynomials for denominator
fl	Lower cutoff frequency [0,fh]
fh	Upper cutoff frequency [0,]
g	Over all filter gain (input and return value) [1]
numt	Polynomials for transformed numerator [n*3]
dent	Polynomials for transformed denominator [n*3]

注解:

This function doubles the order of the filter, numt and dent are therefore twice the length of num and den.

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

it_bilinear()

使用双线性变换设计IIR滤波器。将模拟滤波器转换成数字IIR滤波器。

IIR filter design using bilinear transform. Transform a series of analog 2nd order filters to a series of IIR biquad links using bilinear transformation.

Q is filter quality factor or resonance, in the range of 1 to 1000. The overall filter Q is a product of all 2nd order stages. For example, the 6th order filter (3 stages, or biquads) with individual Q of 2 will have filter . The critical frequency for the analog filter is assumed to be 1Hz.

The parameters for the sections are stored according to:

where is the filter section index. For odd length filters c will be prepended by a first order section.

参数:

n	滤波器阶数
num	S域多项式分子系数
den	S域多项式分母系数
Q	Q value for the filter [1,1000]
fc	临界频率 [0,0.5]
g	增益
c	z域的滤波器系数

注解:

Upon return, the g argument will be set to a value, by which to multiply our actual signal in order for the gain to be the desired passband gain.

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

it_bilinear_adv()

双线性变换的加强版，提供采样频率的参数输入。

IIR filter design using bilinear transform - advanced interface. Transform a series of analog 2nd order filters to a series of IIR biquad links using bilinear transformation.

The parameters for the sections are stored according to:

where is the filter section index. For odd length filters c will be prepended by a first order section.

参数:

n	Filter order
num	s-domain numerator coefficients
den	s-domain denominator coefficients
Q	Q value for the filter [1,1000]
f	s-domain filter critical frequency [0,inf]
fs	Sample frequency [0,inf]
fz	z-domain filter critical frequency [0,fs/2]
g	Filter gain factor, (input and return value, initialise to desired passband gain) [1]
c	Array of z-domain coefficients to be filled in [4*n/2]. The coefficients are ordered according to a1, a2 (denominator) b1, b2 (numerator) b0 and a0 are always 1.

注解:

Upon return, the g argument will be set to a value, by which to multiply our actual signal in order for the gain to be the desired passband gain.

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

it_sos2poly()

数字二阶系统到多项式转换。

Digital second order system to polynomial.

参数:

n

Filter order

c

it_bilinear()输出

g

增益

a

多项式分子系数 [n+1]

b

多项式分母系数 [n+1]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

it_sos2allpass()

Transform second order system to sum of two all-pass filers. Each all-pass filer is represented as a linked list of first and second order all-pass sections.

参数:

n	Filter order
c	Filter coefficients [n*4/2] stored according to The output from it_bilinear()
g	Filter gain
a1	Linked list of all-pass sections number 1
a2	Linked list of all-pass sections number 2

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.

it_quant_allpass()

Quantise the parameters of an all-pass filter to n bits.

参数:

a

Linked list of all-pass sections

n

Number of bits [1,64]

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.
it_allpass2poly()

List of all-pass sections to polynomial. The vectors a and b are allocated by this function and must be freed by caller.

参数:

ap

Linked list of all-pass sections

n

Filter order

a

Denominator polynomial [n+1]

b

Numerator polynomial [n+1]

flags

According to filter.h

返回:

On success 0 is returned, -1 is returned on error and errno is set appropriately.
matrix.h/matrix.c

矩阵代数
mat_calloc()

Allocate space for multidimensional matrix. The allocated matrix will have its elements allocated in one single continous block of data.

wz size in bytes of each element in the allocated matrix n number of dimensions ... size of each inividual dimension [a1, a2, ..., an] in number of elements, size of returned matrix will be [a1, a2, ..., an].

return pointer to matrix data if success, NULL if fail.

参数:

wz

n

…

返回:

return pointer to matrix data if success, NULL if fail.
mat_free()

Free matrix data. This function is recursive.

n number of dimensions m matrix

参数:

n

m

返回:

return void.
mat_conv()

Convolve the two vectors z and y and return the result in z.

参数:

n

m

x

y

z

返回:
mat_cconv()

mat_conv()的复数形式。

参数:

n

m

x

y

z

返回:
mat_flip()

Flip vector.

参数:

n

x

y

返回:
mat_cflip()

mat_flip()的复数形式。

参数:

n

x

y

返回:
mat_mult()

Multiply matricies

参数:

n

m

l

x

y

z

flag

返回:
mat_cmult()

mat_mult()的复数形式。

参数:

n

m

l

x

y

z

flag

返回:
mat_vmult()

Multiply vectors

参数:

n

m

x

y

返回:
mat_cvmult()

mat_vmult()的复数形式。

参数:

n

m

x

y

flag

返回:
mat_add()

Add vectors (and matricies)

参数:

n

m

x

y

z

返回:
mat_cadd()

mat_add()的复数形式。

参数:

n

m

x

y

z

flag

返回:
mat_sub()

Subtract vectors (and matricies) z = x-y

参数:

n

m

x

y

z

返回:
mat_csub()

mat_sub()的复数形式。

参数:

n

m

x

y

z

flag

返回:
mat_omult()

Elementwise vector (and matrix) multiplication

参数:

n

m

x

y

z

返回:
mat_comult()

mat_omult()的复数形式。

参数:

n

m

x

y

z

flag

返回:
mat_odiv()

Elementwise vector (and matrix) division

参数:

n

m

x

y

z

返回:
mat_codiv()

mat_odiv()的复数形式。

参数:

n

m

x

y

z

flag

返回:
mat_real()

Real and complex part of complex vector (or matrix)

参数:

n

x

y

返回:
mat_imag()

Real and complex part of complex vector (or matrix)

参数:

n

x

y

返回:
mat_conj()

Complex conjugate of vector (or matrix)

参数:

n

x

y

返回:
mat_inv()

Matrix inversion

参数:

n

x

y

返回:
mat_cinv()

mat_inv()的复数形式。

参数:

n

x

y

返回:
mat_eye()

Identity matrix

参数:

n

x

返回:
mat_ceye()

mat_eye()的复数形式。

参数:

n

x

返回:
mat_solve()

Solve A*x=b

参数:

n

m

a

b

x

ws

返回:
mat_csolve()

mat_solve()的复数形式。

参数:

n

m

a

b

x

ws

返回:
mat_linspace()

生成一个范围在min和max之间的等间距数组，一般作为组成频率响应图的点的x分量。

Generate n linearly spaced values in the range min < f < max

参数:

n

频率数组点数

min

频率下限

max

频率上限

x

生成的频率数组

返回:

return -1 if fail 0 if success.
mat_clinspace()

Generate n complex linearly spaced values in the range min < f < max

参数:

n

min

max

x

返回:
mat_clinspace_pol()

Generate n linearly spaced complex values specified in polar form

参数:

n

amin

amax

wmin

wmax

x

返回:
Libfilth的windows移植

libfilth安装说明有如下描述：

Dependencies

To compile and run the library you need:

* GSL GNU Scientific Library 1.4

* lp-solve Solve (mixed integer) linear programming problems 4.0

* FFTW Fastest Fourier Transform in the West 3.0

* BLAS Basic Linear Algebra Subroutines 1.1

上面说的是依赖关系，要在linux下编译libfilth需要几个软件库的支持，它们是：

* GSL GNU 科学计算库 1.4

* lp-solve 解决 (混合整数) 线性规划问题 4.0

* FFTW 快速傅立叶变换库 3.0

* BLAS 基本线性算法库（执行向量和矩阵运算的子程序集合） 1.1

这几个依赖的软件库，基本上都是基于linux核GCC开发的。

由于笔者对数字滤波器的相关知识约等于零，所以不去理解libfilth，先想办法编译过去再说，把libfilth拆了，在VS2005里新建工程，一个一个文件往里添。

添加的过程就不说了，列出遇到的问题和解决办法吧。

复数的表达使用问题。gcc支持C语言的复数运算，而VC不提供C语言的复数运算支持，虽然支持complex的变量定义，但不支持这种定义的复数运算，实际上看complex的库定义，它只是一个包含两个实数元素的结构体来代表实部和虚部。如果将程序中的复数全部用这种结构体替代，那么程序中的复数运算将是一个灾难性的工作量。经过几天的研究，发现VC的C++的库中有对复数的支持，于是将libfilth中所有的复数表示替换成c++的声明定义形式，但结果却是无法编译，经研究将libfilth的所有.c文件后缀改为.cpp(.c的后缀VC会进行C编译，因此不识别c++的复数定义)。
lp-solve库的调用。基础功能文件internal.cpp 中有对lp-solve库的调用，而internal.cpp 相当于libfilth内部的函数库，基本所有的文件都要调用其中的函数。在sourceforge找到了lp-solve 5.5版本的源代码，提供对VS2005的支持，集成进本项目，编译后，以静态库的形式调用，成功。
GSL库的调用。同上，也是在internal.cpp中被调用，起初在网上找到了chgsl，不知道是GSL的什么版本，只提供静态库和头文件，不提供源文件，可以集成进项目，但是在编译过程中会出现LNK2005的链接错误，于是寻找了GNU的标准GSL库，它只提供linux下的版本，最终找到非官方的win32下的GSL库源代码，提供对VS2005的支持，并且包含有BASL的支持，而且libfilth也是通过GSL间接使用BLAS的，至此所有依赖库的问题都解决了。
asinh函数的移植，gcc提供反双曲正弦函数，而VC不提供，从网上下载asinh实现源代码集成进项目。
GCC和VC的库函数名称差异。两者库函数有许多地方函数名称存在差异，一般在GCC中支持复数运算的函数都和相通功能的实数运算函数区分开，但是在windows下通用；复数取共轭的表示法不同；使用条件编译等方法使源程序在GCC下和VC下都能够编译，详见filth.h的实现。
I的使用。GCC用I表示单位虚部，但VC不支持，只好类型定义一个单位虚数来代替。
memalign的支持。GCC支持该函数的使用，但VC不提供同样功能的函数，没办法，用malloc代替，还不知道有什么后果。
inline函数的使用。VC中C++程序inline函数的声明和定义必须放在一起(在头文件中写函数体)，否则会出现LNK2001的链接错误。
在调用FFTW函数的时候，libfilth使用了float级精度的调用接口，而在windows下编译的FFTW是double的精度接口，而在编译float精度的接口FFTW库时出错，况且libfilth本身的运算都是double精度的，因此将libfilth的所有调用FFTW的接口改为double精度，并且不会对程序造成任何影响。
编译动态库链接出错。编译静态库通过，而编译动态库出现LNK2001和LNK2005的错误，其中有上面几点的部分原因，另外注意所调用的库编译选项中的运行时库选择要一致，这里选择的是：多线程(/MT)。
C++不支持隐式强制类型转换，必选显示表示。
typeof的支持，GCC支持typeof()扩展，VC有typeid().name()与之对应，但在编译过程中出错，由于有RTTI方面等诸多问题，决定放弃使用。程序使用typeof()也是为了使程序能够使函数兼容实数和复数两种输入参数，那我就把这两种参数拆开变成调用两种不同的内容，虽然麻烦点，但是能用，从这方面也能够看出跨平台的程序绝对不能有编译器的扩展的库函数或方法。
在.cpp文件中不能在函数下方去声明函数参数类型，必须在函数括号里的参数前直接声明，这可能是VC在进行C++编译时不支持这种写法。

解决了如上的诸多问题，libfilth终于能够生成dll文件，至于各功能好不好用，由于libfilth版本号0.4(看着就玄)而且在移植过程中对代码有了较大的调整，所以还有待下一步的检验和测试。

libfilth的测试和使用分析

libfilth已经在VS2005下编译成功了，接下来就是使用测试了，笔者准备了labview来进行这项工作，labview在数字滤波器的功能上是比较全的，而且有着丰富强大的绘图显示功能，把libfith做成dll放在labview中去使用，不但能测试libfilth的所有功能函数，而且能将其功能与labview本身滤波器功能进行详细的对比，在此过程中既能检验libfith的功能，摸清其特点，又能加深对齐功能函数的理解。

libfilth的所有功能接口都被集中在filter.h中，因此分析filter.h的功能函数就可以了。

libfilth给出了函数库的使用说明和示例，下面我们就一边翻译其文档一边熟悉和测试libfilth。
设计IIR滤波器

有限脉冲响应（IIR）滤波器可以由使用双线性变换的模拟滤波器设计。本章描述五类模拟滤波器的设计方法和怎样使用双线性变换把它们装换成IIR滤波器。本章也会涉及如何改变模拟滤波器的截止频率并且怎样把它转换成带通、带阻、高通和低通滤波器。为了方便使用级联biquads处理精确的应用（二阶无限脉冲响应滤波器经常叫作‘biquads‘），这里对滤波器的描述将被分成两节。设滤波器阶数为N。

模拟滤波器的设计

五类经典模拟滤波器是（贝塞尔-汤姆森）Bessel-Thomson, （巴特沃斯）Butterworth, （切比雪夫）Chebyshev, （反切比雪夫）inverse Chebyshev, and （椭圆）elliptic function 滤波器，其描述如下：
1. Bessel-Thomson滤波器在通带上有着近乎恒定的群时延，使它们很适合于音频应用。恒定的群时延牺牲的是一个较宽的过渡带。滤波器的频率响应只由N控制。选择合适的长度没有简单的方法，但是可以通过下面表格的公式选取相近值。
2. Butterworth滤波器在ω = 0和ω = ∞最平滑，并且有着非常平滑的频率和相位响应。滤波器的频率响应仅由N决定，近似的取值为：
(4-1)

滤波器的频率响应介于unity之间，0 < ω < ωp，ωp < 1 且 0 < G < 1。
1. Chebyshev通过参数ε控制通带中的最大纹波。频率响应在ω = ∞最平滑。参数ε可以通过下面的公式获得：
(4-2)

δ1是通带纹波的期望值，a是通带纹波信噪比，滤波器的序列数可有下面的公式获得：

(4-3)

滤波器的频率响应介于unity之间，0 < ω < ωs，ωs < 1 且 0 < G < 1。
1. inverse Chebyshev通过参数k1实现高阻带衰减，频率响应在ω = ∞最平滑。参数k1可以通过下面的公式获得：
(4-4)

δ2是阻带纹波的期望值，b是阻带衰减的信噪比，滤波器的序列数可有下面的公式获得：

(4-5)

滤波器的频率响应介于unity之间，0 < ω < ωp，ωp < 1 且 0 < G < 1。
1. Elliptic Function 滤波器特点由一下四个参数确定：
  1. 通带纹波δ1
  2. 过渡带宽度（1-ωs）
  3. 最大阻带纹波δ2
  4. 滤波器阶数N
给定上面的三个参数，第四个将会很小。设计方法由参数ε 和 k1确定。参数的计算公式如下：

(4-6)

a是通带纹波信噪比，b是阻带衰减的信噪比，对于给定采样次数N的滤波器，使用上述方程过渡区将会很小。

上述设计思想由iiranalog.cpp文件中的函数：id_bessel() ，id_butterworth() ，id_chebyshev() ，id_ichebyshev()，和id_elliptic() 实现。这些功能函数所实现的滤波器被称为原型滤波器。拥有1Hz的截止频率和单位增益的通带。

只有Besel-Thomson和Butterworth滤波在通带上有单位增益，其它三种滤波器的通带增益由不同的参数决定。增益漂移由设计函数进行补偿。

模拟滤波器变换

使用频率变换可以改变原型滤波器的截止频率和响应。下面是四个基本的变换：
1. 低通到低通变换改变原型滤波器的截止频率为F。
(4-7)
1. 低通到高通变换把低通原型滤波器转换为截止频率为F的高通滤波器。
(4-8)
1. 低通到带通变换把低通原型滤波器转换为截止频率为Fl和Fh的带通滤波器。
(4-9)
1. 低通到带阻变换把低通原型滤波器转换为截止频率为Fl和Fh的带阻滤波器。
    (4-10)

这些变换的函数实现由iirtrans.cpp中的it_lp2lp()， it_lp2hp()， it_lp2bp()和it_lp2bs()提供。

双线性变换

双线性变换用于将模拟滤波器变换为相同性能的数字IIR滤波器。把模拟滤波器的拉普拉斯变换的左侧半平面变换为数字滤波器单位圆内的Z变换。其变换结果是一个同模拟滤波器结果相同的稳定的IIR滤波器。使用下面的公式获得该数字IIR滤波器：

    (4-11)

F ∈ [0,∞]是模拟滤波器的临界频率，Fs ∈ [0,∞]是采样频率，Fz ∈[0, Fs/2]是数字IIR滤波器的临界频率。把频率F变换到Fz的公式：

    (4-12)

当滤波器设置了一个以上的临界频率时注意频率扭曲和细致的补偿。如果模拟滤波器的采样频率和临街频率都等于1，公式可以简化为：

    (4-13)

fz ∈ [0, 0.5]。

变换的实现被分成两部分。高级和简易的变换实现对应iirtrans.cpp的it_bilinear_adv()和it_bilinear()两个函数。
处理IIR滤波器设计的频率响应

使用上述模拟滤波器设计的IIR滤波器遵循双线性变换。五种滤波器符合下面的规范：

五个滤波器的必要参数列于图4_1。

滤波器的频率响应见图4_1。贝塞尔-汤姆森滤波器不符合规格，通带的频率响应见图4_1，图中显示除了贝塞尔-汤姆森滤波器的其余4个都符合规格。

图4_1中显示了群延迟。可以看出切比雪夫滤波器有着最差的延迟，巴特沃斯和椭圆滤波器次之。最小的群延迟是反切比雪夫和贝塞尔-汤姆森滤波器。

在l8a.cpp中有实现。图4.2是l8a.cpp生成的五种滤波器频率响应数据由labview绘制的图像，基本与图4.1一致。图4.3是l8a.cpp生成五种滤波器的群延迟数据由labview绘制的图像。

图 41（elliptic低通滤波器的频率响应）

图 42

图 43
处理IIR滤波器设计的阶数

由上述方法设计的的N=10的IIR滤波器，受下面的公式规范：

这五种过滤器所需的设计参数列于图4.4。

图 44(滤波器设计参数)

滤波器的频率响应见图4.5。图4.5显示的是滤波器的通带频率响应，图中除了巴特沃斯和贝塞尔-汤姆森滤波器，剩下的都满足上面公式规范的要求。

贝塞尔-汤姆森滤波器的群延迟最低，椭圆滤波器的群延迟最差。

程序实现见l8b.cpp。图4.6和图4.7分别是程序生成的滤波器的频率响应和群延时数据由labview绘制的图。

图 45(低通IIR滤波器的通带频率响应)

图 46

图 47
IIR滤波器实现

实现高效的IIR滤波器比较困难。图4.9是一个阻带衰减为-40dB，截止频率ωc = 0.4π的8阶滤波器频率响应。

程序实现见l8c.cpp。(原英文文档图片、程序和文字对不上号，严重有问题)

图4-10是正弦波和白噪声的叠加之后通过libfilth设计的椭圆滤波器和labview设计的椭圆滤波器的效果。信号源相同，滤波器设计参数近似(两者的设计给定参数形式不同，得到的的滤波器参数会有出入)，输出数据由labview统一绘制成波形，由图中可以看出labview的效果要好，libfilth效果不如labview，但是也证明了libfilth的IIR滤波器是有效的，而且IIR滤波器的设计还有待进一步理解，所以libfilth还有挖掘的潜力。

l8c.cpp程序中设计了一个IIR滤波器实现的宏：

#define IIR(in,w,q,out) { \

double h0 = (q)[0]; \

double h1 = (q)[1]; \

double hn = (in) - h0 * (w)[0] - h1 * (w)[1]; \

out = hn + h0 * (w)[2] + h1 * (w)[3]; \

(q)[1] = h0; \

(q)[0] = hn; \

}

输入数据和设计好的滤波器通过该宏的运算得到数据输出，具体应用执行请参考该文件。

图 48

图 49

图 410(绿：输入，青：libfilth，红：labview)
模拟滤波器变换

一个5阶的通带波纹1dB、阻带衰减40dB的椭圆滤波器。

l8d.cpp实现生成的数据由labview绘图。

图 411

图 412

图 413

图 414

图 415

图 416

图 417

图 418
IIR设计实现总结

本章介绍了libfilth提供的IIR滤波器的几乎所有功能，示例程序也描述了怎样使用libfilth来构建一个IIR滤波器并进行滤波器分析和怎样使用滤波器进行滤波，下面对这一过程的实现步骤做简单总结：

1.准备好设计滤波器的条件参数，包括：滤波器阶数、通带波纹、阻带衰减、高低截止频率等。

2.将参数带入五种滤波器设计函数，包括：id_bessel()，id_butterworth()，id_chebyshev()，id_ichebyshev()，id_elliptic()，经计算得到滤波器多项式的系数。

3.进行频率变换，上一步得到的滤波器是低通滤波器，需要通过频率变换将其转换成高通、带通和带阻滤波器，包括：it_lp2lp()，it_lp2hp()，it_lp2bp()，it_lp2bs()。

4.将得到的滤波器系数代入it_bilinear()( 它会调用it_bilinear_adv())进行双线性变换，将模拟滤波器转换成数字滤波器，得到Z域的数字滤波器系数。

5.处理信号，将输入信号数据和数字滤波器系数代入IIR滤波器实现宏(见l8c.cpp)，得到滤波后的数据。

6.计算滤波器的频率响应和群延迟，第4步得到的数字滤波器系数转换成多项式系数的形式( it_sos2poly())，然后代入IIR分析函数( ia_response()，期间会调用mat_linspace()来计算输出的响应数据的间隔，即波形图横轴)计算频率响应和群延迟。如果想直接分析双线性变换前的模拟滤波器的频率响应、相位响应等，可以使用模拟滤波器分析函数( ia_sresp())。
使用Remes算法设计线性FIR滤波器

使用Remes算法能够设计线性相位FIR滤波器。本章使用Remes算法来解决各种滤波器的设计问题。下面是算法的执行步骤。
使用Remes变换算法求解多项式

三阶多项式利用二阶多项式求解，使用初始设置T={-0.8,-0.2,0.2,0.8}。首次迭代给出T={-1,-0.3,0.3,1}，二次迭代给出最优解T={-1,-0.5,0.5,1}，三次迭代的误差E(ω)和最优解曲线见图4-19。

图 419

产生的振幅结果连同期望振幅函数见图4-20，实现见Octave程序l4a.m。

图 420
利用Remes算法计算FIR滤波器

手动执行Remes算法设计一个1型5阶的FIR滤波器，期望振幅函数。一个5阶FIR滤波器的振幅函数为。该滤波器设计初始值为T={0,0.35π,0.5π,0.9π}。经首次迭代给出T={0,0.15π,0.5π,0.8π}，二次迭代的最优解是T={0,π/4,π/2,3π/4}，三次迭代的误差函数E(ω)见图4-21。产生的振幅结果连同期望振幅函数见图4-22，实现见Octave程序l4b.m。

图 421

图 422

手动执行Remes算法设计一个1型7阶的FIR滤波器，期望振幅函数:

(4-14)

一个5阶FIR滤波器的振幅函数为。该滤波器设计初始值为T={0.1,0.75,1.5,2.3,3}。经首次迭代给出T={0.5,1.2,1.7,2,2.7}，二次迭代的最优解是T={π/6,π/3,π/2,2π/3, 5π/6}，三次迭代的误差函数E(ω)见图4-23。产生的振幅结果连同期望振幅函数见图4-24，实现见Octave程序l4c.m。

图 423

图 424
利用Remes算法计算FIR低通滤波器

利用Remes算法设计一个1型21阶的FIR低通滤波器，期望振幅函数:

(4-15)

产生振幅见图4-25。实现函数见l4d.cpp，设计中调用了firremez.cpp中的fd_remez()。

图 425

该程序由labview调用生成图4-26。

图 426
利用Remes算法计算FIR高通滤波器

利用Remes算法设计一个1型21阶的FIR高通滤波器，期望振幅函数:

(4-16)

产生振幅见图4-27。实现函数见l4e.cpp，设计中调用了firremez.cpp中的fd_remez()。

图 427

该程序由labview调用生成图4-28。

图 428
Remes算法设计FIR滤波器测试总结

图4-29是正弦波和白噪声的叠加之后通过libfilth使用Remes算法设计的低通滤波器和labview设计的椭圆滤波器的效果，随便选的参数，libfilth的效果不太好。

图 429(绿：输入，青：libfilth，红：labview)

图4-30效果好一点，通过调整输入数据，采样点不变的情况下增加了采样信号正弦量的周期。

图 430(绿：输入，青：libfilth，红：labview)

图4-31为使用Remes算法设计的FIR低通滤波器(通带截止频率：10Hz,阻带截止频率20Hz)的信号处理效果比较。输入信号成分：采样频率100Hz，5Hz正弦(0起始相位，单位幅值)+30Hz正弦(100°起始相位，单位幅值)+高斯白噪声(15%单位幅值)。

图 431

图4-32是将5Hz的输入信号成分改成1Hz的输出效果。

图 432

经过一星期的调试，最后将remes算法的FIR滤波器与labview相应的滤波器参数统一，得到的滤波效果完全一致。
使用最优窗设计FIR滤波器
数字全通自适应递归滤波器

有些滤波器可以理解为两个全通滤波器之和，见图4-31。这些滤波器特别是特别是在通带上，参数量化具有鲁棒性。另一个好处是可以实现很小的乘数。这两种特性使得滤波器适合于硬件实现。

图 433

滤波器的一个充要条件是：

    (4-14)

作为两个全通滤波器和的实现，存在一个无功率滤波器

    (4-15)

使得

    (4-16)

给出特征函数

    (4-17)

一个N阶的IIR滤波器使用下面的步骤变换：
1. 求多项式Q(z)。
  
  由R(z)
    (4-18)

R(z)逆Z变换：

    (4-19)

反过来给我们递归函数

    (4-20)
1. 求两个全通滤波器A1(z) 和 A2(z)。
(4-21)
1. 量化滤波器系数

时间： 2024-11-05 17:23:16

n	过滤阀数，where mod
w	过滤阀
x	输入信号必须是一个顺序索引的环形缓冲区

n	Window length
w	Buffer for window parameters [n]
b	Beta parameter for Kaiser window, Beta >= 1

n	Filter length
w	Buffer for the filter taps [n]
fs	Cutoff frequency [0,0.5]
s	Maximum stop-band ripple > 0
flags	Filter type according to filter.h

n	Number of filter taps
w	Filter taps [n]
fs	Stop-band frequency
e	Error [1] (return value)
flags	Filter type according to filter.h

n	Number of filter taps
in	Input vector [n]
out	Quantised output vector [n]
g	Gain factor (optional set to null if not desired) [1]
type	Output type, see filter.h

n	Filter order
c	it_bilinear()输出
g	增益
a	多项式分子系数 [n+1]
b	多项式分母系数 [n+1]

ap	Linked list of all-pass sections
n	Filter order
a	Denominator polynomial [n+1]
b	Numerator polynomial [n+1]
flags	According to filter.h

n	频率数组点数
min	频率下限
max	频率上限
x	生成的频率数组

n	窗长
w	窗参数缓冲区 [n]

Libfilth(一个滤波器C库)使用

版本历史修改记录

概述

库文件分析

filth.h/filth.c

quantize()

filter_strerror()

filter.h/filter.c

fir()

pfir()

updatepq()

window.h/window.c

wd_boxcar()

wd_triang()

wd_hanning()

wd_hamming()

wd_blackman()

wd_flattop()

wd_kaiser()

filtanalysis.c

fa_ampfunc()

fa_response()

fa_cerror()

ia_response()

ia_sresp()

ia_response_allp()

firwindow.c

fd_window()

firls.c

fd_ls()

fd_cls()

fd_clsgd()

firremez.c

fd_remez()

firminphase.c

fd_minphase()

fd_minphaseminimax()

fd_minphasels()

fd_cepminphase()

firminimax.c

fd_minimax()

fd_cminimax()

fd_cminimaxgd()

firoptwin.c

fd_qpoptwin()

fd_minimaxoptwin()

fa_erroroptwin()

fireq.c

fd_cminimaxeq()

fd_cminimaxgdeq()

firfb.c

fd_halfband()

fd_dftfbadv()

fd_dftfb

ff_dftfb()

fd_pdftfb()

ff_pdftfb()

firtrans.c

ft_parallel()

ft_quantize()

iiranalog.c

id_bessel()

id_butterworth()

id_chebyshev()

id_ichebyshev()

id_elliptic()

iirtrans.c

it_lp2lp()

it_lp2hp()

it_lp2bp()

it_lp2bs()

it_bilinear()

it_bilinear_adv()

it_sos2poly()

it_sos2allpass()

it_quant_allpass()

it_allpass2poly()

matrix.h/matrix.c

mat_calloc()

mat_free()