Discrete Fractional Hilbert Transform Soo-Chang Pei and Min-Hung Yeh
Abstract—The Hilbert transform plays an important role in the theory and practice of signal processing. A generalization of the Hilbert transform, the fractional Hilbert transform, was recently proposed, and it presents physical interpretation in the definition. In this paper, we develop the dis-crete fractional Hilbert transform, and apply the proposed disdis-crete frac-tional Hilbert transform to the edge detection of digital images.
Index Terms—Fractional Hilbert transform, fractioinal Fourier trans-form, Hilbert transform.
I. INTRODUCTION
The Hilbert transform is an important tool for signal processing, and it has been widely used in many areas, such as modulation theory [1], edge detection [2], [3], and so on. Besides the continuous Hilbert trans-form, the discrete Hilbert transform can also be used for digital commu-nication and edge detection of digital images [1]–[3]. A generalization of Hilbert transform, the fractional Hilbert transform, was proposed in [4], and it provides a tool to process signal in the fractional Fourier plane instead of a conventional Fourier plane.
The method for implementing the fractional Hilbert transform in [4] is using optical instruments. The goal of this paper is to develop the discrete fractional Hilbert transform, which can have similar outputs as those of the continuous fractional Hilbert transform.
II. PRELIMINARY
A. The Continuous Hilbert Transform
The conventional Hilbert transform of a continuous signalx(t) is computed as [1]
^x(t) = 1
01
x( )
t 0 d : (1)
The continuous Hilbert transform consists of a=2 radian phase shift (for positive frequencies only) in the frequency domain [6]. Thus the transfer function of Hilbert transform becomes
H1(!) =
j; ! > 0 0; ! = 0 0j; ! < 0:
(2)
B. The Continuous Fractional Hilbert Transform
In [4], two alternative definitions for the continuous fractional Hilbert transform have been developed. One is based upon the mod-ification of spatial filter with a fractional parameter, and its transfer function is defined as
HP(v) = cos H0(v) + sin H1(v) (3)
Manuscript received August 1998; revised July 2000. This paper was recom-mended by Associate Editor V. Madisetti.
S.-C. Pei is with the Department of Electrical Engineering, National Taiwan University, Taipei, Taiwan, R.O.C. (e-mail: [email protected]).
M.-H. Yeh is with the Department of Electronic Engineering, National I-Lan Institute of Technology, I-Lan, Taiwan, R. O. C.
Publisher Item Identifier S 1057-7130(00)09938-9.
Fig. 1. Block diagrams for the different implementations of the fractional Hilbert transform. (a) Spatial filter. (b) FRFT method. (c) Generalized definition.
where = P =2. The above definition of fractional Hilbert transform is a weighted sum of the original signal and its conventional Hilbert transform, and it is based upon modifying the spatial filter with frac-tional parameter.
The other fractional Hilbert transform is based upon the fractional Fourier transform (FRFT) [5]. The FRFT operation indicates a rotation of signal in the time-frequency plane. The transform kernel of FRFT defined in [5] is
K (t; u) = 1 0 j cot 2 ej((t +u )=2) cot 0jut csc (4) where indicates the rotation angle in the time-frequency plane. While = =2, the FRFT will become conventional Fourier transform.
The transfer functionVQfor the other fractional Hilbert transform based upon the FRFT method is defined as [4]
VQ= F0QH1FQ (5)
whereFQis the fractional Fourier transform with fractional orderQ. WhileQ = 1, the FRFT becomes the conventional Fourier transform. The parameterQ defined here is equal to Q=2.
The above two definitions of the fractional Hilbert transform can be merged into a general one [4]. Thus, its transfer function is defined as follows:
HP; Q= F0QHPFQ: (6)
Fig. 1 shows the block diagrams for implementing the fractional Hilbert transform. IfP = 1, the second fractional Hilbert transform which is based upon FRFT is obtained.
C. The Discrete Hilbert Transform
The transfer function of the discrete Hilbert transform is defined as [6], [7] H(!) = j; 0 < ! < 0; ! = 0 and w = 0j; 0 < ! < 0: (7)
Many methods for computing the discrete Hilbert transform have been proposed [6], [7]. Most of them are based upon the transfer func-tion of the Hilbert transform. The method for computing the discrete 1057–7130/00$10.00 © 2000 IEEE
Step 1: Compute the DFT of signalfx[k]g
X[n] = DFT[x[k]] (8)
Step 2: X[n] is multiplied by the mask M1. The maskM1is defined as ifN is even M1= [0; j; j; . . . ; j (N=2)01 ; 0; 0j; 0j; . . . ; 0j (N=2)01 ] (9) ifN is odd M1= [0; j; j; . . . ; j (N01)=2 ; 0j; 0j; . . . ; 0j (N01)=2 ] (10)
Step 3: Compute the inverse DFT to obtain^x[k].
^x[k] = IDFT[X[n]M1[n]]: (11) Then^x[k] will be the discrete Hilbert transform of x[k]. Block diagram for implementing the discrete Hilbert transform is shown in Fig. 2.
III. DEVELOPMENT OF THEDISCRETEFRACTIONALHILBERT TRANSFORM
Similar to the continuous fractional Hilbert transform, a discrete fractional Fourier transform (DFRFT) will be required in the gener-alized discrete Hilbert transform. But in the history of DFRFT devel-opment, the DFRFT has been considered a linear combination of the signal and its spectrum in many documents [8]. Such a definition cannot have similar outputs as those of continuous fractional Fourier transform [9]. We have found that the DFRFT, with discrete Hermite eigenvectors and appropriate eigenvalue assignment rule, can have similar results as those of continuous FRFT [10]. Here we will use our DFRFT compu-tation method to develop the discrete fractional Hilbert transform. The kernel of DFRFT is defined as follows:
FQ= n
e0jnQ(=2)v
nvn3 (12)
wherevnis thenth order DFT Hermite eigenvector which is a DFT eigenvector with similar shape as thenth order Hermite function. The method for finding the nth order DFT Hermite eigenvector can be found in [10] and [11]. The method in [10] is not good enough. Two more accurate methods can be found in [11]. The discrete fractional Hilbert transform can be computed through the following steps.
Step 1: Compute the DFRFT of signalfx[k]g with parameter Q.
XQ[n] = DFRFTQ[x[k]] (13)
Step 2: XQis multiplied by the maskMP. The maskMP is defined as ifN is even MP = [cos ; ej ; ej ; . . . ; ej (N=2)01 ; cos ; e0j ; e0j ; . . . ; e0j (N=2)01 ]: (14) ifN is odd MP = [cos ; ej ; ej ; . . . ; ej (N01)=2 ; e0j ; e0j ; . . . ; e0j (N01)=2 ] (15) where = P =2.
Step 3: Compute the DFRFT with parameter0Q.
^x[k] = DFRFT0Q[XQ[n]MP[n]]: (16) Then^x[k] is the discrete fractional Hilbert transform of x[k]. The responses in the fractional Fourier domain are definedej and e0j for positive and negative transform domains, respectively. The
first and central entries in the maskMPare both equal tocos , which are defined as the middle responses for positive and negative transform domains
cos = ej + e0j
2 (17)
Then it can be easily verified that (9) and (10) are just special cases of (14) and (15). So the proposed discrete fractional Hilbert transform is a generalized version of the conventional discrete Hilbert transform. WhileP = 0, the mask becomes
M0= [1; 1; 1; . . . ; 1] (18) Thus the output of discrete fractional Hilbert transform will be the same as the input signal. In this case, the discrete fractional Hilbert transform will become an identity transform. WhileP = 2, the mask becomes
M2= [01; 01; 01; . . . ; 01]: (19) The output of discrete fractional Hilbert transform becomes the nega-tive value of the input signal.
The block diagram shown in Fig. 2 can also be modified for the implementation of the discrete fractional Hilbert transform. The DFT and IDFT must be changed into DFRFT with parameterQ and 0Q, respectively. And the maskMP in (14) and (15) must be used for the mask block in Fig. 2.
Example 1: The amplitudes of discrete fractional Hilbert transform
for a rectangular window are shown in Figs. 3 and 4. It can be observed that the discrete fractional Hilbert transform of a rectangular function consists of two peaks which mark the edges in the signal. The empha-sizes of positive or negative edges are based upon the parameters for discrete fractional Hilbert transform. While0 < P < 1, the posi-tive edges are emphasized. And the negaposi-tive edges will be emphasized when1 < P < 2. But in the case of Q = 0:5, it can be observed in Fig. 4 that there is no preference either for the negative or for the posi-tive derivaposi-tive.
The results in Example 1 are very similar to those of the continuous fractional Hilbert transform for the results in [4]. This can help us to verify that the proposed discrete fractional Hilbert transform is our de-sired transform.
IV. PROPERTIES OF THEDISCRETEFRACTIONALHILBERTTRANSFORM 1) Periodicity: The period of continuous fractional Hilbert trans-form has the periods 4 for both parameters (P and Q). In the discrete fractional Hilbert transform, this property can also be preserved.
2) Angle Addition: The continuous fractional Hilbert transform has angle addition property for parameterP . Now we will discuss the angle addition property in the discrete case. For a discrete
Fig. 3. The discrete fractional Hilbert transform of a rectangular windowQ = 1.
Fig. 4. The discrete fractional Hilbert transform of a rectangular windowQ = 0:6. signalx = [x0; x1; . . . ; xl01], the dc and ac components of
discrete signalx are defined as follows: xDC= l01 i=0 xi (20) xAC= l01 i=0 (01)ixi: (21)
The angle addition property of discrete fractional Hilbert trans-form can be preserved while the dc and ac components of signal are removed
~x = x 0 xDC0 xAC: (22)
Then the following equation will be satisfied:
HP +P ; 1(~x) = HP ; 1HP ; 1(~x) (23)
whereHP; Qindicates the discrete fractional Hilbert transform with parameterP and Q. The proof of (23) will be straightfor-ward
HP ; 1HP ; 1= (F01MP F)(F01MP F)
= F01M
P MP F:
If the dc and ac components are removed, thenMP MP =
Fig. 5. Results of Example 2: edge detection of a square by the discrete fractional Hilbert transform.Q and Q are equal to 1. P and P are different in each plot.
the outputs of DFRFT are with zero values in the first and central entries.
V. APPLICATIONS OF THEDISCRETEFRACTIONALHILBERT TRANSFORM
The conventional discrete Hilbert transform has been applied to find the edges of digital images [2], [3]. According to the results shown in Example 1, the discrete fractional Hilbert transform can emphasis pos-itive or negative edges for digital signals. In the following example, we will apply the discrete fractional Hilbert transform to detect the edges for digital images. The principle used for edge detection through the discrete fractional Hilbert transform is based upon the idea in [2]. The edges occur in the(m; n) point if the following equation is satisfied:
jhP ; Q (m; n)j2+ jhP ; Q (m; n)j2> threshold (24) wherehP ; Q (m; n) is the output of the discrete fractional Hilbert transform with parametersPxandQxin thex-direction for the point (m; n). hP ; Q (m; n) is the output of the discrete fractional Hilbert
transform with parametersPyandQyin they-direction for the point (m; n). The choices of threshold values are to control the amounts of edges in detection. The selections of parameters (Px,Px,Qx,Qy) are depended upon the desired directional edges in images.Px = 0
indicates no detection in thex-direction; 0 < Px< 1 emphasizes the
positive edges in thex-direction. Px = 1 indicates the edge detection
in thex-direction no matter of positive or negative edges. 1 < Px< 2 emphases the negative edge in thex-direction. These cases are also the same for thePyparameter in they-direction. It must be noted that (24)
use the square of amplitude as measures. It is because transform results of the discrete fractional Hilbert transform are complex numbers.
Example 2: In this example, we will apply the discrete fractional
Hilbert transform for edge detection. The original image is drawn in the upper left corner of Fig. 5, and it is a simple square. The other fifteen images are the detection results. It can be observed that the edges in digital images can be detected through the choices of parameterPxand Py. The parametersQxandQyused in this example are both equal to 1. WhilePx = 0:5 and Py = 0, only the horizontal positive edges
are emphasized. WhenPx = 1:5 and Py = 0, only the horizontal
negative edges are emphasized. These results can be viewed clearly in Fig. 5. The case,Px = 1 and Py = 1, is the conventional discrete Hilbert transform for edge detection, and all directions of edges can be viewed in this case.
From the results shown in Example 2, we know that the positive or negative edges can be obtained through the choices of fractional Hilbert transform parameters. These results cannot be obtained by common edge detectors, but the fractional Hilbert tranform method can achieve them.
VI. CONCLUSION
A method for computing the generalized discrete fractional Hilbert is developed in this brief: forward DFRFT, masking DFRFT, and in-verse DFRFT. Appropriate masks in computing the discrete fractional Hilbert transform are proposed. The proposed discrete fractional Hilbert transform can have similar results as those of continuous fractional Hilbert transform. Moreover, the properties of continuous Hilbert transform can also be preserved. And, the proposed discrete
fractional Hilbert transform can be successfully used in edge and corner detections for digital images.
REFERENCES
[1] R. E. Ziemer and W. H. Tranter, Principles of Communications—Sys-tems, Modulation, and Noise. Boston, MA: Houghton Mifflin, 1990. [2] G. M. Livadas and A. G. Constantinides, “Image edge detection and
segmentation based on the Hilbert transform,” in Proc. IEEE Int. Conf. Acoustics, Speech, and Signal Processing, 1988, pp. 1152–1155. [3] K. Kohlmann, “Corner detection in natural images based on the 2-D
Hilbert transform,” Signal Processing, vol. 48, pp. 225–234, 1996. [4] A. W. Lohmann, D. Mendlovic, and Z. Zalevsky, “Fractional Hilbert
transform,” Opt. Lett., vol. 21, pp. 281–283, Feb 1996.
[5] L. B. Almeida, “The fractional Fourier transform and time-frequency representation,” IEEE Trans. Signal Process., vol. 42, pp. 3084–3091, Nov. 1994.
[6] A. V. Oppenheim, Discrete-Time Signal Processing: Prentice-Hall Inter-national Inc., 1989.
[7] L. B. Jackson, Digital Filters and Signal Processing. Norwell, MA: Kluwer, 1989.
[8] B. Santhanam and J. H. McClellan, “The discrete rotational Fourier transform,” IEEE Trans. Signal Process., vol. 42, pp. 994–998, Apr. 1996.
[9] H. M. Ozaktas, O. Arikan, M. A. Kutay, and G. Bozdagi, “Digital computation of the fractional Fourier transform,” IEEE Trans. Signal Process., vol. 44, pp. 2141–2150, Sept. 1996.
[10] S. C. Pei and M. H. Yeh, “Improved discrete fractional Fourier trans-form,” Opt. Lett., vol. 22, pp. 1047–1049, July 15, 1997.
[11] S. C. Pei, M. H. Yeh, and C. C. Tseng, “Discrete fractional Fourier trans-form based on orthogonal projection,” IEEE Trans. Signal Processing, vol. 47, pp. 1335–1348, May 1999.
Design of Nonuniform Multirate Filter Banks by Semidefinite Programming
Aryan Saadat Mehr and Tongwen Chen
Abstract—In this brief, the design of finite-impulse response (FIR) filter banks by semidefinite programming is discussed. The initial analysis filters are designed according to the characteristics of the input. By the design procedure, for the given set of analysis filters, synthesis filters are found so that the norm of the error system is minimized over all synthesis filters that have a prespecified order. Then, the synthesis filters obtained in the previous step are fixed and the analysis filters are found similarly. By iteration, the norm of the error system decreases until it converges to its final value.
Index Terms—Filter banks, multirate systems, optimization.
I. INTRODUCTION
As mentioned in [5] and [6], it is usually possible to relate a nonuni-form filter bank to a uninonuni-form filter bank with possibly interrelated fil-ters. Thus, the design of a nonuniform filter bank can be converted
Manuscript received September 1999; revised June 2000. This work was sup-ported by the Natural Sciences and Engineering Research Council of Canada. This paper was recommended by A. Skodras.
The authors are with the Department of Electrical and Computer Engi-neering, University of Alberta, Edmonton, AB, T6G 2G7, Canada (e-mail: [email protected]).
Publisher Item Identifier S 1057-7130(00)09936-5.
to the design of a uniform filter bank subject to some structural straints. The design process should be capable of handling these con-straints. For example, methods presented in [3] and [6] are suitable for the cases where no structural constraints are present. In [4], it was shown that filter banks may be designed by model matching, i.e., by minimizing theH1norm of an error system, formed by subtracting the output of a pure delay transfer function from the output of the filter bank. By such a design method, the analysis filters are designed in ad-vance and infinite-impulse response (IIR) synthesis filters are found so that theH1norm of the error system is minimized, then the IIR filters are approximated by finite-impulse response (FIR) filters. The filters are found by solving two Riccati equations. Because of the approxima-tion, the final filters are suboptimal, furthermore this method can not accommodate the structural constraints.
In this brief, we follow an iterative approach. At each iteration, we use semidefinite programming and obtain the FIR synthesis filters for a given set of FIR analysis filters or vise versa. The problem is a convex optimization problem, and since no approximation is involved, at each iteration the solution is optimal, i.e., the FIR synthesis (analysis) filters are optimal for the given analysis (synthesis) filters. Here, we consider theH1norm as the optimality criteria. Thus the designed filter bank is closest to the desired ideal system in the worst case scenario. As we will see, the constraints will not pose any difficulty in the design process.
The semidefinite programming (SDP) problem is the optimization problem of a linear function subject to the constraint that a matrix be positive definite. In other words, the following problem is a semidefi-nite programming problem:
minimize c0x; subject to G(x) > 0 where2 Rmis the variable, and
G(x) = G0+ m i=1
xiGi
and the given matricesG0; 1 1 1 ; Gm 2 Rn2nare symmetric. Here, for real symmetric matricesA and B, A > B whenever A 0 B is positive definite. The inequalityG(x) > 0 is called a linear matrix in-equality (LMI). The SDP problems are convex optimization problems and can be solved using interior point methods. Thus SDP problems are polynomial time solvable, if an a priori bound on their solution is known [1], [2].
This brief is organized as follows. In Section II, we discuss the model-matching formulation for filter banks. In the third section, this problem is then converted to an SDP problem. In Section IV, we give an example for the design of a three channel nonuniform filter bank. The example involves periodic blocks in the synthesis filter bank and frequency selective filters as the analysis filters. Finally, in Section V, we make some concluding remarks.
II. FORMULATION
A nonuniform filter bank as shown in Fig. 1 is considered. In this section, we will discuss how a model matching problem for the design of multirate filter banks can be obtained.
A nonuniform filter bank is a periodic system with period q = lcm(q0; q1; 1 1 1 ; qm01), where qiare the downsampling factors. Therefore, if we block the input and output signals, a multi-input multi-outputq by q LTI system results. In [8], the building blocks of this filter bank are studied and in [3], the transfer matrix of a blocked 1057–7130/00$10.00 © 2000 IEEE
