- Email: [email protected]
A multi-input–multi-output system approach for the computation of discrete fractional Fourier transform
Signal Processing 80 (2000) 1501}1513
A multi-input}multi-output system approach for the computation of discrete fractional Fourier transform Der-Feng Huang, Bor-Sen Chen* Department of Electrical Engineering, National Tsing Hua University, Hsin Chu, 30043, Taiwan, ROC
Abstract In the last decade, there have been several approaches on the topic of computing discrete fractional Fourier transform (DFrFT). The purpose of this work is to provide a uni"ed approach to improve the computation of DFrFT. Based on results by Ozaktas et al., a more general multirate "lter structure is developed for the improvement of computing DFrFT. With introduction of block-input}block-output technique, a speci"c multirate "ltering system can be converted to an equivalent linear time-invariant (LTI) multi-input}multi-output (MIMO) system. Other DFrFT algorithms based on matrix}vector multiplication can also be implemented by the similar structure. Therefore, a uni"ed computation framework is developed for the multirate implementation of the DFrFT via the proposed MIMO system. Finally, for the implementation of DFrFT system, a corresponding fractional Fourier "lter bank (FrFB) is also constructed as an extension of the conventional discrete Fourier transform (DFT) "lter bank system. 2000 Elsevier Science B.V. All rights reserved. Zusammenfassung In der letzten Dekade gab es einige vertiefte UG berlegungen zum Theme der Berechnung der diskreten gebrochenen Fourier-Transformation (DfrFT). Der Zweck dieser Arbeit ist, eine einheitliche Methode zur Verbesserung der Berechnung der DfrDT zu "nden. Basierend auf Resultaten von Ozaktas et al. wird eine allgemeinere, Multiraten"lterstruktur entwickelt, um die Berechnung der DfrFT zu verbessern. Durch EinfuK hrung einer block-input}block-output Technik kann ein spezielles Multiraten"ltersystem in eine aK quivalentes lineares zeitinvariantes (LTI) Multi-Input}Multi-Output System (MIMO) uK berfuK hrt werden. Andere DfrFT Algorithmen, die auf Matrix}Vektor Multiplikationen basieren, koK nnen, durch eine aK hnliche Struktur implementiert werden. Aus diesem Grund wird ein einheitliches Berechnungssystem zur Multiratenimplementierung der DfrDT uK ber das vorgestellte MIMO System entwickelt. Abschlie{end wird zur Implementierung der DfrFT noch eine entsprechende gebrochene Fourier-Filterbank als Erweiterung des konventionellen, diskreten Fourier Transformations (DFT) Filterbanksystems aufgcbaut. 2000 Elsevier Science B.V. All rights reserved. Re2 sume2 Plusieurs approches pour le calcul de la transformeH e de Fourier fractionnelle discre`te (DFrFT) ont eH teH proposeH es durant la dernie`re deH cade. Le but de ce travail est de fournir une approche uni"eH e pour ameH liorer le calcul de celle-ci. Sur la base des reH sultats de Ozaktas et al., une structure de "ltre a` cadences multiples plus geH neH rale est proposeH e pour ce faire. Par introduction d'une technique entreH e par blocs/sortie par blocs, un syste`me de "ltrage a` cadences multiples speH ci"que
* Corresponding author. Fax: 886-03-5715971. E-mail address: [email protected] (B.-S. Chen). 0165-1684/00/$ - see front matter 2000 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 0 ) 0 0 0 5 2 - 9
1502
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
peut e( tre converti en un syste`me multi-entreH es/multi-sorties (MIMO) lineH aire invariant dans le temps (LTI). D'autres algorithmes de DFrFT baseH s sur un produit matrice-vecteur peuvent e`galement e( tre implanteH s avec la me( me structure. De ce fait un cadre uni"eH de calcul est deH veloppeH pour l'implantation multi-cadences de la DFrFT via le syste`me MIMO proposeH . En"n, pour l'implantation du syste`me de DFrFT un banc de "ltres de Fourier fractionnels (FrFB) est eH galement construit comme extension du syste`me de banc de "ltres de transformation de Fourier discre`te (DFT) conventionnel. 2000 Elsevier Science B.V. All rights reserved. Keywords: Discrete fractional Fourier transform; Multi-input}multi-output system; Fractional Fourier "lter bank
1. Introduction In the recent years, the fractional Fourier transform (FrFT) [1,5,6,9,12,13], the generalization of the Fourier transform (FT), had been investigated in many engineering applications. There have been several approaches on the topic of computing DFrFT. For example, the FrFT can be implemented by the optical system as the propagation in a gradient-index medium [7]. However, it can also be simulated with a digital algorithm, a practical way to compute the FrFT was introduced in [11]. Suppose we have a sequence of +s(n), which uniquely characterizes the function s(t), the design objective is to "nd +s?(n),, the samples sequence of s?(t) in the fractional frequency domain. If we let s and s? denote column vectors with M elements containing the samples of s(t) and s?(t), respectively, the DFrFT is an operator that maps the sample vector s of the original function into the sample vector s? in the fractional frequency domain. This mapping may be represented by the matrix F? which we call the DFrFT matrix. In [11], the DFrFT s?(n) of a signal s(n) with sampling frequency f can be computed as (see Fig. 1) s?"F?s"D KK?KJ s, where K? is the DFrFT matrix with entries A K?(m, n)" ( (fe p \@KL(D , 2
(1.1)
"m")M, "n")M,
K is a diagonal matrix with entries (A /2)(fe p ?K (D , "a p/2, A "e\ p (> ( /("sin ", ( ( a"cot , b"csc , D and J denote the 2-decimation and 2-interpolation operations, respectively. Under the same assumption, the second algorithm by Ozaktas et al. [11] is also shown in Fig. 2. This algorithm can also be implemented in a similar form as (1.1) in the following: s?"F?s"D KH KJ s, JN
Fig. 1. The schematic implementation of DFrFT proposed by Ozaktas et al. in [11].
(1.2)
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1503
Fig. 2. The second algorithm proposed by Ozaktas et al. in [11].
where D and J denote the 2-decimation and 2-interpolation operations, respectively, K is the diagonal matrix with entries (A /2)(fe p \@K (D , and H is a convolution operation with transfer function ( JN + h(z)" h(n)z\L, (1.3) L\+ with h(n)"A (
D
1
\D (b
e p e\ pT @e pTLD dv.
(1.4)
More details are referred in [11]. However, we shall elaborate on this form in z-domain and show that the two algorithms can be implemented by a similar system. There are also several algorithms which do not have to perform the multirate operations, for example, the DFrFT proposed by Santhanam et al. [15] is constructed as s?"F?s" p (a)=Gs, 0)a)2p, G G where W is the ordinary DFT matrix and
(1.5)
p (a)"[(1#e ?) cos a],
p (a)"[(1!je ?) sin a], 1 p (a)"[(!1#e ?) cos a], p (a)" [(!1!je ?) sin a]. 2 In particular, this type of DFrFT can be computed by the weighted summation of the signal, the reverse order of the signal, and its discrete Fourier transform (DFT) and the reverse of the DFT. That is, s?(n)"F?s(n)"p (a)s(n)#p (a)s(!n)#p (a)S(n)#p (a)S(!n), where S(n) is the DFT of s(n) and s(!n), S(!n) are the #ip versions of s(n), S(n), respectively. The computation can be also implemented by fast Fourier transform. In [14], Pei et al. proposed a similar matrix}vector multiplication algorithm. The fractional power of matrix can be calculated from the eigendecomposition of Hermite Gauss polynomials as follows: s?"F?s"UK D? pUK 2s, 0)a)2p,
(1.6)
where D? p is a diagonal matrix that corresponding to the fractional power of the FrFT and UK is the matrix containing of the normalized eigenvector corresponding to the mth-order Hermite function. And if the
1504
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
number of sample M is even, UK "[u( "u( "u( "2"u( "u( ], +\ + D? p"diag(e\ , e\ ?, e\ ?,2, e\ +\?, e\ +?) and if the number of sample M is odd, UK "[u( "u( "u( "2"u( "u( ], +\ +\ D? p"diag(e\ , e\ ?, e\ ?,2, e\ +\?, e\ +\?). More details can be found in [14]. The above methods involve multirate operations, matrix}vector multiplication and/or convolution. However, since the computation is performed with one input/one output, it takes much more time to calculate DFrFT. Since there is still no e$cient computation algorithm from practical application point of view, it is more appealing to "nd a more fast algorithm to unify various FrFT computation algorithms for practical implementation and further applications. In this paper, we propose an equivalent multi-input}multi-output (MIMO) linear time-invariant (LTI) system through blocking (s/p) and unblocking (p/s) for these algorithms such that these computations can be improved through an MIMO system with polyphase structure. An MIMO system with M-decimation can operate at a rate which is M times faster than the rate of one-input/output system. The price paid for this is that the implementation of the blocking "lter H(z) of DFrFT requires a larger number of computational units. As a result, we are able to process signals which arrive at M times larger rate. With this structure, an equivalent "lter bank can be designed to construct the discrete fractional Fourier transform with faster speed. In this paper, we propose a uni"ed implementation to improve the computation of the FrFT. Two types of computational algorithms for FrFT are converted to a similar structure of the MIMO LTI system by block-input}block-output mechanism. In Section 2, we introduce a special multirate "lter system and an equivalent block-input}block-output structure for DFrFT. In Section 3, we show that all these algorithms introduced can be implemented by the proposed polyphase system. In Section 4, according to these DFrFT polyphase systems, we construct a corresponding fractional Fourier "lter bank (FrFB) systems as an extension of conventional DFT "lter bank. The last section concludes these methods. In the following sections, the bold-faced letters such as F, s denote matrix and vector, respectively.
2. A multirate 5ltering system for implementation of DFrFT In this section, we propose a more general multirate "ltering system for the DFrFT system (see Fig. 3). Then we convert the multirate DFrFT system via the blocking "ltering technique [2,3,8] to an equivalent MIMO system. Hence, all the algorithms to compute the DFrFT introduced before can be uni"ed in the framework of the proposed MIMO system. The proposed MIMO system can improve the computation rate of DFrFT. In this MIMO system, M-input/M-output data sequences are processed in parallel. This MIMO system consists of three sections: the N-fold interpolation section, the DFrFT kernel section that transforms the data vector to a proposed transformed vector, and the N-fold decimation section to decimate the sample and convert the data vector to a data sequence. Note that the vector-input and vector-output paths are indicated by heavy arrows in the following discussion. The N-fold interpolation section is an operation which interpolates (N!1) values of a signal between each samples. The N-fold interpolation section should consist of an expander (N!) which inserts (N!1) zeros between each pair of samples and an anti-imaging "lter h (z) which is a lowpass "lter with a gain of N and ' a cuto! frequency p/N. We have to point out that in a system for obtaining the interpolated values of a signal
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1505
Fig. 3. The proposed block multirate "ltering system of DFrFT.
s(n) using only discrete-time processing, the frequency response of the output of the expander is a frequencyscaled version of the original input. In this situation, an ideal lowpass "lter is referred to as a bandlimited interpolation, this can be found in most digital signal processing textbooks, for example, [10]. The interpolated signal is then passed through a serial/parallel (s/p) mechanism and the DFrFT transform which is a matrix}vector operation to obtain the signal vector s?(n) The parallel/serial (p/s) mechanism converts the ' transformed signal vector to a serial data sequence. In the N-fold decimation section, the data sequence is then passed through an anti-aliasing "lter h (z). Finally, an N-fold compressor (N ) converts the input " sequence s?(n) to the output sequence s?(n)"s?(Nn). ' ' However, the direct implementation does not illustrate the structure of DFrFT from the system point of view. We will convert the above multirate signal processing in terms of the equivalent MIMO system via the block "ltering technique [2,3,16]. First, we consider the N-fold interpolation section which consists of an anti-imaging "lter h (z) and the serial/parallel (s/p) component which converts the input signal s(n) to the ' upsampled output signal vector of length MN. We can express h (z) with scalar input}output as an equivalent ' DFrFT block "lter H (z) with block-input}block-output (see Fig. 4). More precisely, for any integers M, N ' [16,17] +,\ h (z)" z\JE (z+,), ' J J where E (z)" h (nMN#l)z\L, 0)l)MN!1 is the lth polyphase component of h' (z). Hence, the J L\ ' can be expressed as an equivalent block "lter H (z) with block input}output, that is scalar input}output system ' E (z) E (z) E (z) 2 E (z) +,\ z\E (z) E (z) E (z) 2 E (z) +,\ +,\ (z) z\E (z) E (z) \ E (z) . (2.1) H (z)" z\E +,\ +,\ +,\ ' $ 2 2 \ $ z\E (z)
z\E (z)
z\E (z)
2
E (z)
In fact, the design of h (z) is a standard IIR or FIR lowpass "lter design problem. Any of the techniques ' introduced in the traditional DSP textbooks such as [4,10] can be applied for the design of the lowpass "lter. For the anti-imaging "lter h (z) to be FIR with transfer function h (z)" * h(n)z\L, where ¸"KMN#k, ' ' L 0)k)MN!1, then )+,>I h (z)" h(n)z\L ' L "h(0)#h(1)z\#h(2)z\#h(3)z\#2#h(MN!1)z\+,\ # h(MN)z\+,#h(MN#1)z\+,>#2#h(2MN!1)z\+,\# 222
1506
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
Fig. 4. The equivalent structure of the interpolation section.
#h((K!1)MN)z\)\+,#h((K!1)MN#1)z\)\+,># 2#h(KMN!1)z\)+,\ # h(KMN)z\)+,#h(KMN#1)z\)+,>#2#h(KMN#k)z\)+,>I "E (z +,)#z\E (z +,)#2#z\+,\E (z+,), +,\ where
)\ h(nMN#l)z\L#h(KMN#l)z\J 0)l)k, E (z)" L J )\ h(nMN#l)z\L otherwise. L Let e+"[0,0,21,2,0] be the transpose of ith basis vector in the M-dimensional vector space and G "+ de"ne the operator e+
0 ,\"+ e+ J " 0 , + , ,\"+ +,"+ $ e+ +
0 ,\"+ where 0 is the (N!1);M-dimensional zero matrix. J is the operator which performs the ,\"+ + , expansion N! operation on a M-dimensional vector to produce a MN-dimensional vector. In particular, we get J "I , the M;M identity matrix for any M*1. + +"+ In fact, as the signal passes through the N-fold interpolation section, the blocking mechanism is equivalent to blocking the signal to a vector of length M at "rst, then expanding by (N!1) zeros and "nally passing through the block anti-imaging "lter H (z). We conclude the interpolation section with the following block ' mechanism (see Fig. 4): E (z) z\E (z) +,\ JM (z)"H (z)J " z\E (z) + , ' + , +,\ $ z\E (z)
E (z) , E (z) ,\ E (z) ,\ 2 z\E
,>
E (z) , E (z) ,\ E (z) ,\ 2
(z) z\E
,>
2 2 \ \
E (z) +\, E (z) +\,\ E (z) +\,\ $
(z) 2 z\E (z) +\,>
(2.2)
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1507
which is an MN;M transfer matrix of a discrete-time LTI MIMO system with M inputs and MN outputs. Remark. In particular, for FIR h (z) with length ¸"MN, ' h(0) h(1) h(2) 2 h(MN!1) z\h(MN!1)
h(0)
H (z)" z\h(MN!2) z\h(MN!1) ' $ $ z\h(1)
z\h(2)
h(1)
\
h(MN!2)
h(0)
\
h(MN!3)
$
\
.
$
z\h(3) 2
h(0)
In this situation, (2.2) becomes h(0)
h(N)
h(2N)
2
h((M!1)N)
z\h(MN!1)
h(N!1)
h(2N!1)
\
h((M!1)N!1)
z\h(MN!2) JM " + , z\h(MN!3)
h(N!2)
h(2N!2)
\
h((M!1)N!2)
h(N!3)
h(2N!3)
\
$
$
$
$
\
$
z\h(1)
z\h(N#1)
z\h(2N!1) 2 z\h((M!1)N#1)
which is a MN;M transfer matrix of a discrete-time MIMO system with M inputs and MN outputs. The compression section can be similarly implemented. In the equivalent structure of interpolation section shown in Fig. 4, the data sequence is "rst blocked to a vector of length M, s(nM) s(n)"
s(nM#1) $ s(nM#M!1)
and then passes through the operator J .Then the data sequence passes through the equivalent anti+ , imaging block "lter H (z) of h (z). Any of the techniques introduced in the traditional DSP text books such as ' ' [4,10] can be applied for the design of the lowpass "lter. The z-transform of the output vector at the interpolation section is s (z)"H (z)J s(z) ' ' + , "JM (z)s(z), + , where s (z) and s(z) is the z-transform of the vector sequences s (n) and s(n), respectively. ' ' The kernel section in Fig. 3 consists of three parts: the chirp multiplication, impulse response h(n) and another chirp multiplication. From the perspective of the MN-input vector s (n), Ks (n) is an output of an ' ' LTI MIMO system with transfer matrix K. The convolution transfer function h(z) can also be converted to the equivalent block "ltering form H(z) which is an MIMO system with MN inputs and MN outputs. Hence, the z-transform of the output vector of the kernel section is (see Fig. 5) s?(z)"KH(z)Ks (z)"K ?(z)s (z), ' ' ' where K?(z)"KH(z)K is the transfer matrix of the kernel section.
(2.3)
1508
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
Fig. 5. The equivalent MIMO system of DFrFT in Fig. 3.
Fig. 6. The general scheme of the DFrFT blocking system.
Similarly, the decimation section can be blocked as an equivalent MN;M MIMO transfer matrix, which can be expressed as DM (z)"D H (z), (2.4) + , + , " where D "J2 is a matrix which performs the N-fold decimation operation on the MN-dimensional + , + , vector s?(z) and H (z) is an equivalent block version of the anti-aliasing "lter h (z). Note that the superscript ' " " T denotes the matrix transpose. Combining (2.2)}(2.4), we obtain an equivalent block "ltering system with the following transfer matrix: H ?(z)"DM
(z)K ?(z)JM (z) (2.5) + , + , to compute the DFrFT of a given signal sequence s(n) as shown in Fig. 5. This means that the scale input/output DFrFT "ltering system can be implemented as a parallel MIMO DFrFT "ltering system (2.5). Note that H ?(z) is a transfer matrix with M-input/M-output. The structure of Fig. 6 shows the speed advantage of the blocking. In general, the sampling rate of the input signal s(n) can be M times faster than the speed of the original scale input system. However, since H ?(z) is an M;M system, it requires a larger number of computational units (multipliers and adders) than the original scale input/output system. That is, the computational parallelism is increased in the blocked implementation by increasing the number of computational units. As a result the signals are processed with M times higher rate than the rate that can normally be handled by one computational unit.
3. The DFrFT polyphase system In this section, the proposed general multirate "ltering system (2.5) is applied to the implementation of the known DFrFT algorithms (see Fig. 6). The algorithms involving multirate operations are introduced "rst,
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1509
and the algorithms without multirate operations are then shown to be the special case of the proposed system (2.5). The details are given in the following subsections in order. 3.1. The multirate DFrFT polyphase xltering system With the similar analysis in the previous section, performing the 2! operation is equivalent to performing the operator 1
0
0
$
0
0
0
0
0
$
0
0
0 1 0 J " + $ 2 2
\
\ 0
\
\ $
0
0
0
\
0
1
0
0
0
2
0
0 +"+
and 2 operation is equivalent to performing operator D "J2 . By letting K ?(z)"KK ?K and after some + + manipulations, algorithm (1.1) now can be implemented by an equivalent MIMO signal processing system with the following transfer matrix: DM (z)K? (z)JM (z) a 3[!1.5,!0.5]6[0.5,1.5], + + DM (z)K?\(z)JM (z)= otherwise, + +
H?(z)"
"
D H (z)K? (z)H (z)J a 3[!1.5,!0.5]6[0.5,1.5], + " ' + D H (z)K?\(z)H (z)J = otherwise, + " ' +
(3.1)
where a "(a mod 4)!2, H (z) and H (z) are the block versions of the anti-aliasing "lter h (z) and " ' " anti-imaging "lter h (z), respectively. Similarly, algorithm (1.2) can be implemented by an MIMO system with ' transfer matrix H?(z)"DM (z)KH(z)KJM (z), !1)a)1, + +
(3.2)
where H(z) is the block "lter of the transfer function h(z) (see (1.3)). 3.2. The DFrFT polyphase xltering system without multirate operations If the signal is sampled at twice the Nyquist rate [11], there is no need to perform the multirate operations in algorithms (1.1) and (1.2). In this situation, the DFrFT algorithms (1.1), (1.2), together with (1.5) and (1.6) all involve a matrix}vector multiplication. This means that these algorithms can be implemented by a block "ltering system, hence referred to as an MIMO system with order 0 transfer matrix K ?(z). Therefore, N is set to be 1 in the design of the DFrFT transfer matrix (2.5). Further, the anti-imaging "lter h (z),1 and the ' anti-aliasing "lter h (z),1, hence H (z)"H (z),I, therefore JM (z),I,DM (z). The system in Fig. 3 is " " ' + + now simpli"ed to a system consisting of block mechanism, kernel section and unblocking mechanism shown in Fig. 6. The system with transfer matrix (3.2) is now simpli"ed as H ?(z)"KH(z)K, !1)a)1
(3.3)
1510
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
and the transfer matrix for algorithm (3.1) is K ? (z)
K ?(z)"
a 3[!1.5,!0.5][0.5,1.5],
(3.4)
K?\(z)W otherwise.
The transfer matrix for the rotational Fourier transform (1.5) is K ?(z)" p (a)W G, G G
0)a)2p
(3.5)
and the DFrFT algorithm (1.6) can be implemented by the transfer matrix K ?(z)"UK D? pUK 2, 0)a)2p.
(3.6)
Remark. Note that (3.3) and (3.4) do not obey the additivity of rotations [14], i.e. H?>@(z)"H ?(z)H @(z).
(3.7)
4. Application to fractional Fourier 5lter banks One of the most important classes of "lter banks is the DFT "lter bank [4,16]. In the previous section, we have shown that the above DFrFT algorithms can be implemented by the MIMO transfer matrix H ?(z). In this section, we apply the polyphase systems to a generalization of the DFT "lter bank with parameter a. It is worthwhile pointing out that the block mechanism can be implemented by a delay chain followed by an M-compressor as shown in Fig. 7. Hence, the DFrFT system can be converted to an equivalent M-channel maximally decimated "lter bank system by "rst rearranging the M-compressor and the transfer matrix by noble identity (see Fig. 8) [16]. This DFrFT system can be further converted to the following equivalent "lter bank:
h? (z) h? (z) $
h? (z) +\
z\+\ z\+\
"H?(z+)
$
"
1
H? (z+) H? (z+) $
H? (z+) H? (z+) $
2 2 \
H? (z+) H? (z+) 2 H? (z+) +\ +\ +\ +\
where h? (z) is the transfer function of the kth analysis "lter de"ned as I +\ h? (z)" z\+\\JH? (z+), I I J J
H? (z+) +\ H? (z+) +\ $
k"0,2, M!1
z\+\ z\+\ $
,
1
(4.1)
in which a denotes the fractional parameter of the fractional Fourier "lter bank (aFrFB). Based on di!erent DFrFT polyphase "ltering systems, we can construct the corresponding aFrFB systems. Because the computation of the DFrFT is valid for all the above algorithms, the aFrFB construction is valid as an extension of the DFT "lter bank system.
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1511
Fig. 7. The fractional Fourier "lter bank (aFrFB) system.
Fig. 8. The noble identity for multirate operations.
In particular, when kernel (3.5) is chosen with a"p/2, the aFrFB is the same as the conventional DFT "lter bank which is known as the spectrum analyzer. It provides a unique representation for signal in terms of complex exponentials. Further, when kernel (3.6) is chosen with a"p/2, an equivalent DFT "lter bank is obtained by aFrFB, but expressing the signal in terms of Hermite functions. Apparently, DFrFT systems (1.5) and (3.6) satisfy the unitary property H ?(z)H ?(z)H"H ?(z)H\?(z)"cz\JI,
(4.2)
where superscript * denotes the complex conjugate, c is a constant and l is the system delay (possibly equal to 0). We can easily design a perfect reconstruction aFrFB system with analysis "lter bank (4.1) and synthesis "lter bank (see Fig. 9) with the following transfer functions +\ g? (z)" z\JH\?(z+), k"0,2, M!1. I I J J 5. Conclusion In this paper, we propose a special multirate "ltering system with blocking implementation for fast computation of DFrFT. With the blocking processing we have shown that the DFrFT can be implemented
1512
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
Fig. 9. The scheme of an aFrFB with perfect reconstruction property.
by an MIMO system. As is well known, the structure provides a rather e$cient way to improve the computation of DFrFT. The cost paid for the e$ciency of computation is that the components of the structure are more complex than that of the direct implementation. This blocking structure can also be implemented by an equivalent fractional "lter bank structure. Accordingly, a fractional Fourier "lter bank can be constructed by a corresponding DFrFT polyphase system. The proposed approach can be constructed by known digital signal processing techniques. Therefore, it is more suitable for practical applications.
References [1] L.B. Almeida, The fractional Fourier transform and time-frequency representations, IEEE Trans. Signal Process. 42 (1994) 3084}3091. [2] C.S. Burrus, Block implementation of digital "lters, IEEE Trans. Circuits Theory CT-18 (1971) 697}701. [3] C.S. Burrus, Block realization of digital "lters, IEEE Trans. Audio Electroacoust. AU-20 (1972) 230}235. [4] R.E. Crochiere, L.R. Rabiner, Multirate Digital Signal Processing, Prentice-Hall, Englewood Cli!s, NJ, 1983. [5] A.W. Lohmann, D. Mendlovic, Z. Zalevsky, Fractional transformations in optics, in: Process in Optics, Vol. XXXVIII, Elsevier, Amsterdam, 1998 (Chapter 4). [6] A.C. McBride, F.H. Kerr, On Namias's fractional Fourier transforms, IMA J. Appl. Math. 39 (1987) 159}175. [7] D. Mendlovic, H.M. Ozaktas, Fractional Fourier transforms and their optical implementation: I, J. Opt. Soc. Am. A 10 (9) (1993) 1875}1881. [8] S.K. Mitra, R. Gnanasekaran, Block implementation of recursive digital "lters } new strucrures and properties, IEEE Trans. Circuits Systems CAS-25 (1978) 200}207. [9] V. Namias, The fractional order Fourier transform and its application to quantum mechanics, J. Inst. Math. Appl. 25 (1980) 241}265. [10] A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, Englewood Cli!s, NJ, 1989, pp. 101}111 (Chapter 3). [11] H.M. Ozaktas, O. Arikan, M.A. Kutay, Gozde Bozdagi, Digital computation of the fractional Fourier transform, IEEE Trans. Signal Process. 44 (9) (1996) 2141}2150. [12] H.M. Ozaktas, B. Barshan, D. Mendlovic, L. Onural, Comvolution, "ltering, and multiplexing in fractional Fourier doamins and their relation to chirp and wavelets, J. Opt. Soc. Am. A 11 (1994) 547}559.
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1513
[13] H.M. Ozaktas, M.A. Kutay, D. Mendlovic, Introduction to fractional Fourier transform and applications, in: Advances in Imaging and Electron Physics, Vol. 106, Academic Press, New York, 1999, pp. 263}342 (Chapter 4). [14] S.C. Pei, M.H. Yeh, C.C. Tseng, Discrete fractional Fourier transform based on orthogonal projections, IEEE Trans. Signal Process. 47 (5) (1999) 1335}1348. [15] B. Santhanam, J.M. McCellan, The discrete rotational Fourier transform, IEEE Trans. Signal Process. 44 (4) (1996) 994}998. [16] P.P. Vadyianthan, Multirate Systems and Filter Banks, Prentice-Hall, Englewood Cli!s, NJ, 1993. [17] P.P. Vaidyanathan, S.K. Mitra, Polyphase networks, block digital "ltering, LPTV systems, and alias-free QMF banks: a uni"ed approach based on pseudocirculants, IEEE Trans. Acoust. Speech Signal Process. 36 (3) (March, 1988) 381}391.
A multi-input}multi-output system approach for the computation of discrete fractional Fourier transform Der-Feng Huang, Bor-Sen Chen* Department of Electrical Engineering, National Tsing Hua University, Hsin Chu, 30043, Taiwan, ROC
Abstract In the last decade, there have been several approaches on the topic of computing discrete fractional Fourier transform (DFrFT). The purpose of this work is to provide a uni"ed approach to improve the computation of DFrFT. Based on results by Ozaktas et al., a more general multirate "lter structure is developed for the improvement of computing DFrFT. With introduction of block-input}block-output technique, a speci"c multirate "ltering system can be converted to an equivalent linear time-invariant (LTI) multi-input}multi-output (MIMO) system. Other DFrFT algorithms based on matrix}vector multiplication can also be implemented by the similar structure. Therefore, a uni"ed computation framework is developed for the multirate implementation of the DFrFT via the proposed MIMO system. Finally, for the implementation of DFrFT system, a corresponding fractional Fourier "lter bank (FrFB) is also constructed as an extension of the conventional discrete Fourier transform (DFT) "lter bank system. 2000 Elsevier Science B.V. All rights reserved. Zusammenfassung In der letzten Dekade gab es einige vertiefte UG berlegungen zum Theme der Berechnung der diskreten gebrochenen Fourier-Transformation (DfrFT). Der Zweck dieser Arbeit ist, eine einheitliche Methode zur Verbesserung der Berechnung der DfrDT zu "nden. Basierend auf Resultaten von Ozaktas et al. wird eine allgemeinere, Multiraten"lterstruktur entwickelt, um die Berechnung der DfrFT zu verbessern. Durch EinfuK hrung einer block-input}block-output Technik kann ein spezielles Multiraten"ltersystem in eine aK quivalentes lineares zeitinvariantes (LTI) Multi-Input}Multi-Output System (MIMO) uK berfuK hrt werden. Andere DfrFT Algorithmen, die auf Matrix}Vektor Multiplikationen basieren, koK nnen, durch eine aK hnliche Struktur implementiert werden. Aus diesem Grund wird ein einheitliches Berechnungssystem zur Multiratenimplementierung der DfrDT uK ber das vorgestellte MIMO System entwickelt. Abschlie{end wird zur Implementierung der DfrFT noch eine entsprechende gebrochene Fourier-Filterbank als Erweiterung des konventionellen, diskreten Fourier Transformations (DFT) Filterbanksystems aufgcbaut. 2000 Elsevier Science B.V. All rights reserved. Re2 sume2 Plusieurs approches pour le calcul de la transformeH e de Fourier fractionnelle discre`te (DFrFT) ont eH teH proposeH es durant la dernie`re deH cade. Le but de ce travail est de fournir une approche uni"eH e pour ameH liorer le calcul de celle-ci. Sur la base des reH sultats de Ozaktas et al., une structure de "ltre a` cadences multiples plus geH neH rale est proposeH e pour ce faire. Par introduction d'une technique entreH e par blocs/sortie par blocs, un syste`me de "ltrage a` cadences multiples speH ci"que
* Corresponding author. Fax: 886-03-5715971. E-mail address: [email protected] (B.-S. Chen). 0165-1684/00/$ - see front matter 2000 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 0 ) 0 0 0 5 2 - 9
1502
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
peut e( tre converti en un syste`me multi-entreH es/multi-sorties (MIMO) lineH aire invariant dans le temps (LTI). D'autres algorithmes de DFrFT baseH s sur un produit matrice-vecteur peuvent e`galement e( tre implanteH s avec la me( me structure. De ce fait un cadre uni"eH de calcul est deH veloppeH pour l'implantation multi-cadences de la DFrFT via le syste`me MIMO proposeH . En"n, pour l'implantation du syste`me de DFrFT un banc de "ltres de Fourier fractionnels (FrFB) est eH galement construit comme extension du syste`me de banc de "ltres de transformation de Fourier discre`te (DFT) conventionnel. 2000 Elsevier Science B.V. All rights reserved. Keywords: Discrete fractional Fourier transform; Multi-input}multi-output system; Fractional Fourier "lter bank
1. Introduction In the recent years, the fractional Fourier transform (FrFT) [1,5,6,9,12,13], the generalization of the Fourier transform (FT), had been investigated in many engineering applications. There have been several approaches on the topic of computing DFrFT. For example, the FrFT can be implemented by the optical system as the propagation in a gradient-index medium [7]. However, it can also be simulated with a digital algorithm, a practical way to compute the FrFT was introduced in [11]. Suppose we have a sequence of +s(n), which uniquely characterizes the function s(t), the design objective is to "nd +s?(n),, the samples sequence of s?(t) in the fractional frequency domain. If we let s and s? denote column vectors with M elements containing the samples of s(t) and s?(t), respectively, the DFrFT is an operator that maps the sample vector s of the original function into the sample vector s? in the fractional frequency domain. This mapping may be represented by the matrix F? which we call the DFrFT matrix. In [11], the DFrFT s?(n) of a signal s(n) with sampling frequency f can be computed as (see Fig. 1) s?"F?s"D KK?KJ s, where K? is the DFrFT matrix with entries A K?(m, n)" ( (fe p \@KL(D , 2
(1.1)
"m")M, "n")M,
K is a diagonal matrix with entries (A /2)(fe p ?K (D , "a p/2, A "e\ p (> ( /("sin ", ( ( a"cot , b"csc , D and J denote the 2-decimation and 2-interpolation operations, respectively. Under the same assumption, the second algorithm by Ozaktas et al. [11] is also shown in Fig. 2. This algorithm can also be implemented in a similar form as (1.1) in the following: s?"F?s"D KH KJ s, JN
Fig. 1. The schematic implementation of DFrFT proposed by Ozaktas et al. in [11].
(1.2)
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1503
Fig. 2. The second algorithm proposed by Ozaktas et al. in [11].
where D and J denote the 2-decimation and 2-interpolation operations, respectively, K is the diagonal matrix with entries (A /2)(fe p \@K (D , and H is a convolution operation with transfer function ( JN + h(z)" h(n)z\L, (1.3) L\+ with h(n)"A (
D
1
\D (b
e p e\ pT @e pTLD dv.
(1.4)
More details are referred in [11]. However, we shall elaborate on this form in z-domain and show that the two algorithms can be implemented by a similar system. There are also several algorithms which do not have to perform the multirate operations, for example, the DFrFT proposed by Santhanam et al. [15] is constructed as s?"F?s" p (a)=Gs, 0)a)2p, G G where W is the ordinary DFT matrix and
(1.5)
p (a)"[(1#e ?) cos a],
p (a)"[(1!je ?) sin a], 1 p (a)"[(!1#e ?) cos a], p (a)" [(!1!je ?) sin a]. 2 In particular, this type of DFrFT can be computed by the weighted summation of the signal, the reverse order of the signal, and its discrete Fourier transform (DFT) and the reverse of the DFT. That is, s?(n)"F?s(n)"p (a)s(n)#p (a)s(!n)#p (a)S(n)#p (a)S(!n), where S(n) is the DFT of s(n) and s(!n), S(!n) are the #ip versions of s(n), S(n), respectively. The computation can be also implemented by fast Fourier transform. In [14], Pei et al. proposed a similar matrix}vector multiplication algorithm. The fractional power of matrix can be calculated from the eigendecomposition of Hermite Gauss polynomials as follows: s?"F?s"UK D? pUK 2s, 0)a)2p,
(1.6)
where D? p is a diagonal matrix that corresponding to the fractional power of the FrFT and UK is the matrix containing of the normalized eigenvector corresponding to the mth-order Hermite function. And if the
1504
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
number of sample M is even, UK "[u( "u( "u( "2"u( "u( ], +\ + D? p"diag(e\ , e\ ?, e\ ?,2, e\ +\?, e\ +?) and if the number of sample M is odd, UK "[u( "u( "u( "2"u( "u( ], +\ +\ D? p"diag(e\ , e\ ?, e\ ?,2, e\ +\?, e\ +\?). More details can be found in [14]. The above methods involve multirate operations, matrix}vector multiplication and/or convolution. However, since the computation is performed with one input/one output, it takes much more time to calculate DFrFT. Since there is still no e$cient computation algorithm from practical application point of view, it is more appealing to "nd a more fast algorithm to unify various FrFT computation algorithms for practical implementation and further applications. In this paper, we propose an equivalent multi-input}multi-output (MIMO) linear time-invariant (LTI) system through blocking (s/p) and unblocking (p/s) for these algorithms such that these computations can be improved through an MIMO system with polyphase structure. An MIMO system with M-decimation can operate at a rate which is M times faster than the rate of one-input/output system. The price paid for this is that the implementation of the blocking "lter H(z) of DFrFT requires a larger number of computational units. As a result, we are able to process signals which arrive at M times larger rate. With this structure, an equivalent "lter bank can be designed to construct the discrete fractional Fourier transform with faster speed. In this paper, we propose a uni"ed implementation to improve the computation of the FrFT. Two types of computational algorithms for FrFT are converted to a similar structure of the MIMO LTI system by block-input}block-output mechanism. In Section 2, we introduce a special multirate "lter system and an equivalent block-input}block-output structure for DFrFT. In Section 3, we show that all these algorithms introduced can be implemented by the proposed polyphase system. In Section 4, according to these DFrFT polyphase systems, we construct a corresponding fractional Fourier "lter bank (FrFB) systems as an extension of conventional DFT "lter bank. The last section concludes these methods. In the following sections, the bold-faced letters such as F, s denote matrix and vector, respectively.
2. A multirate 5ltering system for implementation of DFrFT In this section, we propose a more general multirate "ltering system for the DFrFT system (see Fig. 3). Then we convert the multirate DFrFT system via the blocking "ltering technique [2,3,8] to an equivalent MIMO system. Hence, all the algorithms to compute the DFrFT introduced before can be uni"ed in the framework of the proposed MIMO system. The proposed MIMO system can improve the computation rate of DFrFT. In this MIMO system, M-input/M-output data sequences are processed in parallel. This MIMO system consists of three sections: the N-fold interpolation section, the DFrFT kernel section that transforms the data vector to a proposed transformed vector, and the N-fold decimation section to decimate the sample and convert the data vector to a data sequence. Note that the vector-input and vector-output paths are indicated by heavy arrows in the following discussion. The N-fold interpolation section is an operation which interpolates (N!1) values of a signal between each samples. The N-fold interpolation section should consist of an expander (N!) which inserts (N!1) zeros between each pair of samples and an anti-imaging "lter h (z) which is a lowpass "lter with a gain of N and ' a cuto! frequency p/N. We have to point out that in a system for obtaining the interpolated values of a signal
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1505
Fig. 3. The proposed block multirate "ltering system of DFrFT.
s(n) using only discrete-time processing, the frequency response of the output of the expander is a frequencyscaled version of the original input. In this situation, an ideal lowpass "lter is referred to as a bandlimited interpolation, this can be found in most digital signal processing textbooks, for example, [10]. The interpolated signal is then passed through a serial/parallel (s/p) mechanism and the DFrFT transform which is a matrix}vector operation to obtain the signal vector s?(n) The parallel/serial (p/s) mechanism converts the ' transformed signal vector to a serial data sequence. In the N-fold decimation section, the data sequence is then passed through an anti-aliasing "lter h (z). Finally, an N-fold compressor (N ) converts the input " sequence s?(n) to the output sequence s?(n)"s?(Nn). ' ' However, the direct implementation does not illustrate the structure of DFrFT from the system point of view. We will convert the above multirate signal processing in terms of the equivalent MIMO system via the block "ltering technique [2,3,16]. First, we consider the N-fold interpolation section which consists of an anti-imaging "lter h (z) and the serial/parallel (s/p) component which converts the input signal s(n) to the ' upsampled output signal vector of length MN. We can express h (z) with scalar input}output as an equivalent ' DFrFT block "lter H (z) with block-input}block-output (see Fig. 4). More precisely, for any integers M, N ' [16,17] +,\ h (z)" z\JE (z+,), ' J J where E (z)" h (nMN#l)z\L, 0)l)MN!1 is the lth polyphase component of h' (z). Hence, the J L\ ' can be expressed as an equivalent block "lter H (z) with block input}output, that is scalar input}output system ' E (z) E (z) E (z) 2 E (z) +,\ z\E (z) E (z) E (z) 2 E (z) +,\ +,\ (z) z\E (z) E (z) \ E (z) . (2.1) H (z)" z\E +,\ +,\ +,\ ' $ 2 2 \ $ z\E (z)
z\E (z)
z\E (z)
2
E (z)
In fact, the design of h (z) is a standard IIR or FIR lowpass "lter design problem. Any of the techniques ' introduced in the traditional DSP textbooks such as [4,10] can be applied for the design of the lowpass "lter. For the anti-imaging "lter h (z) to be FIR with transfer function h (z)" * h(n)z\L, where ¸"KMN#k, ' ' L 0)k)MN!1, then )+,>I h (z)" h(n)z\L ' L "h(0)#h(1)z\#h(2)z\#h(3)z\#2#h(MN!1)z\+,\ # h(MN)z\+,#h(MN#1)z\+,>#2#h(2MN!1)z\+,\# 222
1506
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
Fig. 4. The equivalent structure of the interpolation section.
#h((K!1)MN)z\)\+,#h((K!1)MN#1)z\)\+,># 2#h(KMN!1)z\)+,\ # h(KMN)z\)+,#h(KMN#1)z\)+,>#2#h(KMN#k)z\)+,>I "E (z +,)#z\E (z +,)#2#z\+,\E (z+,), +,\ where
)\ h(nMN#l)z\L#h(KMN#l)z\J 0)l)k, E (z)" L J )\ h(nMN#l)z\L otherwise. L Let e+"[0,0,21,2,0] be the transpose of ith basis vector in the M-dimensional vector space and G "+ de"ne the operator e+
0 ,\"+ e+ J " 0 , + , ,\"+ +,"+ $ e+ +
0 ,\"+ where 0 is the (N!1);M-dimensional zero matrix. J is the operator which performs the ,\"+ + , expansion N! operation on a M-dimensional vector to produce a MN-dimensional vector. In particular, we get J "I , the M;M identity matrix for any M*1. + +"+ In fact, as the signal passes through the N-fold interpolation section, the blocking mechanism is equivalent to blocking the signal to a vector of length M at "rst, then expanding by (N!1) zeros and "nally passing through the block anti-imaging "lter H (z). We conclude the interpolation section with the following block ' mechanism (see Fig. 4): E (z) z\E (z) +,\ JM (z)"H (z)J " z\E (z) + , ' + , +,\ $ z\E (z)
E (z) , E (z) ,\ E (z) ,\ 2 z\E
,>
E (z) , E (z) ,\ E (z) ,\ 2
(z) z\E
,>
2 2 \ \
E (z) +\, E (z) +\,\ E (z) +\,\ $
(z) 2 z\E (z) +\,>
(2.2)
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1507
which is an MN;M transfer matrix of a discrete-time LTI MIMO system with M inputs and MN outputs. Remark. In particular, for FIR h (z) with length ¸"MN, ' h(0) h(1) h(2) 2 h(MN!1) z\h(MN!1)
h(0)
H (z)" z\h(MN!2) z\h(MN!1) ' $ $ z\h(1)
z\h(2)
h(1)
\
h(MN!2)
h(0)
\
h(MN!3)
$
\
.
$
z\h(3) 2
h(0)
In this situation, (2.2) becomes h(0)
h(N)
h(2N)
2
h((M!1)N)
z\h(MN!1)
h(N!1)
h(2N!1)
\
h((M!1)N!1)
z\h(MN!2) JM " + , z\h(MN!3)
h(N!2)
h(2N!2)
\
h((M!1)N!2)
h(N!3)
h(2N!3)
\
$
$
$
$
\
$
z\h(1)
z\h(N#1)
z\h(2N!1) 2 z\h((M!1)N#1)
which is a MN;M transfer matrix of a discrete-time MIMO system with M inputs and MN outputs. The compression section can be similarly implemented. In the equivalent structure of interpolation section shown in Fig. 4, the data sequence is "rst blocked to a vector of length M, s(nM) s(n)"
s(nM#1) $ s(nM#M!1)
and then passes through the operator J .Then the data sequence passes through the equivalent anti+ , imaging block "lter H (z) of h (z). Any of the techniques introduced in the traditional DSP text books such as ' ' [4,10] can be applied for the design of the lowpass "lter. The z-transform of the output vector at the interpolation section is s (z)"H (z)J s(z) ' ' + , "JM (z)s(z), + , where s (z) and s(z) is the z-transform of the vector sequences s (n) and s(n), respectively. ' ' The kernel section in Fig. 3 consists of three parts: the chirp multiplication, impulse response h(n) and another chirp multiplication. From the perspective of the MN-input vector s (n), Ks (n) is an output of an ' ' LTI MIMO system with transfer matrix K. The convolution transfer function h(z) can also be converted to the equivalent block "ltering form H(z) which is an MIMO system with MN inputs and MN outputs. Hence, the z-transform of the output vector of the kernel section is (see Fig. 5) s?(z)"KH(z)Ks (z)"K ?(z)s (z), ' ' ' where K?(z)"KH(z)K is the transfer matrix of the kernel section.
(2.3)
1508
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
Fig. 5. The equivalent MIMO system of DFrFT in Fig. 3.
Fig. 6. The general scheme of the DFrFT blocking system.
Similarly, the decimation section can be blocked as an equivalent MN;M MIMO transfer matrix, which can be expressed as DM (z)"D H (z), (2.4) + , + , " where D "J2 is a matrix which performs the N-fold decimation operation on the MN-dimensional + , + , vector s?(z) and H (z) is an equivalent block version of the anti-aliasing "lter h (z). Note that the superscript ' " " T denotes the matrix transpose. Combining (2.2)}(2.4), we obtain an equivalent block "ltering system with the following transfer matrix: H ?(z)"DM
(z)K ?(z)JM (z) (2.5) + , + , to compute the DFrFT of a given signal sequence s(n) as shown in Fig. 5. This means that the scale input/output DFrFT "ltering system can be implemented as a parallel MIMO DFrFT "ltering system (2.5). Note that H ?(z) is a transfer matrix with M-input/M-output. The structure of Fig. 6 shows the speed advantage of the blocking. In general, the sampling rate of the input signal s(n) can be M times faster than the speed of the original scale input system. However, since H ?(z) is an M;M system, it requires a larger number of computational units (multipliers and adders) than the original scale input/output system. That is, the computational parallelism is increased in the blocked implementation by increasing the number of computational units. As a result the signals are processed with M times higher rate than the rate that can normally be handled by one computational unit.
3. The DFrFT polyphase system In this section, the proposed general multirate "ltering system (2.5) is applied to the implementation of the known DFrFT algorithms (see Fig. 6). The algorithms involving multirate operations are introduced "rst,
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1509
and the algorithms without multirate operations are then shown to be the special case of the proposed system (2.5). The details are given in the following subsections in order. 3.1. The multirate DFrFT polyphase xltering system With the similar analysis in the previous section, performing the 2! operation is equivalent to performing the operator 1
0
0
$
0
0
0
0
0
$
0
0
0 1 0 J " + $ 2 2
\
\ 0
\
\ $
0
0
0
\
0
1
0
0
0
2
0
0 +"+
and 2 operation is equivalent to performing operator D "J2 . By letting K ?(z)"KK ?K and after some + + manipulations, algorithm (1.1) now can be implemented by an equivalent MIMO signal processing system with the following transfer matrix: DM (z)K? (z)JM (z) a 3[!1.5,!0.5]6[0.5,1.5], + + DM (z)K?\(z)JM (z)= otherwise, + +
H?(z)"
"
D H (z)K? (z)H (z)J a 3[!1.5,!0.5]6[0.5,1.5], + " ' + D H (z)K?\(z)H (z)J = otherwise, + " ' +
(3.1)
where a "(a mod 4)!2, H (z) and H (z) are the block versions of the anti-aliasing "lter h (z) and " ' " anti-imaging "lter h (z), respectively. Similarly, algorithm (1.2) can be implemented by an MIMO system with ' transfer matrix H?(z)"DM (z)KH(z)KJM (z), !1)a)1, + +
(3.2)
where H(z) is the block "lter of the transfer function h(z) (see (1.3)). 3.2. The DFrFT polyphase xltering system without multirate operations If the signal is sampled at twice the Nyquist rate [11], there is no need to perform the multirate operations in algorithms (1.1) and (1.2). In this situation, the DFrFT algorithms (1.1), (1.2), together with (1.5) and (1.6) all involve a matrix}vector multiplication. This means that these algorithms can be implemented by a block "ltering system, hence referred to as an MIMO system with order 0 transfer matrix K ?(z). Therefore, N is set to be 1 in the design of the DFrFT transfer matrix (2.5). Further, the anti-imaging "lter h (z),1 and the ' anti-aliasing "lter h (z),1, hence H (z)"H (z),I, therefore JM (z),I,DM (z). The system in Fig. 3 is " " ' + + now simpli"ed to a system consisting of block mechanism, kernel section and unblocking mechanism shown in Fig. 6. The system with transfer matrix (3.2) is now simpli"ed as H ?(z)"KH(z)K, !1)a)1
(3.3)
1510
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
and the transfer matrix for algorithm (3.1) is K ? (z)
K ?(z)"
a 3[!1.5,!0.5][0.5,1.5],
(3.4)
K?\(z)W otherwise.
The transfer matrix for the rotational Fourier transform (1.5) is K ?(z)" p (a)W G, G G
0)a)2p
(3.5)
and the DFrFT algorithm (1.6) can be implemented by the transfer matrix K ?(z)"UK D? pUK 2, 0)a)2p.
(3.6)
Remark. Note that (3.3) and (3.4) do not obey the additivity of rotations [14], i.e. H?>@(z)"H ?(z)H @(z).
(3.7)
4. Application to fractional Fourier 5lter banks One of the most important classes of "lter banks is the DFT "lter bank [4,16]. In the previous section, we have shown that the above DFrFT algorithms can be implemented by the MIMO transfer matrix H ?(z). In this section, we apply the polyphase systems to a generalization of the DFT "lter bank with parameter a. It is worthwhile pointing out that the block mechanism can be implemented by a delay chain followed by an M-compressor as shown in Fig. 7. Hence, the DFrFT system can be converted to an equivalent M-channel maximally decimated "lter bank system by "rst rearranging the M-compressor and the transfer matrix by noble identity (see Fig. 8) [16]. This DFrFT system can be further converted to the following equivalent "lter bank:
h? (z) h? (z) $
h? (z) +\
z\+\ z\+\
"H?(z+)
$
"
1
H? (z+) H? (z+) $
H? (z+) H? (z+) $
2 2 \
H? (z+) H? (z+) 2 H? (z+) +\ +\ +\ +\
where h? (z) is the transfer function of the kth analysis "lter de"ned as I +\ h? (z)" z\+\\JH? (z+), I I J J
H? (z+) +\ H? (z+) +\ $
k"0,2, M!1
z\+\ z\+\ $
,
1
(4.1)
in which a denotes the fractional parameter of the fractional Fourier "lter bank (aFrFB). Based on di!erent DFrFT polyphase "ltering systems, we can construct the corresponding aFrFB systems. Because the computation of the DFrFT is valid for all the above algorithms, the aFrFB construction is valid as an extension of the DFT "lter bank system.
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1511
Fig. 7. The fractional Fourier "lter bank (aFrFB) system.
Fig. 8. The noble identity for multirate operations.
In particular, when kernel (3.5) is chosen with a"p/2, the aFrFB is the same as the conventional DFT "lter bank which is known as the spectrum analyzer. It provides a unique representation for signal in terms of complex exponentials. Further, when kernel (3.6) is chosen with a"p/2, an equivalent DFT "lter bank is obtained by aFrFB, but expressing the signal in terms of Hermite functions. Apparently, DFrFT systems (1.5) and (3.6) satisfy the unitary property H ?(z)H ?(z)H"H ?(z)H\?(z)"cz\JI,
(4.2)
where superscript * denotes the complex conjugate, c is a constant and l is the system delay (possibly equal to 0). We can easily design a perfect reconstruction aFrFB system with analysis "lter bank (4.1) and synthesis "lter bank (see Fig. 9) with the following transfer functions +\ g? (z)" z\JH\?(z+), k"0,2, M!1. I I J J 5. Conclusion In this paper, we propose a special multirate "ltering system with blocking implementation for fast computation of DFrFT. With the blocking processing we have shown that the DFrFT can be implemented
1512
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
Fig. 9. The scheme of an aFrFB with perfect reconstruction property.
by an MIMO system. As is well known, the structure provides a rather e$cient way to improve the computation of DFrFT. The cost paid for the e$ciency of computation is that the components of the structure are more complex than that of the direct implementation. This blocking structure can also be implemented by an equivalent fractional "lter bank structure. Accordingly, a fractional Fourier "lter bank can be constructed by a corresponding DFrFT polyphase system. The proposed approach can be constructed by known digital signal processing techniques. Therefore, it is more suitable for practical applications.
References [1] L.B. Almeida, The fractional Fourier transform and time-frequency representations, IEEE Trans. Signal Process. 42 (1994) 3084}3091. [2] C.S. Burrus, Block implementation of digital "lters, IEEE Trans. Circuits Theory CT-18 (1971) 697}701. [3] C.S. Burrus, Block realization of digital "lters, IEEE Trans. Audio Electroacoust. AU-20 (1972) 230}235. [4] R.E. Crochiere, L.R. Rabiner, Multirate Digital Signal Processing, Prentice-Hall, Englewood Cli!s, NJ, 1983. [5] A.W. Lohmann, D. Mendlovic, Z. Zalevsky, Fractional transformations in optics, in: Process in Optics, Vol. XXXVIII, Elsevier, Amsterdam, 1998 (Chapter 4). [6] A.C. McBride, F.H. Kerr, On Namias's fractional Fourier transforms, IMA J. Appl. Math. 39 (1987) 159}175. [7] D. Mendlovic, H.M. Ozaktas, Fractional Fourier transforms and their optical implementation: I, J. Opt. Soc. Am. A 10 (9) (1993) 1875}1881. [8] S.K. Mitra, R. Gnanasekaran, Block implementation of recursive digital "lters } new strucrures and properties, IEEE Trans. Circuits Systems CAS-25 (1978) 200}207. [9] V. Namias, The fractional order Fourier transform and its application to quantum mechanics, J. Inst. Math. Appl. 25 (1980) 241}265. [10] A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, Englewood Cli!s, NJ, 1989, pp. 101}111 (Chapter 3). [11] H.M. Ozaktas, O. Arikan, M.A. Kutay, Gozde Bozdagi, Digital computation of the fractional Fourier transform, IEEE Trans. Signal Process. 44 (9) (1996) 2141}2150. [12] H.M. Ozaktas, B. Barshan, D. Mendlovic, L. Onural, Comvolution, "ltering, and multiplexing in fractional Fourier doamins and their relation to chirp and wavelets, J. Opt. Soc. Am. A 11 (1994) 547}559.
D.-F. Huang, B.-S. Chen / Signal Processing 80 (2000) 1501}1513
1513
[13] H.M. Ozaktas, M.A. Kutay, D. Mendlovic, Introduction to fractional Fourier transform and applications, in: Advances in Imaging and Electron Physics, Vol. 106, Academic Press, New York, 1999, pp. 263}342 (Chapter 4). [14] S.C. Pei, M.H. Yeh, C.C. Tseng, Discrete fractional Fourier transform based on orthogonal projections, IEEE Trans. Signal Process. 47 (5) (1999) 1335}1348. [15] B. Santhanam, J.M. McCellan, The discrete rotational Fourier transform, IEEE Trans. Signal Process. 44 (4) (1996) 994}998. [16] P.P. Vadyianthan, Multirate Systems and Filter Banks, Prentice-Hall, Englewood Cli!s, NJ, 1993. [17] P.P. Vaidyanathan, S.K. Mitra, Polyphase networks, block digital "ltering, LPTV systems, and alias-free QMF banks: a uni"ed approach based on pseudocirculants, IEEE Trans. Acoust. Speech Signal Process. 36 (3) (March, 1988) 381}391.