# Design of IIR Digital Filters with Arbitrary Flatness Using Iterative Quadratic Programming

Authors: Yasunori Sugita
Journal: Journal of Electrical and Computer Engineering
 (2015)

---

## Abstract

This paper presents a design method of Chebyshev-type and inverse-Chebyshev-type infinite impulse response (IIR) filters with an approximately linear phase response. In the design of Chebyshev-type filters, the flatness condition in the stopband is preincorporated into a transfer function, and an equiripple characteristic in the passband is achieved by iteratively solving the QP problem using the transfer function. In the design of inverse-Chebyshev-type filters, the flatness condition in the passband is added to the constraint of the QP problem as the linear matrix equality, and an equiripple characteristic in the stopband is realized by iteratively solving the QP problem. To guarantee the stability of the obtained filters, we apply the extended positive realness to the QP problem. As a result, the proposed method can design the filters with more high precision than the conventional methods.
The effectiveness of the proposed design method is illustrated with some examples.


---

## Body

## 1. Introduction

Digital filters are needed in many applications in signal processing. Infinite impulse response (IIR) filters are useful for high-speed processing, and IIR filters with lower order can be realized that are comparable to the magnitude responses of finite impulse response (FIR) filters. However, its design is more difficult than FIR filters because IIR filters have a rational transfer function and are not always stable. For that reason, the study on IIR filter design is still a hot topic in the area of digital signal processing and many design methods on the IIR filters have been presented [1–8]. In this paper, we will treat two types of filters: Chebyshev type which is equiripple in the passband and flat in the stopband, Inverse Chebyshev type which is equiripple in the stopband and flat in the stopband. These filters are needed in many image processing applications to suppress ringing and chessboard distortion and are used in each situation appropriately. Moreover, in many signal processing applications, such as wavelet transform and channel equalization, the complex coefficient digital filters are required to meet some specifications that cannot be achieved by real coefficient filters, such as asymmetric spectral response [9–11].In [12, 13], the design methods based on the Remez algorithm were proposed for the Chebyshev-type and inverse-Chebyshev-type IIR filters with an approximately linear phase response. These methods can design the filters with small computational complexity. However, the filters that can be designed using these methods are restricted greatly because of a condition imposed on setting the initial value. Moreover, these methods cannot guarantee the stability of the filter obtained. By using the linear semi-infinite programming and the extended positive realness, the design method of stable inverse Chebyshev-type IIR filters with an approximately linear phase response has been proposed [14]. This method needs to set an initial value appropriately and the performance of the filter obtained is dependent on the initial value.In this paper, a design method based on quadratic programming (QP) is presented for stable IIR filters and for FIR filters with prescribed flatness in the passband or stopband and an approximately linear phase response in the passband. To guarantee the stability of the filter obtained, we apply the extended positive realness to the QP problem. As a result, the proposed method can design the filters that cannot be designed by the conventional methods. This paper is organized as follows: in Section2, the flatness conditions in the passband and stopband are described. In Section 3, the design algorithm for inverse-Chebyshev-type filters is described. In Section 4, the design algorithm for Chebyshev-type filters is described. To verify the effectiveness of the proposed method, several design examples are given in Section 5. Section 6 is the conclusions of this work.


## 2. IIR Digital Filters and Flatness Conditions

The frequency responseH
(
e
j
ω
) of an IIR digital filter is defined as
(1)
H
e
j
ω
=
A
(
e
j
ω
)
B
e
j
ω
=
∑
n
=
0
N
‍
a
n
e
-
j
n
ω
∑
m
=
0
M
‍
b
m
e
-
j
m
ω
,
where N and M are the orders of the numerator and denominator, respectively.  a
n
=
a
r
n
+
j
a
i
n and b
m
=
b
r
m
+
j
b
i
m are the filter coefficients and b
0
=
1 in general. The desired frequency response H
d
(
e
j
ω
) can be expressed as
(2)
H
d
(
e
j
ω
)
=
e
-
j
τ
d
ω
,
passband
,
0
,
stopband
,
where τ
d is a desired group delay.Then, the flatness conditions of the magnitude and group delay atω
=
ω
u in the passband are given as follows:
(3)
∂
i
H
e
j
ω
∂
ω
i
ω
=
ω
u
=
1
,
i
=
0
,
0
,
i
=
1,2
,
…
,
K
p
-
1
,
∂
i
τ
(
ω
)
∂
ω
i
ω
=
ω
u
=
τ
d
,
i
=
0
,
0
,
i
=
1,2
,
…
,
K
p
-
2
,
where K
p is a parameter expressing the flatness in the passband. The magnitude flatness condition at ω
=
ω
t in the stopband is
(4)
∂
i
H
e
j
ω
∂
ω
i
ω
=
ω
t
=
0
,
i
=
0,1
,
…
,
K
s
-
1
,
where K
s is a parameter expressing the flatness in the stopband.


## 3. Design of Inverse-Chebyshev-Type Filters

### 3.1. Problem Formulation

LetH
~
(
e
j
ω
) be a noncausal shifted version of H
(
e
j
ω
);
(5)
H
~
e
j
ω
=
H
(
e
j
ω
)
e
j
τ
d
ω
=
∑
n
=
0
N
‍
a
n
e
-
j
n
-
τ
d
ω
∑
m
=
0
M
‍
b
m
e
-
j
m
ω
=
A
~
e
j
ω
B
e
j
ω
.
With (5), the flatness conditions in (3) become
(6)
∂
i
H
~
e
j
ω
∂
ω
i
ω
=
ω
u
=
1
,
i
=
0
,
0
,
i
=
1,2
,
…
,
K
p
-
1
.
It is shown in Appendix that (6) is equivalent to (3). Equation (6) is the same as follows [15]:
(7)
∂
i
A
~
e
j
ω
∂
ω
i
ω
=
ω
u
=
∂
i
B
e
j
ω
∂
ω
i
ω
=
ω
u
,
(
i
=
0,1
,
…
,
K
p
-
1
)
.
Consequently, we can get the linear equation in matrix form:
(8)
U
h
=
V
,
where U
=
[
U
1
U
2
], V
=
[
1,0
,
…
,
0
]
T, and h
=
[
a
r
0
,
…
,
a
r
N
,
a
i
0
,
…
,
a
i
N
,
b
r
1
,
…
,
b
r
M
,
b
i
1
,
…
,
b
i
M
]
T and (
·
)
T denotes the transpose of (
·
). U
1 and U
2 are in (9) and (10), respectively:(9)
U
1
=
0
-
τ
d
0
cos
⁡
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
0
cos
⁡
⁡
(
(
N
-
τ
d
)
ω
u
)
0
-
τ
d
0
sin
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
0
sin
⁡
(
(
N
-
τ
d
)
ω
u
)
⋮
⋱
⋮
⋮
⋱
⋮
0
-
τ
d
K
p
-
1
cos
⁡
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
K
p
-
1
cos
⁡
⁡
(
(
N
-
τ
d
)
ω
u
)
0
-
τ
d
K
p
-
1
sin
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
K
p
-
1
sin
⁡
(
(
N
-
τ
d
)
ω
u
)
-
0
-
τ
d
0
sin
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
-
N
-
τ
d
0
sin
⁡
(
(
N
-
τ
d
)
ω
u
)
0
-
τ
d
0
cos
⁡
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
0
cos
⁡
⁡
(
(
N
-
τ
d
)
ω
u
)
⋮
⋱
⋮
⋮
⋱
⋮
-
0
-
τ
d
K
p
-
1
sin
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
-
N
-
τ
d
K
p
-
1
sin
⁡
(
(
N
-
τ
d
)
ω
u
)
0
-
τ
d
K
p
-
1
cos
⁡
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
K
p
-
1
cos
⁡
⁡
(
(
N
-
τ
d
)
ω
u
)
,
(10)
U
2
=
-
1
0
cos
⁡
(
1
ω
u
)
⋯
-
M
0
cos
⁡
(
M
ω
u
)
-
1
0
sin
(
1
ω
u
)
⋯
-
M
0
sin
(
M
ω
u
)
⋮
⋱
⋮
⋮
⋱
⋮
-
1
K
p
-
1
cos
⁡
(
1
ω
u
)
⋯
-
M
K
p
-
1
cos
⁡
(
M
ω
u
)
-
1
K
p
-
1
sin
(
1
ω
u
)
⋯
-
M
K
p
-
1
sin
(
M
ω
u
)
1
0
sin
(
1
ω
u
)
⋯
M
0
sin
(
M
ω
u
)
-
1
0
cos
⁡
(
1
ω
u
)
⋯
-
M
0
cos
⁡
(
M
ω
u
)
⋮
⋱
⋮
⋮
⋱
⋮
1
K
p
-
1
sin
(
1
ω
u
)
⋯
M
K
p
-
1
sin
(
M
ω
u
)
-
1
K
p
-
1
cos
⁡
(
1
ω
u
)
⋯
-
M
K
p
-
1
cos
⁡
(
M
ω
u
)
.LetH
~
d
(
e
j
ω
) be the desired magnitude response; that is,
(11)
H
~
d
(
e
j
ω
)
=
1
,
passband
,
0
,
stopband
.
Using (5) and (11), the weighted least squares design problem is
(12)
min
⁡
h
⁡
J
=
∑
l
=
1
L
‍
W
ω
l
H
~
e
j
ω
l
-
H
~
d
e
j
ω
l
2
,
where L is the total number of grid points in the passband and stopband, W
(
ω
l
) is the weighting function, and ω
l  (
l
=
1
,
…
,
L
) are the discrete frequency points used in the calculation. However, it is difficult to solve (12) directly because H
~
(
e
j
ω
l
) is a rational function. Thus, we use the following iterative design formula:
(13)
min
⁡
h
J
=
∑
l
=
1
L
‍
W
ω
l
A
~
e
j
ω
l
-
H
~
d
e
j
ω
l
B
e
j
ω
l
2
B
k
-
1
e
j
ω
l
2
,
where k is the number of the iterations.After some manipulation, (13) can be formulated as the following QP problem:
(14)
min
⁡
h
⁡
h
T
P
R
T
W
P
R
+
P
I
T
W
P
I
h
-
2
d
R
T
W
P
R
-
d
I
T
W
P
I
h
,
where P
R
=
[
(
e
N
R
+
e
N
I
)
,
-
(
D
R
e
M
R
+
D
I
e
M
I
)
], P
I
=
[
(
e
N
I
-
e
N
R
)
,
(
D
I
e
M
R
-
D
R
e
M
I
)
], e
N
R
=
[
1
,
…
,
cos
⁡
⁡
(
ω
l
(
N
-
τ
d
)
)
,
0
,
…
,
sin
⁡
(
ω
l
(
N
-
τ
d
)
)
], e
N
I
=
[
0
,
…
,
-
sin
⁡
(
ω
l
(
N
-
τ
d
)
)
,
1
,
…
,
cos
⁡
⁡
(
ω
l
(
N
-
τ
d
)
)
], e
M
R
=
[
cos
⁡
⁡
(
ω
l
)
,
…
,
cos
⁡
⁡
(
ω
l
M
)
,
sin
⁡
(
ω
l
)
,
…
,
sin
⁡
(
ω
l
M
)
], e
M
I
=
[
-
sin
⁡
(
ω
l
)
,
…
,
-
sin
⁡
(
ω
l
M
)
,
cos
⁡
⁡
(
ω
l
)
,
…
,
cos
⁡
⁡
(
ω
l
M
)
], D
=
diag
(
d
)
=
D
R
-
j
D
I, d
=
[
H
~
d
(
e
j
ω
1
)
,
…
,
H
~
d
(
e
j
ω
L
)
]
=
d
R
-
j
d
I, and W
=
diag
W
(
ω
1
)
/
|
B
k
-
1
(
ω
0
)
|
2
,
…
,
W
(
ω
L
)
/
|
B
k
-
1
(
ω
L
)
|
2.


### 3.2. Update of the Weighting FunctionW
(
ω
)

It has been well known that the filters obtained under weighted least square criterion have a large ripple near the band edges. So in order to realize the equiripple characteristics in the passband or stopband or both, the weighting function used at every iteration is adjusted using the modified Lawson’s method [16] and the QP problem is solved to obtain the coefficients. In this paper, the weighting function W
(
ω
) in kth iteration step is updated as follows:
(15)
W
k
+
1
ω
=
W
k
ω
β
k
ω
1
/
L
∑
l
=
1
L
‍
W
k
ω
l
β
k
ω
l
,
where the envelope function β
k
(
ω
) is given as the function of straight line formed by joining together all the extremal points of the same frequency band of interest on the error function which is expressed as
(16)
E
k
ω
=
H
~
e
j
ω
-
H
~
d
e
j
ω
.
Using the extremal points ω
^
i of E
k
(
ω
), β
k
(
ω
) can be calculated by
(17)
β
k
ω
=
ω
-
ω
^
i
ω
^
i
+
1
-
ω
^
i
E
k
ω
^
i
+
1
+
ω
^
i
+
1
-
ω
ω
^
i
+
1
-
ω
^
i
E
k
ω
^
i
,
for
ω
^
i
<
ω
<
ω
^
i
+
1
,
where ω
^
i denotes the ith extremal frequency of the error function E
k
(
ω
).


### 3.3. Stability Constraint

To obtain the stable IIR filters, the stability condition based on a positive realness has been applied to many design methods. However, the use of this condition may exclude the candidate for the transfer function with excellent performance because this condition is a sufficient condition to assure the stability and is often too restrictive. In [17], an iterative method for the stability guarantee based on the positive realness was proposed in order to get a better transfer function. In this method, a stability condition is given by
(18)
Re
B
e
j
ω
≥
δ
,
-
π
≤
ω
<
π
,
where δ
<
1. If the maximum pole radius P
max
⁡ of the filter obtained using a given δ is greater than prescribed maximum allowable pole radius r
m, δ is updated and redesign is carried out using the updated δ. The update of δ and the redesign are repeated until satisfying |
P
max
⁡
-
r
m
|
≤
ϵ
R, where ϵ
R is a positive small value. The update procedure of δ is described in Section 4.2.Using the discrete angular frequencyω
i  (
i
=
1
,
…
,
R
), (18) can be expressed as the linear matrix inequality
(19)
Γ
h
≥
λ
,
where
(20)
Γ
=
0
⋯
0
cos
⁡
⁡
(
ω
1
)
⋯
cos
⁡
⁡
(
M
ω
1
)
⋮
⋱
⋮
⋮
⋱
⋮
0
⋯
0
cos
⁡
⁡
(
ω
R
)
⋯
cos
⁡
⁡
(
M
ω
R
)
,
λ
=
δ
-
1
,
…
,
δ
-
1
T
.Thus, the design problem in which the stability constraint and the update of the weighting function were considered becomes a standard QP problem as below:(21)
min
⁡
h
k
h
k
T
P
R
T
W
k
P
R
+
P
I
T
W
k
P
I
h
k
-
2
d
R
T
W
k
P
R
-
d
I
T
W
k
P
I
h
k
,
sub
.
to
Γ
h
k
≥
λ
,
U
h
k
=
V
,
where W
k
=
diag
(
[
W
k
(
ω
1
)
/
|
B
k
-
1
(
ω
0
)
|
2
,
…
,
W
k
(
ω
L
)
/
|
B
k
-
1
ω
L
|
2
]
).This problem can be solved using a powerful QP tool, such asq
u
a
d
p
r
o
g in MATLAB.


## 3.1. Problem Formulation

LetH
~
(
e
j
ω
) be a noncausal shifted version of H
(
e
j
ω
);
(5)
H
~
e
j
ω
=
H
(
e
j
ω
)
e
j
τ
d
ω
=
∑
n
=
0
N
‍
a
n
e
-
j
n
-
τ
d
ω
∑
m
=
0
M
‍
b
m
e
-
j
m
ω
=
A
~
e
j
ω
B
e
j
ω
.
With (5), the flatness conditions in (3) become
(6)
∂
i
H
~
e
j
ω
∂
ω
i
ω
=
ω
u
=
1
,
i
=
0
,
0
,
i
=
1,2
,
…
,
K
p
-
1
.
It is shown in Appendix that (6) is equivalent to (3). Equation (6) is the same as follows [15]:
(7)
∂
i
A
~
e
j
ω
∂
ω
i
ω
=
ω
u
=
∂
i
B
e
j
ω
∂
ω
i
ω
=
ω
u
,
(
i
=
0,1
,
…
,
K
p
-
1
)
.
Consequently, we can get the linear equation in matrix form:
(8)
U
h
=
V
,
where U
=
[
U
1
U
2
], V
=
[
1,0
,
…
,
0
]
T, and h
=
[
a
r
0
,
…
,
a
r
N
,
a
i
0
,
…
,
a
i
N
,
b
r
1
,
…
,
b
r
M
,
b
i
1
,
…
,
b
i
M
]
T and (
·
)
T denotes the transpose of (
·
). U
1 and U
2 are in (9) and (10), respectively:(9)
U
1
=
0
-
τ
d
0
cos
⁡
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
0
cos
⁡
⁡
(
(
N
-
τ
d
)
ω
u
)
0
-
τ
d
0
sin
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
0
sin
⁡
(
(
N
-
τ
d
)
ω
u
)
⋮
⋱
⋮
⋮
⋱
⋮
0
-
τ
d
K
p
-
1
cos
⁡
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
K
p
-
1
cos
⁡
⁡
(
(
N
-
τ
d
)
ω
u
)
0
-
τ
d
K
p
-
1
sin
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
K
p
-
1
sin
⁡
(
(
N
-
τ
d
)
ω
u
)
-
0
-
τ
d
0
sin
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
-
N
-
τ
d
0
sin
⁡
(
(
N
-
τ
d
)
ω
u
)
0
-
τ
d
0
cos
⁡
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
0
cos
⁡
⁡
(
(
N
-
τ
d
)
ω
u
)
⋮
⋱
⋮
⋮
⋱
⋮
-
0
-
τ
d
K
p
-
1
sin
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
-
N
-
τ
d
K
p
-
1
sin
⁡
(
(
N
-
τ
d
)
ω
u
)
0
-
τ
d
K
p
-
1
cos
⁡
⁡
(
(
0
-
τ
d
)
ω
u
)
⋯
N
-
τ
d
K
p
-
1
cos
⁡
⁡
(
(
N
-
τ
d
)
ω
u
)
,
(10)
U
2
=
-
1
0
cos
⁡
(
1
ω
u
)
⋯
-
M
0
cos
⁡
(
M
ω
u
)
-
1
0
sin
(
1
ω
u
)
⋯
-
M
0
sin
(
M
ω
u
)
⋮
⋱
⋮
⋮
⋱
⋮
-
1
K
p
-
1
cos
⁡
(
1
ω
u
)
⋯
-
M
K
p
-
1
cos
⁡
(
M
ω
u
)
-
1
K
p
-
1
sin
(
1
ω
u
)
⋯
-
M
K
p
-
1
sin
(
M
ω
u
)
1
0
sin
(
1
ω
u
)
⋯
M
0
sin
(
M
ω
u
)
-
1
0
cos
⁡
(
1
ω
u
)
⋯
-
M
0
cos
⁡
(
M
ω
u
)
⋮
⋱
⋮
⋮
⋱
⋮
1
K
p
-
1
sin
(
1
ω
u
)
⋯
M
K
p
-
1
sin
(
M
ω
u
)
-
1
K
p
-
1
cos
⁡
(
1
ω
u
)
⋯
-
M
K
p
-
1
cos
⁡
(
M
ω
u
)
.LetH
~
d
(
e
j
ω
) be the desired magnitude response; that is,
(11)
H
~
d
(
e
j
ω
)
=
1
,
passband
,
0
,
stopband
.
Using (5) and (11), the weighted least squares design problem is
(12)
min
⁡
h
⁡
J
=
∑
l
=
1
L
‍
W
ω
l
H
~
e
j
ω
l
-
H
~
d
e
j
ω
l
2
,
where L is the total number of grid points in the passband and stopband, W
(
ω
l
) is the weighting function, and ω
l  (
l
=
1
,
…
,
L
) are the discrete frequency points used in the calculation. However, it is difficult to solve (12) directly because H
~
(
e
j
ω
l
) is a rational function. Thus, we use the following iterative design formula:
(13)
min
⁡
h
J
=
∑
l
=
1
L
‍
W
ω
l
A
~
e
j
ω
l
-
H
~
d
e
j
ω
l
B
e
j
ω
l
2
B
k
-
1
e
j
ω
l
2
,
where k is the number of the iterations.After some manipulation, (13) can be formulated as the following QP problem:
(14)
min
⁡
h
⁡
h
T
P
R
T
W
P
R
+
P
I
T
W
P
I
h
-
2
d
R
T
W
P
R
-
d
I
T
W
P
I
h
,
where P
R
=
[
(
e
N
R
+
e
N
I
)
,
-
(
D
R
e
M
R
+
D
I
e
M
I
)
], P
I
=
[
(
e
N
I
-
e
N
R
)
,
(
D
I
e
M
R
-
D
R
e
M
I
)
], e
N
R
=
[
1
,
…
,
cos
⁡
⁡
(
ω
l
(
N
-
τ
d
)
)
,
0
,
…
,
sin
⁡
(
ω
l
(
N
-
τ
d
)
)
], e
N
I
=
[
0
,
…
,
-
sin
⁡
(
ω
l
(
N
-
τ
d
)
)
,
1
,
…
,
cos
⁡
⁡
(
ω
l
(
N
-
τ
d
)
)
], e
M
R
=
[
cos
⁡
⁡
(
ω
l
)
,
…
,
cos
⁡
⁡
(
ω
l
M
)
,
sin
⁡
(
ω
l
)
,
…
,
sin
⁡
(
ω
l
M
)
], e
M
I
=
[
-
sin
⁡
(
ω
l
)
,
…
,
-
sin
⁡
(
ω
l
M
)
,
cos
⁡
⁡
(
ω
l
)
,
…
,
cos
⁡
⁡
(
ω
l
M
)
], D
=
diag
(
d
)
=
D
R
-
j
D
I, d
=
[
H
~
d
(
e
j
ω
1
)
,
…
,
H
~
d
(
e
j
ω
L
)
]
=
d
R
-
j
d
I, and W
=
diag
W
(
ω
1
)
/
|
B
k
-
1
(
ω
0
)
|
2
,
…
,
W
(
ω
L
)
/
|
B
k
-
1
(
ω
L
)
|
2.


## 3.2. Update of the Weighting FunctionW
(
ω
)

It has been well known that the filters obtained under weighted least square criterion have a large ripple near the band edges. So in order to realize the equiripple characteristics in the passband or stopband or both, the weighting function used at every iteration is adjusted using the modified Lawson’s method [16] and the QP problem is solved to obtain the coefficients. In this paper, the weighting function W
(
ω
) in kth iteration step is updated as follows:
(15)
W
k
+
1
ω
=
W
k
ω
β
k
ω
1
/
L
∑
l
=
1
L
‍
W
k
ω
l
β
k
ω
l
,
where the envelope function β
k
(
ω
) is given as the function of straight line formed by joining together all the extremal points of the same frequency band of interest on the error function which is expressed as
(16)
E
k
ω
=
H
~
e
j
ω
-
H
~
d
e
j
ω
.
Using the extremal points ω
^
i of E
k
(
ω
), β
k
(
ω
) can be calculated by
(17)
β
k
ω
=
ω
-
ω
^
i
ω
^
i
+
1
-
ω
^
i
E
k
ω
^
i
+
1
+
ω
^
i
+
1
-
ω
ω
^
i
+
1
-
ω
^
i
E
k
ω
^
i
,
for
ω
^
i
<
ω
<
ω
^
i
+
1
,
where ω
^
i denotes the ith extremal frequency of the error function E
k
(
ω
).


## 3.3. Stability Constraint

To obtain the stable IIR filters, the stability condition based on a positive realness has been applied to many design methods. However, the use of this condition may exclude the candidate for the transfer function with excellent performance because this condition is a sufficient condition to assure the stability and is often too restrictive. In [17], an iterative method for the stability guarantee based on the positive realness was proposed in order to get a better transfer function. In this method, a stability condition is given by
(18)
Re
B
e
j
ω
≥
δ
,
-
π
≤
ω
<
π
,
where δ
<
1. If the maximum pole radius P
max
⁡ of the filter obtained using a given δ is greater than prescribed maximum allowable pole radius r
m, δ is updated and redesign is carried out using the updated δ. The update of δ and the redesign are repeated until satisfying |
P
max
⁡
-
r
m
|
≤
ϵ
R, where ϵ
R is a positive small value. The update procedure of δ is described in Section 4.2.Using the discrete angular frequencyω
i  (
i
=
1
,
…
,
R
), (18) can be expressed as the linear matrix inequality
(19)
Γ
h
≥
λ
,
where
(20)
Γ
=
0
⋯
0
cos
⁡
⁡
(
ω
1
)
⋯
cos
⁡
⁡
(
M
ω
1
)
⋮
⋱
⋮
⋮
⋱
⋮
0
⋯
0
cos
⁡
⁡
(
ω
R
)
⋯
cos
⁡
⁡
(
M
ω
R
)
,
λ
=
δ
-
1
,
…
,
δ
-
1
T
.Thus, the design problem in which the stability constraint and the update of the weighting function were considered becomes a standard QP problem as below:(21)
min
⁡
h
k
h
k
T
P
R
T
W
k
P
R
+
P
I
T
W
k
P
I
h
k
-
2
d
R
T
W
k
P
R
-
d
I
T
W
k
P
I
h
k
,
sub
.
to
Γ
h
k
≥
λ
,
U
h
k
=
V
,
where W
k
=
diag
(
[
W
k
(
ω
1
)
/
|
B
k
-
1
(
ω
0
)
|
2
,
…
,
W
k
(
ω
L
)
/
|
B
k
-
1
ω
L
|
2
]
).This problem can be solved using a powerful QP tool, such asq
u
a
d
p
r
o
g in MATLAB.


## 4. Design of Chebyshev-Type Filters

### 4.1. Problem Formulation

In this section, we consider the filters with a flatness in the stopband. With a noncausal shifted versionH
~
(
e
j
ω
) in (5), the flatness condition in (4) becomes
(22)
∂
i
H
~
e
j
ω
∂
ω
i
ω
=
ω
t
=
0
,
(
i
=
0,1
,
…
,
K
s
-
1
)
.
In order to meet the flatness condition of (22), it is necessary to place K
s multiple zeros at ω
=
ω
t. Hence, the frequency response H
^
(
e
j
ω
) can be expressed as
(23)
H
^
e
j
ω
=
A
^
e
j
ω
B
e
j
ω
=
1
+
e
-
j
ω
K
s
∑
n
=
0
N
-
K
s
‍
c
n
e
-
j
(
n
-
τ
d
)
ω
1
+
∑
m
=
1
M
‍
b
m
e
-
j
m
ω
,
where c
n
=
c
r
n
+
j
c
i
n. To obtain the equiripple responses in the passband, we consider the following iterative design formula:
(24)
min
⁡
c
,
b
J
=
∑
l
=
1
L
‍
W
(
ω
l
)
B
k
-
1
e
j
ω
l
2
A
^
e
j
ω
l
-
H
d
e
j
ω
l
B
e
j
ω
l
2
.
Thus, the design problem in which the stability constraint and the update of the weighting function were considered becomes a standard QP problem as below:
(25)
min
⁡
h
^
k
h
^
k
T
P
^
R
T
W
k
P
^
R
+
P
^
I
T
W
k
P
^
I
h
^
k
-
2
d
R
T
W
k
P
^
R
-
d
I
T
W
k
P
^
I
h
^
k
,
sub
.
to
Γ
h
^
k
≥
λ
,
where h
^
=
[
c
r
0
,
…
,
c
r
(
N
-
K
s
)
,
c
i
0
,
…
,
c
i
(
N
-
K
s
)
,
b
r
1
,
…
,
b
r
M,b
i
1
,
…
,
b
i
M
]
T, P
^
R
=
[
(
G
R
e
^
N
R
+
G
I
e
^
N
I
)
,
-
(
D
R
e
M
R
+
D
I
e
M
I
)
], P
^
I
=
[
(
G
R
e
^
N
I
-
G
I
e
^
N
R
)
,
(
D
I
e
M
R
-
D
R
e
M
I
)
], G
=
diag
(
(
1
-
e
-
j
(
ω
l
-
ω
t
)
)
K
s
)
=
G
R
-
j
G
I, e
^
N
R
=
[
1
,
…
,
cos
⁡
(
ω
l
(
N
-
K
s
)
)
,
0
,
…
,
sin
(
ω
l
(
N
-
K
s
)
)
], e
^
N
I
=
[
0
,
…
,
-
sin
(
ω
l
(
N
-
K
s
)
)
,
1
,
…
,
cos
⁡
(
ω
l
(
N
-
K
s
)
)
], e
M
R
=
[
cos
⁡
(
ω
l
)
,
…
,
cos
⁡
(
ω
l
M
)
,
sin
(
ω
l
)
,
…
,
sin
(
ω
l
M
)
], e
M
I
=
[
-
sin
(
ω
l
)
,
…
,
-
sin
(
ω
l
M
)
,
cos
⁡
(
ω
l
)
,
…
,
cos
⁡
(
ω
l
M
)
], D
=
diag
(
d
)
=
D
R
-
j
D
I, d
=
[
H
d
(
e
j
ω
1
)
,
…
,
H
d
(
e
j
ω
L
)
]
=
d
R
-
j
d
I, and W
k
=
diag
W
k
(
ω
1
)
/
|
B
k
-
1
(
ω
1
)
|
2
,
…
,
W
k
(
ω
L
)
/
|
B
k
-
1
(
ω
L
)
|
2.


### 4.2. Design Procedure

The design procedure of the proposed method is summarized as follows.Step 0. Set the design specifications.Step 1. Solve the QP problem to obtain the filter coefficient h
k.Step 2. If sum
(
|
h
k
-
h
k
-
1
|
)
/
sum
(
|
h
k
|
)
≤
ϵ and P
max
⁡
<
r
m, stop; if sum
(
|
h
k
-
h
k
-
1
|
)
/
sum
(
|
h
k
|
)
≤
ϵ and P
max
⁡
>
r
m, go to Step 4; otherwise, go to Step 3.Step 3. Update the weighting function W
(
ω
) using (15)–(17) and go back to Step 1.Step 4. Calculate δ
′
=
min
⁡
all
ω
Re
{
B
(
ω
)
} from the obtained filter and then set to δ
l
=
δ
′ and δ
u
=
1.Step 5. Calculate δ
=
(
δ
l
+
δ
u
)
/
2.Step 6. Solve the QP problem using the updated δ.Step 7. If sum
(
|
h
k
-
h
k
-
1
|
)
/
sum
(
|
h
k
|
)
≤
ϵ and |
P
max
⁡
-
r
m
|
≤
ϵ
R, stop; if sum
(
|
h
k
-
h
k
-
1
|
)
/
sum
(
|
h
k
|
)
≤
ϵ and |
P
max
⁡
-
r
m
|
>
ϵ
R, go to Step 9; otherwise then go to Step 8.Step 8. Update the weighting function W
(
ω
) using (15)–(17) and go back to Step 6.Step 9. If P
max
⁡
<
r
m, set to δ
u
=
δ; if P
max
⁡
>
r
m, set to δ
l
=
δ, and then go back to Step 5.


## 4.1. Problem Formulation

In this section, we consider the filters with a flatness in the stopband. With a noncausal shifted versionH
~
(
e
j
ω
) in (5), the flatness condition in (4) becomes
(22)
∂
i
H
~
e
j
ω
∂
ω
i
ω
=
ω
t
=
0
,
(
i
=
0,1
,
…
,
K
s
-
1
)
.
In order to meet the flatness condition of (22), it is necessary to place K
s multiple zeros at ω
=
ω
t. Hence, the frequency response H
^
(
e
j
ω
) can be expressed as
(23)
H
^
e
j
ω
=
A
^
e
j
ω
B
e
j
ω
=
1
+
e
-
j
ω
K
s
∑
n
=
0
N
-
K
s
‍
c
n
e
-
j
(
n
-
τ
d
)
ω
1
+
∑
m
=
1
M
‍
b
m
e
-
j
m
ω
,
where c
n
=
c
r
n
+
j
c
i
n. To obtain the equiripple responses in the passband, we consider the following iterative design formula:
(24)
min
⁡
c
,
b
J
=
∑
l
=
1
L
‍
W
(
ω
l
)
B
k
-
1
e
j
ω
l
2
A
^
e
j
ω
l
-
H
d
e
j
ω
l
B
e
j
ω
l
2
.
Thus, the design problem in which the stability constraint and the update of the weighting function were considered becomes a standard QP problem as below:
(25)
min
⁡
h
^
k
h
^
k
T
P
^
R
T
W
k
P
^
R
+
P
^
I
T
W
k
P
^
I
h
^
k
-
2
d
R
T
W
k
P
^
R
-
d
I
T
W
k
P
^
I
h
^
k
,
sub
.
to
Γ
h
^
k
≥
λ
,
where h
^
=
[
c
r
0
,
…
,
c
r
(
N
-
K
s
)
,
c
i
0
,
…
,
c
i
(
N
-
K
s
)
,
b
r
1
,
…
,
b
r
M,b
i
1
,
…
,
b
i
M
]
T, P
^
R
=
[
(
G
R
e
^
N
R
+
G
I
e
^
N
I
)
,
-
(
D
R
e
M
R
+
D
I
e
M
I
)
], P
^
I
=
[
(
G
R
e
^
N
I
-
G
I
e
^
N
R
)
,
(
D
I
e
M
R
-
D
R
e
M
I
)
], G
=
diag
(
(
1
-
e
-
j
(
ω
l
-
ω
t
)
)
K
s
)
=
G
R
-
j
G
I, e
^
N
R
=
[
1
,
…
,
cos
⁡
(
ω
l
(
N
-
K
s
)
)
,
0
,
…
,
sin
(
ω
l
(
N
-
K
s
)
)
], e
^
N
I
=
[
0
,
…
,
-
sin
(
ω
l
(
N
-
K
s
)
)
,
1
,
…
,
cos
⁡
(
ω
l
(
N
-
K
s
)
)
], e
M
R
=
[
cos
⁡
(
ω
l
)
,
…
,
cos
⁡
(
ω
l
M
)
,
sin
(
ω
l
)
,
…
,
sin
(
ω
l
M
)
], e
M
I
=
[
-
sin
(
ω
l
)
,
…
,
-
sin
(
ω
l
M
)
,
cos
⁡
(
ω
l
)
,
…
,
cos
⁡
(
ω
l
M
)
], D
=
diag
(
d
)
=
D
R
-
j
D
I, d
=
[
H
d
(
e
j
ω
1
)
,
…
,
H
d
(
e
j
ω
L
)
]
=
d
R
-
j
d
I, and W
k
=
diag
W
k
(
ω
1
)
/
|
B
k
-
1
(
ω
1
)
|
2
,
…
,
W
k
(
ω
L
)
/
|
B
k
-
1
(
ω
L
)
|
2.


## 4.2. Design Procedure

The design procedure of the proposed method is summarized as follows.Step 0. Set the design specifications.Step 1. Solve the QP problem to obtain the filter coefficient h
k.Step 2. If sum
(
|
h
k
-
h
k
-
1
|
)
/
sum
(
|
h
k
|
)
≤
ϵ and P
max
⁡
<
r
m, stop; if sum
(
|
h
k
-
h
k
-
1
|
)
/
sum
(
|
h
k
|
)
≤
ϵ and P
max
⁡
>
r
m, go to Step 4; otherwise, go to Step 3.Step 3. Update the weighting function W
(
ω
) using (15)–(17) and go back to Step 1.Step 4. Calculate δ
′
=
min
⁡
all
ω
Re
{
B
(
ω
)
} from the obtained filter and then set to δ
l
=
δ
′ and δ
u
=
1.Step 5. Calculate δ
=
(
δ
l
+
δ
u
)
/
2.Step 6. Solve the QP problem using the updated δ.Step 7. If sum
(
|
h
k
-
h
k
-
1
|
)
/
sum
(
|
h
k
|
)
≤
ϵ and |
P
max
⁡
-
r
m
|
≤
ϵ
R, stop; if sum
(
|
h
k
-
h
k
-
1
|
)
/
sum
(
|
h
k
|
)
≤
ϵ and |
P
max
⁡
-
r
m
|
>
ϵ
R, go to Step 9; otherwise then go to Step 8.Step 8. Update the weighting function W
(
ω
) using (15)–(17) and go back to Step 6.Step 9. If P
max
⁡
<
r
m, set to δ
u
=
δ; if P
max
⁡
>
r
m, set to δ
l
=
δ, and then go back to Step 5.


## 5. Design Examples

In this section, some design examples are given to illustrate the effectiveness of the proposed method. In all the following examples,ϵ
=
1
0
-
7, ϵ
R
=
1
0
-
5, and the initial value of δ is -
1
0
2. The program of the proposed design algorithm was coded by MATLAB and was run on the PC which has Core i5-M520 processor and 2 GB memory.


### 5.1. Example 1

To compare with the conventional method [12], we design the Chebyshev-type IIR filters which have an equiripple characteristic in passband and a flat characteristic in stopband. The filter specifications are as follows:
(26)
N
=
15
,
M
=
6
,
K
s
=
11
,
H
d
(
e
j
ω
)
=
e
-
j
τ
d
ω
,
0
≤
ω
≤
0.3
π
,
0
,
ω
=
π
,
τ
d
=
9,10,11,12,13
.
In the proposed method, we set the maximum allowable pole radius r
m to 0.98, and the number of grid points is L
=
500 and R
=
1000. Note that the conventional method [12] cannot guarantee the stability of the obtained filter.The performance of the obtained filters is shown in Table1 which are obtained by the conventional method [12] based on Remez algorithm. In Table 1, R
p is the maximum magnitude error in passband, G
err
⁡ is the maximum group delay error in passband, and P
max
⁡ is the maximum pole radius of the obtained filter. From Table 1, it is confirmed that the proposed method needs more computational cost (CPU time and iteration number) than the conventional method, but the characteristic of the obtained filters by the proposed method is about the same as or better than them by the conventional method. Particularly, the filter obtained by the conventional method is unstable when τ
d
=
9. The stability of the filters obtained by the method which cannot guarantee the stability depends on a given group delay specification. In fact, when the group delay τ
d is less than 9.5, all filters obtained by conventional method were unstable. In contrast, the proposed method can always obtain the stable filters because the maximum pole radius of the filter obtained can be prespecified. Therefore, the proposed method can design the filters that cannot be designed by the conventional method.Table 1
Comparison with [12] in Example 1.
τ
d
R
p
G
err
⁡
P
max
⁡
CPU time (s)
Iteration number
Proposed
9
5.622 × 10−5
3.618 × 10−3
0.980
1.334
16
10
1.248 × 10−5
7.741 × 10−4
0.952
0.587
10
11
4.416 × 10−5
6.501 × 10−3
0.890
0.537
10
12
1.098 × 10−4
3.224 × 10−2
0.858
0.510
8
13
2.060 × 10−4
1.146 × 10−1
0.847
0.525
9
Reference  [12]
9
2.281 × 10−6
2.657 × 10−5
1.051
0.121
6
10
1.269 × 10−5
7.257 × 10−4
0.953
0.134
6
11
4.608 × 10−5
6.499 × 10−3
0.891
0.124
6
12
1.168 × 10−4
3.238 × 10−2
0.858
0.124
6
13
3.129 × 10−4
1.122 × 10−1
0.846
0.130
6


### 5.2. Example 2

Next we design the inverse Chebyshev-type IIR filters with the real coefficients in order to compare with the conventional method [13]. The filter specifications are as follows:
(27)
N
=
12
,
M
=
5
,
K
p
=
10
,
H
d
(
e
j
ω
)
=
e
-
j
τ
d
ω
,
ω
=
0
,
0
,
0.5
π
≤
ω
≤
π
,
τ
d
=
7,8
,
9,10,11
.
In the proposed method, we set the maximum allowable pole radius r
m to 0.98, and the number of grid points is L
=
500 and R
=
1000. Note that the conventional method [13] cannot guarantee the stability of the obtained filter.The numerical performance of the obtained filters is shown in Table2 which are obtained by the conventional method [13] based on Remez algorithm. In Table 2, R
s is the minimum stopband attenuation in dB. From Table 2, it is confirmed that the proposed method needs more computational cost than the conventional method but gives slightly better characteristics. Moreover, when the group delay τ
d is less than 7.2, the filters obtained by conventional method were unstable. In contrast, the proposed method can always obtain the stable filters. Therefore, the proposed method can design the filters that cannot be designed by the conventional method.Table 2
Comparison with [13] in Example 2.
τ
d
R
s (dB)
P
max
⁡
CPU time (s)
Iteration number
Proposed
7
27.85
0.980
1.812
45
8
35.02
0.928
0.512
16
9
41.38
0.850
0.372
14
10
46.61
0.784
0.385
12
11
50.89
0.746
0.351
11
Reference  [13]
7
25.78
1.002
0.212
7
8
34.07
0.924
0.205
7
9
40.46
0.850
0.202
7
10
45.76
0.784
0.207
7
11
50.10
0.746
0.198
6Next, we designed the filter withτ
d
=
9 for comparison with [14]. The minimum stopband attenuation of the filter obtained was 41.38dB in the proposed method and 40.04dB in [14]. It is thought that this difference is because the performance of the filter obtained by the conventional method is dependent on the given initial value.


### 5.3. Example 3

Next, as the example of the filter with the complex coefficients, we design the inverse Chebyshev-type IIR filter with the following specifications:(28)
N
=
12
,
M
=
5
,
H
d
(
e
j
ω
)
=
e
-
j
9
ω
,
ω
=
0.1
π
,
0
,
-
π
≤
ω
≤
-
0.6
π
,
0.7
π
≤
ω
<
π
,
K
p
=
9,10,11
.
In this example, we use r
m
=
1.0, L
=
500, and R
=
1000.The proposed method took an average of 0.78 seconds for convergence. The magnitude response and group delay response of the obtained filter are shown in Figures1, 2, and 3. From these figures, it is confirmed that the magnitude and group delay responses have a flat characteristic at ω
=
0.1
π and the magnitude response in stopband is equiripple. Table 3 is the numerical performance of the obtained filter.Table 3
The performance of the filters in Example 3.
K
p
R
s (dB)
P
max
⁡
9
94.43
0.902
10
69.92
0.865
11
48.03
0.856Figure 1
Overall magnitude response of the proposed IIR filters in Example 3.Figure 2
Passband magnitude response of the proposed IIR filters in Example 3.Figure 3
Group delay response of the proposed IIR filters in Example 3.


### 5.4. Example 4

Finally, we show the design examples of inverse-Chebyshev-type FIR filters. The filter specifications are as follows:(29)
N
=
1999
,
M
=
0
,
H
d
(
e
j
ω
)
=
e
-
j
900
ω
,
ω
=
0
,
0
,
0.0080
π
≤
ω
≤
π
,
K
p
=
10,12,14
.
In this example, we set the number of the grit points L to 10,000.The proposed method took an average of 2,794 seconds for convergence. The magnitude response and group delay response of the obtained filter are shown in Figures4, 5, and 6. From these figures, it is confirmed that the magnitude and group delay responses have a flat characteristic at ω
=
0 and the magnitude response in stopband is equiripple. The minimum stopband attenuation of the obtained filters was 72.98 (dB) for K
p
=
10, 59.11 (dB) for K
p
=
12, and 47.19  (dB) for K
p
=
14. We confirmed that these filters cannot be designed by the conventional method [13]. Hence, the proposed method can also design the high-order filters that cannot be designed by the conventional method.Figure 4
Overall magnitude response of the proposed FIR filters in Example 4.Figure 5
Passband magnitude response of the proposed FIR filters in Example 4.Figure 6
Group delay response of the proposed FIR filters in Example 4.


## 5.1. Example 1

To compare with the conventional method [12], we design the Chebyshev-type IIR filters which have an equiripple characteristic in passband and a flat characteristic in stopband. The filter specifications are as follows:
(26)
N
=
15
,
M
=
6
,
K
s
=
11
,
H
d
(
e
j
ω
)
=
e
-
j
τ
d
ω
,
0
≤
ω
≤
0.3
π
,
0
,
ω
=
π
,
τ
d
=
9,10,11,12,13
.
In the proposed method, we set the maximum allowable pole radius r
m to 0.98, and the number of grid points is L
=
500 and R
=
1000. Note that the conventional method [12] cannot guarantee the stability of the obtained filter.The performance of the obtained filters is shown in Table1 which are obtained by the conventional method [12] based on Remez algorithm. In Table 1, R
p is the maximum magnitude error in passband, G
err
⁡ is the maximum group delay error in passband, and P
max
⁡ is the maximum pole radius of the obtained filter. From Table 1, it is confirmed that the proposed method needs more computational cost (CPU time and iteration number) than the conventional method, but the characteristic of the obtained filters by the proposed method is about the same as or better than them by the conventional method. Particularly, the filter obtained by the conventional method is unstable when τ
d
=
9. The stability of the filters obtained by the method which cannot guarantee the stability depends on a given group delay specification. In fact, when the group delay τ
d is less than 9.5, all filters obtained by conventional method were unstable. In contrast, the proposed method can always obtain the stable filters because the maximum pole radius of the filter obtained can be prespecified. Therefore, the proposed method can design the filters that cannot be designed by the conventional method.Table 1
Comparison with [12] in Example 1.
τ
d
R
p
G
err
⁡
P
max
⁡
CPU time (s)
Iteration number
Proposed
9
5.622 × 10−5
3.618 × 10−3
0.980
1.334
16
10
1.248 × 10−5
7.741 × 10−4
0.952
0.587
10
11
4.416 × 10−5
6.501 × 10−3
0.890
0.537
10
12
1.098 × 10−4
3.224 × 10−2
0.858
0.510
8
13
2.060 × 10−4
1.146 × 10−1
0.847
0.525
9
Reference  [12]
9
2.281 × 10−6
2.657 × 10−5
1.051
0.121
6
10
1.269 × 10−5
7.257 × 10−4
0.953
0.134
6
11
4.608 × 10−5
6.499 × 10−3
0.891
0.124
6
12
1.168 × 10−4
3.238 × 10−2
0.858
0.124
6
13
3.129 × 10−4
1.122 × 10−1
0.846
0.130
6


## 5.2. Example 2

Next we design the inverse Chebyshev-type IIR filters with the real coefficients in order to compare with the conventional method [13]. The filter specifications are as follows:
(27)
N
=
12
,
M
=
5
,
K
p
=
10
,
H
d
(
e
j
ω
)
=
e
-
j
τ
d
ω
,
ω
=
0
,
0
,
0.5
π
≤
ω
≤
π
,
τ
d
=
7,8
,
9,10,11
.
In the proposed method, we set the maximum allowable pole radius r
m to 0.98, and the number of grid points is L
=
500 and R
=
1000. Note that the conventional method [13] cannot guarantee the stability of the obtained filter.The numerical performance of the obtained filters is shown in Table2 which are obtained by the conventional method [13] based on Remez algorithm. In Table 2, R
s is the minimum stopband attenuation in dB. From Table 2, it is confirmed that the proposed method needs more computational cost than the conventional method but gives slightly better characteristics. Moreover, when the group delay τ
d is less than 7.2, the filters obtained by conventional method were unstable. In contrast, the proposed method can always obtain the stable filters. Therefore, the proposed method can design the filters that cannot be designed by the conventional method.Table 2
Comparison with [13] in Example 2.
τ
d
R
s (dB)
P
max
⁡
CPU time (s)
Iteration number
Proposed
7
27.85
0.980
1.812
45
8
35.02
0.928
0.512
16
9
41.38
0.850
0.372
14
10
46.61
0.784
0.385
12
11
50.89
0.746
0.351
11
Reference  [13]
7
25.78
1.002
0.212
7
8
34.07
0.924
0.205
7
9
40.46
0.850
0.202
7
10
45.76
0.784
0.207
7
11
50.10
0.746
0.198
6Next, we designed the filter withτ
d
=
9 for comparison with [14]. The minimum stopband attenuation of the filter obtained was 41.38dB in the proposed method and 40.04dB in [14]. It is thought that this difference is because the performance of the filter obtained by the conventional method is dependent on the given initial value.


## 5.3. Example 3

Next, as the example of the filter with the complex coefficients, we design the inverse Chebyshev-type IIR filter with the following specifications:(28)
N
=
12
,
M
=
5
,
H
d
(
e
j
ω
)
=
e
-
j
9
ω
,
ω
=
0.1
π
,
0
,
-
π
≤
ω
≤
-
0.6
π
,
0.7
π
≤
ω
<
π
,
K
p
=
9,10,11
.
In this example, we use r
m
=
1.0, L
=
500, and R
=
1000.The proposed method took an average of 0.78 seconds for convergence. The magnitude response and group delay response of the obtained filter are shown in Figures1, 2, and 3. From these figures, it is confirmed that the magnitude and group delay responses have a flat characteristic at ω
=
0.1
π and the magnitude response in stopband is equiripple. Table 3 is the numerical performance of the obtained filter.Table 3
The performance of the filters in Example 3.
K
p
R
s (dB)
P
max
⁡
9
94.43
0.902
10
69.92
0.865
11
48.03
0.856Figure 1
Overall magnitude response of the proposed IIR filters in Example 3.Figure 2
Passband magnitude response of the proposed IIR filters in Example 3.Figure 3
Group delay response of the proposed IIR filters in Example 3.


## 5.4. Example 4

Finally, we show the design examples of inverse-Chebyshev-type FIR filters. The filter specifications are as follows:(29)
N
=
1999
,
M
=
0
,
H
d
(
e
j
ω
)
=
e
-
j
900
ω
,
ω
=
0
,
0
,
0.0080
π
≤
ω
≤
π
,
K
p
=
10,12,14
.
In this example, we set the number of the grit points L to 10,000.The proposed method took an average of 2,794 seconds for convergence. The magnitude response and group delay response of the obtained filter are shown in Figures4, 5, and 6. From these figures, it is confirmed that the magnitude and group delay responses have a flat characteristic at ω
=
0 and the magnitude response in stopband is equiripple. The minimum stopband attenuation of the obtained filters was 72.98 (dB) for K
p
=
10, 59.11 (dB) for K
p
=
12, and 47.19  (dB) for K
p
=
14. We confirmed that these filters cannot be designed by the conventional method [13]. Hence, the proposed method can also design the high-order filters that cannot be designed by the conventional method.Figure 4
Overall magnitude response of the proposed FIR filters in Example 4.Figure 5
Passband magnitude response of the proposed FIR filters in Example 4.Figure 6
Group delay response of the proposed FIR filters in Example 4.


## 6. Conclusion

In this paper, a design method based on quadratic programming has been proposed for approximately linear phase IIR filters and FIR filters with prescribed flatness in passband or stopband. To guarantee the stability of the filter obtained, we applied the extended positive realness to the QP problem. With this method, stable IIR filters and FIR filters that are Chebyshev type and inverse Chebyshev type can be easily designed. In the design examples, we showed that the performance of the filter obtained by the proposed method is about the same as or better than the conventional method, and the filters which cannot be designed by the conventional method can also be designed.



---

*Source: 304174-2015-02-25.xml*