Libfilth使用說明
winshton
2009年2月
（*本文大部分翻譯自libfilth,還有一部分是個人使用實踐
*時間水平均有限，翻譯的不完整，尤其第二章可以忽略）
版本歷史修改記錄
| 版本 |
作者 |
日期 |
備注 |
| V1.0 |
winshton |
2009-2-1 |
創建 |
-
目 錄
2.8.2. fd_minphaseminimax() 16
2.16.30. mat_clinspace_pol() 47
-
概述
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 [3*n/2] if n is even [3*n/2-1] if n is odd
den
Denominator filter taps [3*n/2] if n is even [3*n/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.
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 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
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 濾波器，其描述如下：
-
Bessel-Thomson濾波器在通帶上有着近乎恆定的群時延，使它們很適合於音頻應用。恆定的群時延犧牲的是一個較寬的過渡帶。濾波器的頻率響應只由N控制。選擇合適的長度沒有簡單的方法，但是可以通過下面表格的公式選取相近值。
-
Butterworth濾波器在ω = 0和ω = ∞最平滑，並且有着非常平滑的頻率和相位響應。濾波器的頻率響應僅由N決定，近似的取值為：
(4-1) 濾波器的頻率響應介於unity之間，0 < ω < ωp，ωp < 1 且 0 < G < 1。
-
Chebyshev通過參數ε控制通帶中的最大紋波。頻率響應在ω = ∞最平滑。參數ε可以通過下面的公式獲得：
(4-2) δ1是通帶紋波的期望值，a是通帶紋波信噪比，濾波器的序列數可有下面的公式獲得：
(4-3) 濾波器的頻率響應介於unity之間，0 < ω < ωs，ωs < 1 且 0 < G < 1。
-
inverse Chebyshev通過參數k1實現高阻帶衰減，頻率響應在ω = ∞最平滑。參數k1可以通過下面的公式獲得：
(4-4) δ2是阻帶紋波的期望值，b是阻帶衰減的信噪比，濾波器的序列數可有下面的公式獲得：
(4-5) 濾波器的頻率響應介於unity之間，0 < ω < ωp，ωp < 1 且 0 < G < 1。
-
Elliptic Function 濾波器特點由一下四個參數確定：
-
通帶紋波δ1
-
過渡帶寬度（1-ωs）
-
最大阻帶紋波δ2
-
濾波器階數N
-
給定上面的三個參數，第四個將會很小。設計方法由參數ε 和 k1確定。參數的計算公式如下：
(4-6) a是通帶紋波信噪比，b是阻帶衰減的信噪比，對於給定采樣次數N的濾波器，使用上述方程過渡區將會很小。
上述設計思想由iiranalog.cpp文件中的函數：id_bessel() ，id_butterworth() ，id_chebyshev() ，id_ichebyshev()，和id_elliptic() 實現。這些功能函數所實現的濾波器被稱為原型濾波器。擁有1Hz的截止頻率和單位增益的通帶。
只有Besel-Thomson和Butterworth濾波在通帶上有單位增益，其它三種濾波器的通帶增益由不同的參數決定。增益漂移由設計函數進行補償。
模擬濾波器變換
使用頻率變換可以改變原型濾波器的截止頻率和響應。下面是四個基本的變換：
-
低通到低通變換改變原型濾波器的截止頻率為F。
(4-7) -
低通到高通變換把低通原型濾波器轉換為截止頻率為F的高通濾波器。
(4-8) -
低通到帶通變換把低通原型濾波器轉換為截止頻率為Fl和Fh的帶通濾波器。
(4-9) -
低通到帶阻變換把低通原型濾波器轉換為截止頻率為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濾波器使用下面的步驟變換：
-
求多項式Q(z)。
由R(z)
(4-18) R(z)逆Z變換：
(4-19) 反過來給我們遞歸函數
(4-20) -
求兩個全通濾波器A1(z) 和 A2(z)。
(4-21) -
量化濾波器系數
-