Xiaofeng Liu

Youjian Zhao
email: zhaoyoujian@tsinghua.edu.cn

Yajuan Wu
email: scwuyajuan@yahoo.com.cn

An Efficient Crosstalk-Free Routing Algorithm Based on Permutation Decomposition for Optical Multi-log 2 N Switching Networks

Keywords: Permutation, Optical switching, Multi-log2N networks, Optical crosstalk

Optical switching networks (OSN) based on optical directional couplers (DC) may be the most promising candidate to provide a high switching rate when the speed mismatch problem between links (optical fibers) and switches is increasingly serious. Although such switches have many advantages, the DC suffers from an inherent crosstalk problem that can greatly aggravate the switching performance. Based on semipermutations, a parallel decomposition algorithm,which is called multidecomposition, is proposed in this paper for solving the optical crosstalk problem in optical multi-log2N switching networks. According to the number of planes in a multi-log2N network, the multi-decomposition is performed in parallel to partition a permutation into several subpermutations, each of which is established without crosstalk within each plane. We demonstrate that our algorithm can completely remove the crosstalk in optical multi-log2N networks when n is even, and that it may be generated only in the stage (n-1)/2 (i.e., the middle stage) when n is odd, but the corresponding probability of generating crosstalk is to be less than or equal to 1 2 (n+1)/2 -1 . In addition, our algorithm can achieve a low complexity for decomposition a permutation due to its parallelism so that any permutations can be realized in multi-log2N networks under the constraint of avoiding crosstalk.

Introduction

The Internet is an important product of the information age. From a high-level perspective, the entire Internet architecture consists of two parts: communication links and switching nodes. At present, the capacity of these two parts has an enormous difference. The speed of communication links has been drastically increased with the advent of dense wavelength-division multiplexing (DWDM) technologies, but the progress made in the switches has lagged relatively behind. In order to solve the speed mismatch problem, the optical switches have been widely explored in recent years.

A multi-log 2 N network [START_REF] Lea | Multi-log2N networks and their applications in high-speed electronic and photonic switching systems[END_REF], composed of several single-log 2 N networks, possesses many characteristics which are helpful for photonic switching system. As an optical directional-coupler (DC) is introduced into this switching network to replace the old (electronic) switching element (SE), the transmission rate of signals can achieve several T bps if the state (cross or bar) of the DC has been set up properly. Meanwhile, A blocking, which is called crosstalk-blocking [START_REF] Vaez | Strictly nonblocking directional-coupler-based switching networks under crosstalk constraint[END_REF] and differs from link-blocking, is introduced into this switching system by DCs. This blocking will occur when two light signals pass through the same DC at the same time. The crosstalk-blocking limits the scalability of switching networks so that it is not easy to use DCs to construct a larger switching network [START_REF] Hinton | An introduction to photonic switching fabrics[END_REF] because two signals passing through the same DC interact with each other. Therefore, the crosstalk-blocking is an important factor to affect the quality of communication in an optical switching network (OSN). Eliminating the crosstalk-blocking in an OSN is just the main objective of this paper.

In this paper, we propose an efficient algorithm to eliminate the crosstalkblocking (abbreviated as crosstalk) of optical multi-log 2 N switching networks. The central idea of the proposed algorithm is based on the concept of a semipermutation in the literature [START_REF] Yang | Permutation capability of optical multistage interconnection networks[END_REF], but we extend this concept to multiple decomposition of a permutation for the crosstalk-free routing in optical multilog 2 N switching networks, i.e., a permutation is first partitioned into two semipermutations, and is further divided into four quarter-permutations and so on. By so doing, any two different optical signals traversing down two node-disjoint (or DC-disjoint) paths cannot generate the optical crosstalk problem within each plane. Thus this multiple decomposition algorithm is named multi-decomposition in the remainder of the paper.

This paper is followed by four sections. Section 2 describes the crosstalk and its common solutions. A basic network model and all related preliminaries associated with this work are illustrated in section 3. Section 4 presents the multi-decomposition in detail and then conclusions are given in section 5.

Crosstalk and Related Researches

Crosstalk and signal attenuation are the major problems that have been hindering the development of OSNs all the time. The signal attenuation problem can be solved by a semiconductor optical amplifier (SOA). However, the crosstalk cannot be removed by the SOA because the crosstalk is also amplified when the desirable signal is amplified. Thus we must find the other methods to remove the crosstalk in an OSN. Indeed,it is almost impossible to eliminate the optical crosstalk totally unless only one optical signal passes through a DC at any given time, as shown in Fig. 1. Traversing ways in a DC Fig. 1. Both space domain approach [START_REF] Padmanabhan | dilated networks for photonic switching[END_REF] and time domain approach [START_REF] Qiao | A time domain approach for avoiding crosstalk in optical blocking multistage interconnection networks[END_REF] are based on this idea. Essentially, the former adds the number of DCs, while at most only one of the two inputs and outputs is active at a time in a DC. Therefore, hardware cost has been sacrificed to trade the crosstalk-free routing. For the latter, a permutation is decomposed into several sub-permutations such that all connections of each sub-permutation can be established simultaneously without crosstalk, i.e., each time slot is required to route each sub-permutation and all the connections corresponding the permutation are established within several time slots in a time-division multiplexed fashion. As such, this approach uses time cost to exchange the crosstalk-free routing. In addition to these two approaches above, wavelength dilation [START_REF] Sharony | The wavelength dilation concept in lightwave networks -Implementation and system considerations[END_REF] is the other one that eliminates the crosstalk in an OSN. A technology of wavelength grouping method (WGM) [START_REF] Sharony | The wavelength dilation concept in lightwave networks -Implementation and system considerations[END_REF] is proposed in this approach, in which the wavelengths are partitioned into several groups so that the wavelengths in each group are widely separated, with nearby wavelengths placed in different groups, i.e., the wavelengths are selected should be far enough apart (a few nanometers) so as not to interact with each other. A drawback of this approach is that a specific structure is required to select and separate the wavelength, and then result in increasing the hardware cost. This method is adopted in [START_REF] Zheng | A parallel self-routing rearrangeable nonblocking multi-log2N photonic switching networks[END_REF] [START_REF] Wong | Crosstalk reduction through wavelength assignment in WDM photonic switching networks[END_REF] to remove the crosstalk of optical multi-log 2 N networks.

The existing methods above removing the crosstalk focus usually on the single-log 2 N network ( [START_REF] Wong | Crosstalk reduction through wavelength assignment in WDM photonic switching networks[END_REF] is an exception), and they sacrifice either time cost or hardware cost. Moreover, due to the unique path property of a log 2 N network, some sub-permutations can be realized within a log 2 N network in a single pass whereas the others cannot [START_REF] Yang | Permutation capability of optical multistage interconnection networks[END_REF]. In fact, the spirit of all of these approaches mentioned above is to avoid two light signals with the same wavelength passing through a common DC at the same time, and the proposed approach in this paper is no exception. However, our algorithm uses the idea of multiple decomposition of any permutation to realize the crosstalk-free routing within each plane of a multi-log 2 N network and does not increase time cost and hardware cost for optical multi-log 2 N switching networks.

Basic Network Model and Preliminaries

Multi-log 2 N Network Model

Multi-log 2 N networks are vertically stacked with multiple log 2 N networks, N demultiplexers and N multiplexers. The log 2 N network has been composed of N (= 2n) inputs and outputs and n(= log 2 N ) stages. Each stage consists of N /2 2 × 2 DCs and any adjacent stages are connected by N interstage fiber links. Planes are vertically stacked to N demultiplexers (resp. multiplexers) in input (resp. output) stage. N denotes the number of source inputs and destination outputs labeled by 0, 1, • • • , N -1 from top to bottom, n(= log 2 N ) denotes the number of stages numbered by 0, 1, • • • , n -1 from left to right, and m indicates the number of planes contained in a multi-log 2 N network. Since all log 2 N networks including baseline, omega, and banyan-type [START_REF] Goke | Banyan networks for partitioning multiprocessor systems[END_REF] have the topologically equivalent feature [START_REF] Wu | On a class of multistage interconnection networks[END_REF], we use an N × N (N = 2n) baseline network as the representative of routing planes in our work. An example of an 8 × 8 baseline network is shown in Fig. 2, and the corresponding multi-log 2 N network is illustrated in Fig. 3.

(0)000 Stage 0 1 2
(1)001 ( 2)010 ( 3)011 ( 4)100 ( 5)101 ( 6)110 [START_REF] Sharony | The wavelength dilation concept in lightwave networks -Implementation and system considerations[END_REF]111 000(0) 001( 1) 010( 2) 011( 3) 100( 4) 101( 5) 110( 6) 111( 7)

Related Preliminaries

For a log 2 N network, let u and v be any two SEs in stage i, S j (u) and S j (v) are the two sets of SEs to which u and v can reach in stage

(i + j) (0 ≤ i < n -1, 1 ≤ j ≤ n -i -1)
. Sets S j (u) and S j (v) certainly can satisfy one of the following properties:S j (u) ∩ S j (v) = ∅ or S j (u) = S j (v) [START_REF] Zheng | A parallel self-routing rearrangeable nonblocking multi-log2N photonic switching networks[END_REF]. For S j (u) and S j (v),the equality S j (u) = S j (v) holds, which implies that SEs u and v are sure to share the same SE in stage (i + j), hence light signals passing through SEs u and v must generate the crosstalk phenomenon at the shared SE in stage (i + j). Therefore, we call SEs u and v possessing (i, j)-buddy if the equality S j (u) = S j (v) holds.

A theorem can be obtained immediately from the description of (i, j)-buddy above.

Theorem 1. Let u and v be two different DCs in stage i. Two optical signals going through u and v do not generate crosstalk in stage

(i+j) (0 ≤ i < n-1, 1 ≤ j ≤ n -i -1) if u and v do not have the (i, j)-buddy property.
Proof. if u and v don't have the (i, j)-buddy property, then optical signals going through u and v do not share any SEs in stage (i + j), so there is no reason to occur crosstalk in stage (i + j).

⊓ ⊔

Multi-decomposition Algorithm of a Permutation

Decomposition of a Permutation

A permutation is usually adopted to describe a mapping between inputs and outputs for a switching network. For an N × N network, suppose there is a permutation P which maps input x i to output y i , i.e., P

(x i ) = y i where x i , y i ∈ {0, 1, • • • , N -1} for 0 ≤ i < N -1.
We use the representation as equation ( 1) to denote this permutation.

P = ( x 0 x 1 • • • x N -1 y 0 y 1 • • • y N -1 ) (1)
The semi-permutation [START_REF] Yang | Permutation capability of optical multistage interconnection networks[END_REF] is the basis of the multi-decomposition algorithm, so we proceed with the concept of semi-permutation first.

Definition 1 ( [4]). For any permutation

P of {0, 1, • • • , N -1}, a partial per- mutation ( x i1 x i2 • • • x i N 2 y i1 y i2 • • • y i N 2 ) is referred to as a semi-permutation, if {⌊ xi 1 2 ⌋, ⌊ xi 2 2 ⌋, • • • , ⌊ xi N/2 2 ⌋} = {0, 1, • • • , N/2 -1} and {⌊ yi 1 2 ⌋, ⌊ yi 2 2 ⌋, • • • , ⌊ yi N/2 2 ⌋} = {0, 1, • • • , N/2 -1}. Since 2(k -1) ≤ x i k , y i k ≤ 2k -1, k = 1, 2, • • • , N/2,
a semi-permutation ensures that the crosstalk can be removed in the first stage and last stage of a log 2 N network, but the crosstalk cannot be always removed in the intermediate stages. Besides, this approach needs several passes to establish a permutation in a single-log 2 N network, and hence the routing time is increased.

To eliminate the crosstalk of each stage in an optical multi-log 2 N switching network, we give the definition of the multi-decomposition based on semipermutation (i.e., Definition 1) as follows.

Definition 2. For any permutation

P of {0, 1, • • • , N -1}, let ( x i1 x i2 • • • x i N m y i1 y i2 • • • y i N m )
be a partial permutation of P, where m is an integral power of 2. {⌊

xi 1 m ⌋, ⌊ xi 2 m ⌋, • • • , ⌊ xi N/m m ⌋} = {0, 1, • • • , N/m -1} and {⌊ yi 1 m ⌋, ⌊ yi 2 m ⌋, • • • , ⌊ yi N/m m ⌋} = {0, 1, • • • , N/m -1}
are used to partition this partial permutation such that the resulting sub-permutation is smaller than semi-permutation in size.

This decomposition method is referred to as multi-decomposition for the multiple decompositions. The resulting partial permutation of multi-decomposition of a permutation is still called sub-permutation from now on.

Note that m denotes the number of sub-permutations into which the full permutation is divided, and each sub-permutation contains N/m connections of the full permutation. An example is given to understand this partition process as follows. )

Step 1. Run the multi-decomposition algorithm (m = 2) to get two semi-permutations P 1 and P 2 . )

P 2 =
( 1 2 4 7 9 11 13 14 0 9 6 13 2 4 11 14

)
Step 2. For all semi-permutations, run the multi-decomposition algorithm (m = 4) again in parallel to get four quarter-permutations, P 11 , P

)
Step 3. Multi-decomposition ends when the number of sub-permutations is m, or else goes to step 2. Since m is equal to 4 in this example, the multidecomposition algorithm should be stop.

Fig. 4 illustrates the switching of these four quarter-permutations in a multilog 2 N network (N = 16, m = 4).

Now, we answer the following two questions about the multi-decomposition. The first one is how many times the multi-decomposition algorithm should be performed when a full permutation needs to be decomposed, i.e., what is the proper value of m? The other is the effectiveness of this algorithm. The proper value of m is determined by the number of copies in a multi-log 2 N network. Since the rearrangeable nonblocking (RNB) network is considered in our work, the number of copies which is needed to build a rearrangeable nonblocking multilog 2 N network is 2 ⌊n/2⌋ [START_REF] Lea | Multi-log2N networks and their applications in high-speed electronic and photonic switching systems[END_REF]. Therefore, the value of m is taken 2 ⌊n/2⌋ in this paper. For the second one, we use the following theorem to give an answer.

Theorem 2. For a multi-log 2 N network, any permutation P can be partitioned into m sub-permutations, where n = log 2 N and m = 2 ⌊ n 2 ⌋ .

Proof. The multi-decomposition is implemented by an iterative manner. A permutation is first partitioned into two semi-permutations, and then each semipermutation is further partitioned into its own two semi-permutations and so on. This partitioning process cannot stop until the number of sub-permutations is equal to m. According to this partitioning process, this theorem is essentially equivalent to a proposition that "Any permutation can be decomposed into two semi-permutations." This proposition can be easily proved by using the P. Hall's distinct system representative [START_REF] Hall | On representatives of subsets[END_REF]. The detailed proof can be found in [START_REF] Yang | Permutation capability of optical multistage interconnection networks[END_REF], so we can use a method similar to one in [START_REF] Yang | Permutation capability of optical multistage interconnection networks[END_REF] to prove this theorem. Due to space limitation, the detailed proof is omitted here.

⊓ ⊔

Theorem 2 guarantees the feasibility of the multi-decomposition of permutations. Next, we prove that the multi-decomposition is effective in removing the crosstalk problem. We first give the following concepts of input set (IS) and output set (OS). Definition 3. For an N × N (N = 2n) multi-log 2 N network, there are n(= log 2 N ) stages.

(1) If n is even, the N inputs (resp. outputs) of the multi-log 2 N network are divided equally into i = 2 n/2 input (resp. output) sets from top to bottom. These sets are referred to as

I 0 , I 1 , • • • , I i-1 and O 0 , O 1 , • • • , O i-1 , respectively.
As shown in Fig. 5(a).

(2) If n is odd, the N inputs (resp. outputs) of the multi-log 2 N network are divided equally into i = 2 (n+1)/2 input (resp. output) sets from top to bottom. These sets are also referred to as Note that the order of input and output ports in each set is consecutive and ascending, and that the number of ports is equal. An important fact obtained from Definition 2 and 3 is that different connections of a sub-permutation always start from different ISs and end in different OSs. Thus every sub-permutation does not have any crosstalk in the input and output stage. The following theorem tells us that the optical crosstalk does not occur in the intermediate stages, either.

I 0 , I 1 , • • • , I i-1 and O 0 , O 1 , • • • , O i-1 ,
Theorem 3. For an N × N (N = 2n) multi-log 2 N network, n = log 2 N , and m = 2 ⌊n/2⌋ . Multi-decomposition guarantees that [START_REF] Lea | Multi-log2N networks and their applications in high-speed electronic and photonic switching systems[END_REF] The crosstalk never occurs in the optical multi-log 2 N network when n is even;

(2) The crosstalk may be possible only in the center stage (i.e., stage-(n-1)/2) when n is odd, but the probability of generating crosstalk is less than or equal to

1 2 (n+1)/2 -1 .
Proof.

(1) M is equal to 2 n/2 if n is even. According to Definition 3, N input and output ports are divided evenly into ISs and OSs, and every sub-permutation contains N/m(= 2 n/2 ) connections. As stated previously, the crosstalk is never generated in the input stage and output stage (by Definition 2 and 3). We now prove that the crosstalk does not occur in the intermediate stages. We first consider the first n/2 stages. If the crosstalk has occurred in stage j (1 ≤ j ≤ n 2 -1) , there are at least two different DCs u and v in the same IS of a plane satisfying the (0, j)-buddy (by Theorem 1), but only one connection in a sub-permutation starts from an IS within each plane at any time. This contradiction proves our conclusion that no crosstalk occurs in the first n/2 stages. In the last n/2 stages, the crosstalk can also be avoided, which can be proved by contradiction. Once the crosstalk is generated, some connections in the same plane will be routed to the same OS, but this case is also impossible. In fact, these two cases are entirely symmetrical.

(2) If n is odd, m = 2 (n-1)/2 . The crosstalk does not occur in the first (n -1)/2 stages and last (n -1)/2 stages, but it possibly emerges in the center stage. Take the Fig. 5(b) as an example, two connections (2 → 3) and (7 → 5) share a common DC at the middle stage, now the crosstalk has occurred, but no DC can be shared in the first (n-1)/2 stages and last (n-1)/2 stages. Following a similar argument to the first case (n is even), we can prove that when n is odd, the crosstalk is never generated in the first (n -1)/2 stages and last (n -1)/2 stages. Here we fix our attention on the proof of the probability of generating crosstalk in the middle stage.

⊓ ⊔

For brevity and clarity, we introduce the following definitions of buddy input set (BIS) and buddy output set (BOS). The set BIS is comprised of all input sets I k sharing a common DC at the center stage, and the set BOS is comprised of all output sets O k sharing the same DC at the center stage. For example, I 0 and I 1 belong to BIS 0 because they share the same DC a in the center stage, i.e., BIS 0 = {I 0 , I 1 }; similarly, BOS 0 = {O 0 , O 1 }. As shown in Fig. 6. According to the previous definitions, some good characteristics of the multidecomposition are summarized as follows:

0
I 0 I 1 I 2 I 3 I 4 I 5 I 6 I 7 O 0 O 1 O 2 O 3 O 4 O 5 O 6 O 7 Stage 0 1 2 3 4 BIS 0 BIS 1 BIS 2 BIS 3 0 1 2 3
-It is impossible that two (or more) connections in each plane start from the same IS and head for the same OS (by Definition 3); -Two connections sharing the same DC in the center stage should come from the same BIS and head for the same BOS; -According to the distinct BISs, all DCs in the center stage are partitioned into i/2 sections (as the dashed boxes in Fig. 6) such that at most one DC possibly generates the crosstalk within each section at any time.

These characteristics imply that all DCs in the same section are mutually exclusive in generating crosstalk, and that the corresponding DCs among distinct sections are mutually exclusive in generating crosstalk as well. Thus, these i/2 sections possess the identical statistical nature, we may take the first DC a (the shadow box in Fig. 6) in the first section as the representative to discuss the probability of generating crosstalk.

Any inputs in each IS can be connected to arbitrary OS with equivalent probability, i.e., P {I k → O j } = 1/i(0 ≤ k, j < i), so the probability of generating crosstalk in SE a is

P crosstalk (a) = P {(I 0 → O 0 , I 1 → O 1 ) ∪ (I 0 → O 1 , I 1 → O 0 )} = P {(I 0 → O 0 , I 1 → O 1 ) + P {(I 0 → O 1 , I 1 → O 0 )} = P {I 0 → O 0 }P {I 1 → O 1 |I 0 → O 0 } + P {I 0 → O 1 }P {I 1 → O 0 |I 0 → O 1 } = 1 i × 1 i -1 + 1 i × 1 i -1 = 2 i(i -1) (2)
On the other hand, at most one DC generates crosstalk possibly in each section in any given time. Then the probability of generating crosstalk in each section, which is denoted by P crosstalk (section), is equal to P crosstalk (a). Therefore the probability of generating crosstalk in entirely plane is

P crosstalk (plane) = P {∪ i/2 k=0 P crosstalk (section k )} ≤ ∪ i/2 k=0 P crosstalk (section k ) = ∪ i/2 k=0 P crosstalk (a) = 1 i -1 (3)
Now, substituting i = 2 (n+1)/2 into expression (3), we obtain P crosstalk (plane) ≤ 1 2 (n+1)/2 -1 .

Implementation of the Multi-decomposition and Its Analysis

In this subsection, we will discuss the implementation of the multi-decomposition algorithm of permutations and its time complexity.

There is a permutation as the equation ( 1), in which input x i is mapped to output y i (0 ≤ i ≤ N -1). Based on the semi-permutation [START_REF] Yang | Permutation capability of optical multistage interconnection networks[END_REF], we design the decomposition algorithm of a permutation as Algorithm 1 in order to design our multi-decomposition algorithm. The Algorithm 1 contains three parameters that are P, N and k. These parameters denote the permutation to be decomposed, the number of requests and the decomposition levels, respectively.

Algorithm 1: Dichotomy-of-permutation (P, N, k) Input: Any permutation P contains N request pairs. Output: Two sub-permutations P1 and P2 3. Traverse the bipartite graph and color the two adjacent edges of the same vertex with different colors, then all edges with the same color are grouped into forming a sub-permutation. Since the chromatic number of the bipartite graph is two [START_REF] Yang | Permutation capability of optical multistage interconnection networks[END_REF], the permutation can be partitioned into two sub-permutations P1 and P2.

Two vertex sets VI

= {Ai 0 , Ai 1 , • • • , Ai N/2 k -1 } and VO = {Ao 0 , Ao 1 , • • • , Ao N/2 k -1 } are built from the per- mutation P, where Ai j = {x j•2 k , x j•2 k +1 , • • • , x j•2 k +2 k -1 }, AO j = {y j•2 k , y j•2 k +1 , • • • , y j•2 k +2 k -1 }, 0 ≤ j < N/2 k
The multi-decomposition algorithm can be implemented by calling Algorithm 1 repeatedly. Let m be the integral power of 2, i.e., m = 2 k (1 ≤ k ≤ n). The multi-decomposition algorithm becomes semi-permutation decomposition in [START_REF] Yang | Permutation capability of optical multistage interconnection networks[END_REF] when k is 1. One exactly connection is established in each plane when k equals to n, thus the crosstalk and blocking can be avoided undoubtedly. In the other cases, Algorithm 1 is called repeatedly k rounds and is executed 2 i-1 (1 ≤ i ≤ k) times in parallel in each round. The corresponding algorithm is demonstrated as Algorithm 2. As the Example 1 mentioned earlier, N = 16, m = 4(= 2 2 ), and the Algorithm 1 has been called 2 rounds. Two semi-permutations are obtained after the first round, and four quarter-permutations are got after the second round.

A permutation is partitioned into m sub-permutations, and the decomposition process is similar to that of building a binary tree with m sub-permutations as leaves. Furthermore, this process is carried out in parallel, so the time complexity of Algorithm 2 is O(log 2 m). On the other hand, the Algorithm1 has Algorithm 2. multi-decomposition of permutations Input: N requests contained by a permutation P ; Output: M sub-permutations, m = 2 k ,k is a natural number. For i = 1 to k do For j = 1 to 2 i-1 do in parallel {P [j, 1], P [j, 2]}=Dichotomy-of-permutation (P [j], N, i); End parallel End for O(N ) time complexity, so the overall complexity of Algorithm 2 is O(N log 2 m). We will further improve Algorithm 1 by parallelism in the future so that the time cost can be reduced much.

Conclusions

In this paper, we proposed an efficient algorithm called multi-decomposition to remove the crosstalk in optical multi-log 2 N switching networks. A permutation is partitioned into several sub-permutations by using our decomposition algorithm, which ensures that each sub-permutation can be connected without crosstalk within each plane of optical multi-log 2 N networks. Compared with other approaches, our algorithm does not increase the time cost and hardware cost. Although semi-permutation [START_REF] Yang | Permutation capability of optical multistage interconnection networks[END_REF] is the foundation of our algorithm, we have extended the concept of semi-permutation to multiple decomposition of a permutation and successfully solved the crosstalk problem of optical multi-log 2 N switching networks.

We have proved the validity of the multi-decomposition algorithm. The crosstalk can be entirely removed in optical multi-log 2 N networks when n(= logN ) is even. If n is odd, crosstalk may be possible only in the middle stage (i.e., stage-(n-1)/2), but the probability of generating crosstalk is proved to be less than or equal to 1 2 (n+1)/2 -1 . What's more, our algorithm has low time complexity to decompose a permutation due to its parallelism so that any permutations can be routed without crosstalk in an optical multi-log 2 N switching network.

Fig. 2 .Fig. 3 .

 23 Fig. 2. An 8 × 8 baseline network

Example 1 .

 1 For N = 16, m = 4, and a permutation

  (a) Switching of P 11 . (b) Switching of P 21 . (c) Switching of P 22 . (d) Switching of P 12 .

Fig. 4 .

 4 Fig. 4. Switching of four quarter-permutations in a multi-log2N network (N = 16, m = 4)

  respectively. As shown in Fig. 5(b).

Fig. 5 .

 5 Fig. 5. The diagram of input and output sets

Fig. 6 .

 6 Fig. 6. The baseline network when n is odd (n = 5)

2 .

 2 Construct a bipartite graph G = (VI , VO, E) based on VI and VO. A edge e ∈ E is associated with a request pair (xj, yj), where xj ∈ Ai j and yj ∈ Ao j .

  12 , P 21 and P 22 .

	P 11 =	(	0 5 8 12 1 8 15 7	)	, P 12 =	(	5 10 3 12 3 6 10 15	)
	P 21 =	(	1 7 11 13 0 13 4 11	)	, P 22 =	(	2 4 9 14 9 6 2 14

Acknowledgments

. This work is supported by National Natural Science Foundation of China (No.61073167),National High Technology Development Program of China (No.2011AA010704),Projects in the National Science & Technology Pillar Program(2011BAK08B05-02) and Major Training Program of China West Normal University (09A003).


Peng Liu
email: liupeng@zju.edu.cn

Chunming Huang

Jun Guo

Yang Geng

Weidong Wang
email: wdwang@zju.edu.cn

Mei Yang
email: meiyang@unlv.edu

Scalable-Grain Pipeline Parallelization Method for Multi-core Systems

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

With the continuous advance of semiconductor technology, the enormous number of transistors available on a single chip enables the integration of tens or hundreds of processing cores on a multiprocessor system-on-chips (MPSoCs). These processing cores could be homogeneous or heterogeneous, such as processors, digital signal processor cores, memory blocks, etc. To efficiently utilize these parallel resources available on an MPSoC, one challenge is how to parallelize the legacy sequential programs. However, most research and development efforts in MPSoCs are on the hardware architectural side. Research in application program parallelization and parallel programming for MPSoCs is far more behind.

Existing efforts to exploit pipeline parallelism in C programs are mostly finegrained [START_REF] Ottoni | Automatic Thread Extraction with Decoupled Software Pipelining[END_REF][START_REF] Rangan | Decoupled Software Pipelining with the Synchronization Array[END_REF]. Some approaches partition individual instructions across processors, such as the decoupled software pipeline (DSWP) method [START_REF] Ottoni | Automatic Thread Extraction with Decoupled Software Pipelining[END_REF], while others are dedicated to parallelizing scientific and numerical applications, such as DOALL [START_REF] Allen | Optimizing Compilers for Modern Architectures: A Dependence-based Approach[END_REF] and DOACROSS [START_REF] Cytron | Doacross: Beyond Vectorization for Multiprocessors[END_REF]. The HELIX project [START_REF] Campanoni | HE-LIX: Automatic Parallelization of Irregular Programs for Chip Multiprocessing[END_REF] is a generation of the DOACROSS scheme. HELIX satisfies only the necessary loop-carried data dependences. The synchronization required for loop-carried dependences is implemented using a per-thread memory area which resides in the system's shared memory. For streaming applications, these techniques are not sufficient, because the pipeline parallelism in streaming applications is coarse-grained and more complex than that in scientific/numerical applications. To resolve intricate dependency, the fine-grained methods may usually merge those tasks which take part in a dependence cycle. Other work on integration of parallelization techniques [START_REF] Huang | Decoupled Software Pipelining Creates Parallelization Opportunities[END_REF][START_REF] Raman | Parallel-stage Decoupled Software Pipelining[END_REF][START_REF] Vachharajani | Speculative Decoupled Software Pipelining[END_REF] includes the speculation DSWP [START_REF] Vachharajani | Speculative Decoupled Software Pipelining[END_REF], which mainly focuses on cutting dependences caused by loop conditional statements but cannot cut the inter-iteration dependences, and the parallel-stage DSWP [START_REF] Raman | Parallel-stage Decoupled Software Pipelining[END_REF] which integrates DSWP with DOALL. However, these approaches are not suitable for exploiting course-grained parallelism and cannot manage multi-layer loop structures.

Unfortunately, currently no method can solve the problem of extracting coarse-grained pipeline parallelism well [START_REF] Ryoo | Automatic Discovery of Coarse-grained Parallelism in Media Applications[END_REF]. The method to exploiting coarsegrained pipeline parallelism in C programs proposed in [START_REF] Thies | A Practical Approach to Exploiting Coarse-grained Pipeline Parallelism in c Programs[END_REF] is a language-extension approach, which imposes the burden of parallelism extraction on programmers. The Paralax [START_REF] Vandierendonck | The Paralax Infrastructure: Automatic Parallelization with a Helping Hand[END_REF] is a compiler-based parallelization framework which focuses on the dependence analysis but lacks transformation for multi-layer loops. MAPS [START_REF] Ceng | MAPS: An Integrated Framework for MPSoC Application Parallelization[END_REF] is a framework for semi-automatic parallelism extraction from sequential legacy code and extended to support parallel dataflow programming.

Programmers have been familiar with sequential programming languages like C for a long time. Instead of using a language/model, an evolutionary parallelization methodology for sequential codes will be more significant. Streaming applications represent a large set of MPSoC applications, such as video, audio, cryptographic, wireless baseband processing, etc. These programs are characterized by heavy use of pointers, multi-layer loop structures, and streaming data input. The parallelism within these programs usually is implicit [START_REF] Ryoo | Automatic Discovery of Coarse-grained Parallelism in Media Applications[END_REF]. Extracting the pipeline parallelism hidden in loops becomes very critical for the rich loop control-flow constructions in the C programs of streaming applications. Our focus is on exploiting coarse-grained pipeline parallelism in the C programs of streaming applications, which are characterized by multi-layer loop structures with intricate dependence relations and fixed data flow.

In this paper, a scalable-grain pipeline parallelization (SPP) method is proposed, which exploits the coarse-grained pipeline parallelism hidden in multilayer loops existing in streaming applications. These applications are described by a block diagram with a fixed flow of data and the regular communication pattern. We exploit coarse-grained pipeline parallelism by using the source program transformation for embedded applications that can overcome the traditional barriers. The proposed method first resolves intricate dependency by transforming the task dependence graph, which is then scheduled to the target MPSoC to minimize the maximal execution time of a pipeline stage. The experimental results of the parallel programs of the applications on an eight-core platform justify that this method is efficient in parallelizing C sequential programs.

The rest of this paper is organized as follows. The problem to be solved is formulated in Section 2. Section 3 presents the framework of the proposed method. Section 4 presents the experimental result. Section 5 concludes.

Problem Formulation

The MPSoC under consideration is modeled as an architectural characterization graph. For a streaming application, its task dependence graph is extracted from the source code of the application. Definition 1. An MPSoC architectural characterization graph (ACG), ACG = (P, L), is a undirect graph, where a vertex p i ∈ P represents one processing element (PE) in MPSoC, with m(p i ) denoting the memory space of p i and t p (p j ) denoting the type of p i , an arc l ij ∈ L represents the link between p i and p j , with r(l ij ) denoting the link bandwidth. An MPSoC architecture has r(l ij ) = R, ∀l ij ∈ L, and m(p i ) = M , ∀p i ∈ P . Definition 2. A task dependence graph (TDG), T DG = (V, E, W v , W e ), is a directed graph, which represents the dependence relation of a program, where V represents the set of nodes, each representing a task of the program, E represents the set of edges, each representing the dependence relation between two tasks. The sets W v and W e represent the properties of nodes in V and edges in E, respectively. For each n i ∈ V , the predecessor node set of n i is defined as pre(n i ), and the successor node set of n i is defined as succ(n i ). For an edge e ij = (n i , n j ) ∈ E, then n i ∈ pre(n j ) and n j ∈ succ(n i ). There are two types of dependences: data and control dependences. Data dependences are caused by data transfer. Control dependences are evoked by the conditional statements or loop statements.

These are four types of nodes in set V: Control node set V C , Branch node set V B , Loop node set V L , and Ordinary node set V O . The edges in E are sorted into three groups: Control edge set E C , Inter-iteration edge set E I , and Ordinary edge set E O . The properties of W v and W e are listed in Table 1.

Definition 3. The dependence distance d dep indicates the number of iterations between two loop nodes forming the inter-iteration dependences. ∀e ij ∈ E I , n i , n j ∈ V L , if n i is the kth iteration of the loop and n j is the lth iteration of the loop, then d dep (e ij ) = l -k. Definition 4. A strongly-connected component (SCC) of a graph is a maximal strongly connecter sub-graph, in which there exists a path from each vertex to every other vertex in the sub-graph. A directed graph is acyclic if and only if it has no SCC. Given two nodes n i and n j in a TDG, the SCC distance d s (n i , n j ) indicates the number of SCCs between SCCs. Definition 5. A TDG after transformation, T DG = (V , E , W v , E e ), is an acyclic directed graph, in which each SCC in the TDG is merged into a single vertex and the control dependences are removed. V is the union of vertex sets V B , V L , and V O . The edge set E does not include any branch/loop control edge, which is the union of edge sets E I and E O . The node set W v and edge set W e properties are defined as the same as the set W v and W e in TDG.

Using these definitions, the problem to be solved can be formulated as below. Given an application T DG = (V, E, W v , W e ) and a target MPSoC platform ACG = (P, L), find a transformation T : T DG → T DG (V , E , W v , W e ) and a scheduling function S : V → P , so that the total execution time of the parallel program is minimized. As streaming applications are executed in pipelined way, the pipeline execution time is bounded by the slowest stage. Hence, the objective is minimizing the maximum runtime (including the execution and communication time) among all pipeline stages, i.e., min(max(

i 1 t(n i 1 , p1) + i 1 n k 1 ∈succ(n i 1 ) (c(e i 1 k 1 )/R), i 2 t(n i 2 , p2) + i 2 n k 2 ∈succ(n i 2 ) (c(e i 2 k 2 )/R), ..., i |p| t(n i |P | , p |P | ) + im n k |P | ∈succ(n i |P | ) (c(e i |P | k |P | )/R))) (1)
subject to

S(n i j ) = pj (2) ∀pj ∈ P, ∀n k j , n i j ∈ V , ∀e i j k j ∈ E , S m(n i j ) ≤ mp(pj ) ( 3
)
where R is the bandwidth of each link between processing elements in P, each sum term calculates the execution time of pipeline stage i, including the execution time of the tasks scheduled on processor p i and communication time for the traffic send to the successor tasks. Condition 2 restricts that one task is only scheduled to one processing element. Condition 3 ensures that the total memory consumption of all tasks that are assigned to processing elements that should not exceed the memory space of p j .

Given the C program of a streaming application, the following steps as shown in Figure 1 will be performed under the framework of parallelizing sequential program:

Step1: Select the computing hot-spot region (CHR), which is the procedure with the maximum execution time, of the whole program through dynamic profiling. The source code in computing hot-spot region is scanned to find the relations of parameters, pointers with data structure, data dependence and control dependence between statements through a top-down method using in-house developed tools. Then build the T DG of the computing hot-spot region. The main procedure of a program is usually chosen as the initial computing hot-spot region. Step2: Transform the T DG to a directed acyclic graph T DG to eliminate the control dependences and inter-iteration data dependences (P1). Step3: Allocate and schedule each node in the T DG to a proper processing element in the ACG to form a thread (P2). A pipeline stage consists of one or more threads. Step4: Apply split/merge parallelizing technique to balance the workload among all pipeline stages (P3). Step5: Evaluate the execution time of the scheduled pipeline. The objective of scheduling is to minimize the maximum runtime among all pipeline stages. To further reduce the execution time, choose the bottleneck stage with the largest runtime to be the new computing hot-spot region, then jump to step 1. Step6: Based on the final scheduling result, generate the parallel program for the application and compile the codes to get the executable files for the hardware platform.

The sub-problems dependence transformation (P1), task scheduling (P2), and workload balancing (P3) will be defined and the solutions to these subproblems will be discussed in the next subsections.

Dependence Transformation

The sub-problem P1 is defined as: Given the TDG of an application, find a transformation method T : T DG → T DG , which removes the redundant dependences to make the output graph satisfying the following requirements: 1) No control dependence edge. 2) No dependence cycle, i.e. no SCC composed of more than one node. Thus T DG is an acyclic directed graph with only data dependences. path is selected to execute, the other branch path will not be executed. For every branch task, the branch exclusive array is updated according to the following steps.

Branch Control Dependence

First, scan the branch information i b of every branch node of which the type label l(n

) is n∈V C , or n∈V C ∩V B , or n∈V B ∩V L . Assume n j ∈ V B , i b (n j ).branch level = N . Search branch label[i] and branch condition[i] of i b (n j ) starting from i = 0.
Second, traverse the ith level branch control nodes outgoing control edges and find the other nodes controlled by it. Check the branch label and branch condition of these nodes with those of node n j to determine whether they are mutually exclusive. If so, add these nodes into the branch exclusive array of node n j .

Third, the ith level branch control node is combined with the branch node n j . Let i = i + 1, repeat steps 2 and 3 until i = N .

Finally, after finishing process all the branch nodes, delete all the branch control nodes in T DG.

Through this transformation, each branch control node is merged with its successive branch nodes in different branch paths to form new nodes.

Loop Control Dependence Due to their repeating characteristics, loops are the important structures to explore for parallelizing a program. According to the conventional methods [START_REF] Ceng | MAPS: An Integrated Framework for MPSoC Application Parallelization[END_REF][START_REF] Rangan | Decoupled Software Pipelining with the Synchronization Array[END_REF][START_REF] Sinnen | Task Scheduling for Parallel Systems[END_REF] , all tasks in a loop should be merged into one large task to eliminate the loop control dependences. Indubitably this will impede exploiting the parallelism of the tasks in the loop. In addition, the large task is likely to become the bottleneck of parallelization.

We apply the speculation technique to remove the loop control dependences. Provided that the loop runs fixed times or many iterations, the loop is regarded as biased. As a matter of fact, through profiling, many loop behaviors of streaming applications are highly predictable. The following steps are performed to resolve loop control dependences.

First, check the values of loop level, loop label, and loop num to see whether this loop is biased, i.e., if there exists loop num[i] > 0 (ith level loop is biased). The loop control dependences of this type of loop are seen as highly predictable, and they are chose to speculate.

Second, remove the selected ith level loop control edges. Insert code to detect the mis-speculation.

-Insert an unconditional branch statement, such as while (true), in all the loop tasks that are dependent on this loop control task. If loop num[i] is known by profiling, a counter is inserted in every loop task. -A mis-speculation procedure needs to be inserted in the last task of this loop, when loop num[i] is a variable. The mis-speculation detection is achieved by copying and updating the computing results of last iteration of the current iteration which is running. If the predicted loop is not taken, the misspeculation procedure is responsible for recovering data and jumping to the loop exit path.

Resolving Inter-iteration Dependence The mechanism for resolving interiteration dependences is as follows. First, on the transformed T DG, provisionally remove all the inter-iteration dependence edges with their dependence distance greater than 1. Next, merge tasks which belong to a SCC in the transformed T DG into one SCC node. Then, identify the SCC distance d s of each task. We can identify the set of inter-iteration edges that satisfy any one of the following rules. Assume an inter-iteration dependence edge e ij exists, e ij representing the dependence from task n i to task n j .

-If d dep (e ij ) > d s (n i , n j ) ≥ 1, the inter-iteration dependences can be ignored by the pipeline stages built on current SCCs. When the next stage is going to be executed, the inter-iteration dependence relation it relied on has been already satisfied. Thus, removing e ij does not affect the parallel scheduling on the current SCCs.

-If d dep (e ij ) = 1, check the d s (n i , n j ) of this edge. If d s (n i , n j ) = -1, it means
task n i and n j are in the same SCC and e ij can be ignored. If d s (n i , n j ) = 0, actually this is the special case, e ij can be removed.

The inter-iteration dependence edges that satisfy any of the above rules are regarded as redundant and are allowed to be removed. All other provisionally inter-iteration dependences are not allowed to be removed and are added back to the T DG before invoking the transformation.

Finally, grouping the tasks and data dependence edges which are involved in a dependence cycle, then the acyclic directed graph T DG has been built.

Scheduling

The scheduling sub-problem is to find a scheduling function S : V → P with the objective of minimizing the maximum run time among all pipeline stages. The constraint of the memory space at each processor need be be satisfied. A heuristic scheduling approach is used here which allocates the tasks in T DG to the processors such that the workload on each processor is kept balanced.

Two-level priorities are defined to indicate the order of scheduling. A task batch is the group of tasks which can be executed in parallel according to their dependence relations. Each task batch has a queue structure to load the tasks, and a corresponding batch priority to indicate the execution order of the task batch. Every task in the task batch has a task priority to indicate its order in the same task batch in case of resource confliction.

The batch priority (TP) of a task batch is set according to the dependence relations through a breadth first search method. The smaller its batch priority value, the earlier the task batch can be executed. Once a task is assigned to a processing element, the batch priority values of its successor tasks will be updated. Thus, the batch priority can be seen as a dynamic priority.

The task priority of a task n i in a task batch queue is defined as a linear function of three major factors. The larger its task priority value, the higher priority of the task is in the task batch. Given T DG = (V , E , W v , W e ), ACG = (P, L), ∀n i ∈ V

T P (n i ) = α × bl(n i ) + β × DM em(n i ) + γ × SM em(n i ) ( 4
)
where the major factors are defined below.

bl(n i ): the bottom level of n i is the length of the longest path starting from n i [START_REF] Sinnen | Task Scheduling for Parallel Systems[END_REF]. If the bl(n i ) is high, it implies that n i is a critical task and should be given a high priority corresponding to a larger task priority. -DM em(n i ): the consumption of the communication buffers for task n i . The DM em

(n i ) = j c(e ij ),
where edge e ij originates from n i , c(e ij ) represents the communication traffic of the dependence edge e ij . If DM em(n i ) is high, it means that the inter-processor communication traffic may be large. The task with grater value of dynamic memory should be allocated to the processor with a higher priority. -SM em(n i ): the memory requirements of instruction and static data, which are obtained at the profiling time. -The scale coefficients of α, β, and γ are used to normalize the elements. The three coefficients are defined as:

α = 1/len(cp), β = 1/(R × min(t(n i , pj ))), γ = 1/M
ALGORITHM 1: Parallel Scheduling (T DG , ACG) where len(cp) indicates the length of the critical path (cp), which is longest path in T DG , R represents the link bandwidth between processors, M represents the memory space of processing element.

Input: Graph T DG (V , E , W v , W e ),
The scheduling algorithm is divided into two parts: 1) sort the tasks in the T DG in task batches and calculate the batch priorities and task priorities, and 2) assign each task to a proper processing element to form a thread according to both batch priority and task priority. Each processing element is responsible for one thread at one time. This approach focuses on good workload balance, and also takes into account mutually exclusive branch tasks identified in earlier phases and data locality optimization. Algorithm 1 outlines the scheduling algorithm. For each task it is necessary to determine which processing element the task should be scheduled to and the time slot the task will be execute on the processing element. The principle of processor selection is according to [START_REF]Specification for the Advanced Encryption Standard[END_REF].

∀n i ∈ V , ∀pj ∈ P AvailableF actor(n i , pj ) = λ × DL(n i , pj ) -max(DRT (n i ), P EAT (n i , pj )) -t(n i , pj ) (5)
where DL(n i , p j ) represents the data locality factor which indicates the data reuse time and reduction of communication time, when assigning task n i to processing element p j , λ represents the proportion of DL adjusting factor, DRT (n i ) represents the data ready time of task n i , P EAT (n i , p j ) represents the avail-able time of the processing element p j for task n i , and t(n i , p j ) represents the execution time of task n i on the chosen processing element p j .

In selecting the suitable processing element for a task, the processing element with the maximum value of AvailableFactor is chosen. The larger this value is, the higher available level the processing element has. As described above, the principle for selecting a proper processing element also takes into account the mutual exclusive of the branch tasks. The main difference between scheduling of the branch tasks and that of the ordinary tasks is the computation of the AvailableFactor as mentioned before. The computing about P EAT (n i , p j ) and t(n i , p j ) for the branch tasks should take more attention on the mutual exclusive property.

Workload Balancing

In this step, the basic software pipeline technique is applied in conjunction with the split and merge technique to further balance the workload. Given the initial scheduling result, split and merge the stages to further reduce the execution time of the computing hot-spot region. It is achieved by assigning more threads to large stages and merging small stages into one stage so that the workload balance and efficiency to the processing elements can be improved.

If the outer-most loop of a program as a computing hot-spot region is split into several task sets in a pipeline style, then each task set is called a stage. Usually one stage can be assigned to one or more threads. If the inner loop of a program or a stage of a computing hotspot region is spilt task sets in a pipeline style, then each task set is called a sub-stage. Usually the stage consists of one or more sub-stages. If the minimum execution cost among processing elements account of less than 50% of the maximum one, or the estimated pipeline execution time does not satisfy the user's requirement, the scheduling result is regarded as workload imbalanced, then the granularity of the pipeline stages need be adjusted.

A heuristic approach is used here which operates in three steps. First, check whether the loop in the largest stage can be split into several independent iterations without further profiling. If so, split the largest stage and apply DOALL to assign the iterations to different threads. Second, merge those small stages into one processing element or insert them into the spare time of other working processing elements. Third, after the simple split/merge processing, the stage with largest runtime is selected as the new CHR. Then repeat the steps of transformation and scheduling as shown in Figure 1.

Experimental Results

To evaluate speedups of our method, the parallelization framework in Figure 1 is applied to the sequential C programs to produce the parallel codes. Table 2 lists the software tool chains that we have used through the parallelization. Our experiments are conducted on the multi-FPGA-based networks-on-chip emulation [START_REF] Nethercote | Valgrind: A Framework for Heavyweight Dynamic Binary Instrumentation[END_REF] get information about function call and computing cost of sequential code gcov [START_REF]Gcov -Using the GNU Compiler Collection (GCC)[END_REF] get the branch selecting information Dependence Analysis-developed in-house VarAnalyzer analysis the information of variables in a function, including type, size, defineuse chain, and life time. DepViz analysis the data dependences and control dependences among functions. Cross Compiler mipsisa32el-gcc compiler for MIPS32 compatible architecture compiler option -O0 -march=4rkc -nostdinc -g -fno-delayed-branch platform [START_REF] Liu | Building a Multi-FPGA-based Emulation Framework to Support Networks-on-Chip Design and Verification[END_REF]. Eight 32-bit compatible MIPS4Kc RISC cores were instantiated on the platform, which can be configured to 2, 4, and 8 cores. Each processor core is attached to the advanced microcontroller bus architecture bus in order to connect with peripheral memory, communication, and debugging interfaces. A 3x3 mesh of routers is used to interconnect RISC cores, which implemented with deterministic routing algorithm. Each router consists of five input/output ports, 16-depth first-in first-out buffers, and two virtual channels for each port.

The proposed method is applied to two realistic streaming programs: Advanced Encryption Standard (AES) [START_REF]Specification for the Advanced Encryption Standard[END_REF] and T264 decoder [15]. The characteristics of programs are shown in Table 3.

T264 decoder [15] is one of the open source video codecs based on H.264 standard. The T264 decoder carries out the complementary processes of decoding, inverse transform, and reconstruction to produce a decoded video sequence. It processes frames of video sequence in units of a macroblock and each macroblock of 16x16 pixels. The outer program loop is responsible for decoding each frame of the video sequence, while the inner loop deals with every macroblock of a frame. T264 dec parse slice function is responsible for decoding the frames. It is also the computing hot-spot region. In the T264 dec parse slice function, there is an inner-loop to do the decoding work of macroblocks in the frame. Actually there are complicated dependences among macorblocks and frames. The TDG and T DG of T264 dec parse slice function are shown in Figure 2. The arcs for inter-iteration dependences are represented with red solid lines. Control dependences are represented with dashed lines. Data dependences are represented with black solid lines. Every node is recorded the global variable, local variable, and extern parameters and their sizes. The node for stands for the loop control node, which computation and storage cost are not accounted. If we use conventional parallelization method to deal with the multi-layer loop structure which makes the inner-loop one large SCC or one large task, the parallelism in the loop cannot be exploited fully. And when the tasks are assigned to stages, the workload balance will be very poor as shown in Figure 3. The T264 dec parse slice is split into two stages, each is assigned to one thread. The first stage is the inner loop which does the decoding work for macro blocks, while the second stage extracts the prediction information from the current frame for the later frames.

Using the proposed method, the intricate dependences are resolved and the inner-loop is split into 4 sub-stages. At the same time, we partition the data transferred from stage 1 to stage 2 to make sure that each thread of stage 2 is only given the ownership of a dedicated block of data. Each stage 2 thread follows one part of the thread of stage 1 as shown in Figure 4. With one thread assigned to stage 1 and 3 threads assigned to stage 2, the speedup can be improved further. Through analyzing the inner loop we find it that the frame-decoding operation can be partitioned into blocks which consist of several macroblocks. So we gather the dependent macroblocks into one thread and schedule them on the eight-core platform. We assign 8 threads to T264 decoder program which are in a pipeline style, and each thread is mapped to a processing element.

Figure 5 shows that through the proposed method, a speedup can be achieved on the two case studies. The baseline is the conventional parallelization method which merges control dependences into one task and only explores the parallelism of the outer program loops. The 4-core and 8-core represent the parallelization result on a four cores and eight cores platform respectively. Under the proposed parallelization method, the speedup is 5.48x for T264 decocder and 5.12x for AES program respectively. As shown in the figure, the parallel scheme on four cores platform makes full use of the hardware of which the efficiency is more than 80%. Since the workload balance of the parallelization scheme is affected by the characteristic of application, the efficiency of processors on eight cores platform is smaller than on four cores platform.

Conclusions

In this article, we have proposed method for embedded system applications to exploit the coarse-grained pipeline parallelism hidden in multi-layer loops. The method first resolves intricate dependency by transforming the task dependence graph, which is then scheduled to the target MPSoC in parallel to minimize the maximal execution time of a pipeline stage. The parallel scheme is adjusted in a heuristic way to further improve performance. The experimental results of two typical applications confirm the efficiency of the method in practical. The method will be applied to parallelizing other embedded applications on multicore embedded systems. The sequential program parallelization needs to be integrated with these techniques, such as eliminating redundant dependences, task scheduling, and independent multi-threading to extracting the pipeline parallelism.

Fig. 1 .

 1 Fig. 1. Framework of parallelizing sequential code.

Fig. 2 .

 2 Fig. 2. TDG and T DG of T264 dec parse slice function.

Fig. 3 .

 3 Fig. 3. Scheduling for T264 dec parsc slice function.

Fig. 4 .Fig. 5 .

 45 Fig. 4. (a) Scheduling for T264 decoder at inner loop level, (b) Mapping result for stages to processing element. The SUB2-1 represents the first thread of stage2.

Table 1 .

 1 Properties of Node and Edge C , c(eij )denotes the data traffic amount of edge eij d dep ∀eij ∈ E I , d dep (eij ) denotes the dependence distance of edge eij

	Wv
	l	∀n ∈ V , l(n) denotes the type label of node n
	m ∀n ∈ V , ∀p ∈ P , m(n) denotes the memory requirement of task n
	s	∀n ∈ V , s(n) denotes the number of strongly-connected components that task n belongs to
	i b	∀n ∈ V B , i b (n) denotes the branch information of task n
	i l	∀n ∈ V L , i l (n) denotes the loop information of task n
	t	∀n ∈ V , ∀p ∈ P , t(n, pi) denotes the execution time of task n on the PE pi. If pi, pj ∈ P
		satisfy tp(pi) = tp(pj ), then t(n, pi) = t(n, pj )
	We
	c	∀eij ∈ E ∩ eij / ∈ E

  For branch control dependence, it is important to focus on the mutually exclusive branch tasks. The tasks in different branch paths controlled by conflicting conditions are mutually exclusive. Once a branch

	Architectural Characterization Graph (ACG) P2: Task Scheduling		Source code Profiling	Program Analysis
	Assign TDG' to processing
	elements to form threads			Computing Hotspot
					Region (CHR)
	Scheduling Result
			P3: Workload	Dependence Analyzing
			Balancing
	Workload balance?	No	Split/Merge	Task Dependence Graph(TDG)
				Workload
	Yes			balance?	Transformation for TDG
	Split the sequential code to parallel parts.			No	Branch control dependences Loop control dependences Inter-iteration dependences
	Parallel code			P1: Dependence
					Transformation
	Compile the codes with cross-compiler			TDG'
	Code Generation			Choose a new
					CHR
	Exec file for	Exec file for	. . .	Exec file for
	PE 0	PE 1			PE n-1
				Code Execution
	Real-Time Operating System
	Hardware paltform

  and target MPSoC ACG(P, L). Output: Schedule tasks to processors in ACG. Part One: Calculate the batch priority of each node in V and insert each node into the associated task batch queue through a breadth first search procedure; for each task batch bj ∈ B do Calculate the TP for each task in bi according to Equation 4; Sort tasks in a task queue in the decreasing order of TP as n π 0 , n π 1 ,..., n π |b j |-1 ; end Sort the task batches in a task batch queue in the increasing order of batch priority as b π 0 , b π 1 ,..., b π |b j |-1 ;

	Part Two:

for (k = 0; k < |B|; k + +) do for (l = 0; l < |b k |; l + +) do PESelect(task n π l in bπ k , P ); Delete the task n π from V ; for each task n i ∈ succ(n π l ) do Reassign n i into a proper task batch bj ∈ B, according to dependence relations in graph T DG ; end end end

Table 2 .

 2 Software Tool

	Program Profiling
	valgrind

Table 3 .

 3 Characteristics of Benchmark

	AES
	Characteristic	128-bit plaintext and 128-bit encrypt key
	T264 Decoder
	Version	0.14
	Sequence	forman, akiyo, container
	Input Size	QCIF, 176x144, 300 frames
	Characteristic	99 macroblocks/frame, two B frames between P frames, no rate control, deblock,
		CAVLC entropy coding

This work is supported in part by National Natural Science Foundation of China under the grants 60873112 and 61028004, National High Technology Research and Development Program of China under the grant 2009AA01Z109, and the Huawei Innovation Research Program under the grant YJCB2011033RE.


Guangchao Yao

Yao Zheng

Limin Xiao

Li Ruan
email: ruanli@buaa.edu.cn

Zhen Lin

Junjie Peng

An Effective Approach for Vocal Melody Extraction from Polyphonic Music on GPU

Melody extraction from polyphonic music is a valuable but difficult problem in music information retrieval. The extraction incurs a large computational cost that limits its application. Growing processing cores and increased bandwidth have made GPU an ideal candidate for the development of finegrained parallel algorithms. In this paper, we present a parallel approach for salience-based melody extraction from polyphonic music using CUDA. For 21 seconds of polyphonic clip, the extraction time is cut from 3 seconds to 33 milliseconds using NVIDIA GeForce GTX 480 which is up to 100 times faster. The increased performance allows the melody extraction to be carried out for real-time applications. Furthermore, the evaluation of the extraction on huge datasets is also possible. We give insight into how such significant speed gains are made and encourage the development and adoption of GPU in music information retrieval field.

Introduction

In the field of music information retrieval (MIR), the extraction of acoustic features is always the first step. After extraction, the features are then utilized by further applications. Melody, the fundamental frequency (F0) contour of the polyphonic music's lead vocal [START_REF] Poliner | Melody transcription from music audio: Approaches and Evaluation[END_REF], has been a remarkable feature owing to its numerous applications in the past few years. Having melodies, we can use them in many ways: the most attractive one may arise from the field of query by humming [START_REF] Dannenberg | A comparative evaluation of search techniques for query-by-humming using the MUSART testbed[END_REF], where the melody fragment of humming as a query will be fuzzily searched in the feature database. Apart from that, they are also frequently used in singing voice separation, singer identification and extraction of musical structure, etc.

Unfortunately, the promising applications cannot hide the difficulty of the melody extraction from polyphonic music. The complexity of the task is twofold -firstly, due to the superposition of all instruments which play simultaneously, the accuracy of the extraction is still hovering at a relatively low level. This can be seen from the results of the melody extraction of Music Information Retrieval Evaluation eXchange (MIREX) [START_REF] Downie | The music information retrieval evaluation exchange 2005-2007: A window into music information retrieval research[END_REF] in recent years. Secondly, most methods submitted to MIREX are very time-consuming [START_REF] Poliner | Melody transcription from music audio: Approaches and Evaluation[END_REF]. The inherent high computational complexity of audio signal processing is the chief culprits. At the same time, for better evaluation of the algorithm, the corpora is becoming larger and larger, which may contain as many as thousands (or even millions, in the case of some commercial databases) of audio files. All these factors hinder the development and application of the melody. The progressing rate of the extraction could be greatly improved if the execution time was reduced sharply enough to allow a quicker evaluation and tuning of algorithm parameters.

The solution to the above problem relies on the improvement of the computational capability. Some researchers [START_REF] Schmidt | Efficient Acoustic Feature Extraction for Music Information Retrieval Using Programmable Gate Arrays[END_REF] use FPGA to accelerate the extraction of the acoustic features. However, this method features a long developing cycle and is very difficult. Clusters are another common way to solve this problem. Marsyas is used to implement efficient distributed audio analysis algorithms [START_REF] Bray | Distributed audio feature extraction for music[END_REF]. But the construction of a cluster is expensive and inefficient. Besides the distributed parallel method, there are some on-chip parallel ones, i.e. multi-core and many-core methods. Owing to much more cores on many-core architecture than that on multi-core one, the former can provide a tremendous amount of computational capability than the latter. A typical many-core is GPU. The latest NVIDIA card, e.g. GTX 690 has up to 3072 cores. The availability of programming models such as CUDA has also made GPU a strong candidate for performing many computation intensive tasks.

The adoption of GPU has permeated almost all areas which require significant computational resources. Many applications ranging from general signal processing or physics simulation to computational finance or computational biology can be accelerated by GPU. Battenberg [START_REF] Battenberg | Accelerating nonnegative matrix factorization for audio source separation on multi-core and many-core architectures[END_REF] uses the multi-core and many-core to accelerate the non-negative matrix factorization for audio source separation, and the result reveals that the many-core architecture has more advantages than the multi-core's. Specific to the extraction of acoustic feature field, Schmä decke [START_REF] Schmä Decke | GPU-based acoustic feature extraction for electronic media processing[END_REF] accelerates six different exemplary features, which are among the time domain, frequency domain, and on the autocorrelated signal. However, all these features are the most basic ones of music, not including the melody feature. Although some features can be used to compute the melody, they can only work for pure voice instead of polyphonic music.

Melody extraction from polyphonic music, as mentioned above, remains a challenging and unsolved task, the overall accuracy is around 70 % which is lower compared with melody extraction from MIDI files. Since classic features cannot estimate the melody accurately, more complicated approaches are proposed. There are some designs based on the source/filter model [START_REF] Ozerov | Adaptation of bayesian models for single-channel source separation and its application to voice/music separation in popular songs[END_REF], which is sufficiently flexible to capture the variability of the singing voice and the accompaniment in terms of pitch range and timbre. Some use the salience-based method [START_REF] Guangchao | Efficient Vocal Melody Extraction from Polyphonic Music Signals[END_REF] to extract the melody. As for accuracy, the salience-based approaches give a better performance. But no matter which method is adopted, significant computational resources are necessarily needed. These arise from the nature of audio signal processing because the audio signal is usually partitioned into a large number of frames and the computation cost at every frame is high.

In this paper, we will use GPU to present a salience-based melody extraction approach to demonstrate the dramatic speedup achieved by GPU and to encourage MIR researchers to develop and reuse high performance parallel implementations of im-portant MIR procedures. Actually, the system won't be built from scratch. For example, the extraction methods using Expectation Maximization (EM) can benefit from Kumar's work [START_REF] Kumar | Fast parallel expectation maximization for gaussian mixture models on gpus using cuda[END_REF]. The methods using Support Vector Machine (SVM) will find it useful of Catanzaro's work [START_REF] Catanzaro | Fast Support Vector Machine Training and Classification on Graphics Processors[END_REF].

In section 2, the salience-based melody extraction method will be described in more details. The GPU's hardware implementation, thread hierarchy and its different kinds of memories will be stated in section 3. Section 4 will introduce our melody extraction approach on CUDA. Section 5 presents the experimentation results on different GPUs and different aspects of melody extraction. Finally we propose concluding remarks and future work to be pursued in Section 6.

Salience-Based Melody Extraction

The salience-based design has a common structure: at first, get the spectral representation of the signal. The most popular technique is the short time Fourier transform (STFT). Secondly, use the spectral representation to compute the F0 candidates.

There exist many different strategies to compute the candidates, [START_REF] Klapuri | Multiple fundamental frequency estimation by summing harmonic amplitudes[END_REF] uses the harmonic summation of the spectral peaks with assigned weights, while [START_REF] Goto | A real-time music-scene-description system: predominant-f0 estimation for detecting melody and bass lines in real-world audio signals[END_REF] lets the possible F0 to compete for harmonics based on expectation-maximization (EM) model. At last, the melody is chosen from the candidate F0 using different methods.

The procedure to compute the candidates presents significant diversities, and our approach is based on the generation of melody contours [START_REF] Guangchao | Efficient Vocal Melody Extraction from Polyphonic Music Signals[END_REF]: firstly, the pitch salience at every frame will be computed by utilizing the spectral peaks, which is called multipitch extraction. Secondly, a set of pitch contours are created using the salient candidate pitches. Then the corresponding contour characteristics will be defined, which can be used to discriminate whether the contour belongs to the melody at the melody identification stage. At the last stage -post-processing stage, vocal melody is chosen out of all contours in a three-step singing voice detection stage with the help of contour characteristics. Fig. 1 outlines these three stages through a block diagram. More details can be found in our previous work [START_REF] Guangchao | Efficient Vocal Melody Extraction from Polyphonic Music Signals[END_REF].

For better comparison with the parallel melody extraction and identifying any potential performance bottlenecks, it's necessary to perform an intensive analysis of the serial approach. We will use a polyphonic music clip with duration 21 s to verify the necessity of acceleration. Table 1 shows the execution time of the main parts of the algorithm using Intel(R) Core(TM) i5 CPU 750. From the table, it can be easily seen that the extraction of melody from polyphonic music is so time-consuming that it can hardly be used in real-time applications because the extraction will be finished in seconds. For another, the parts occupying the majority of execution time can be easily found-sinusoid extraction and pitch salience on which most attention will be paid when achieving the parallel approach. The slowness of the sinusoid extraction arises from the relative high time complexity. By contrast, the calculation of pitch salience runs slowly because of its intensive floating point operation.

GPU Programming Model

For higher throughput and better organization of so many cores, NVIDIA has done a lot of innovations on its GPU. The first to be mentioned will be its hierarchy programming model, which provides more choice at choosing different parallel granularities. Additionally, the diverse device memories are also very important in GPU programming. Each has its own special property and if used properly, they will bring an obvious promotion to the performance.

Programming Model

The NIVIDA GPU hardware is implemented as a set of streaming multiprocessors (SM), which manage a set of scalar processor cores (SP). The SPs at the same SM are equal peers and share some important memories. CUDA has a special programming model catering the hardware-Grid-Block-Thread model, i.e. the simple instruction multiple thread (SIMT) model. Users must group threads into blocks and construct a grid of some number of thread blocks. A block will be executed on a single SM. If all the threads in one block finish their job, a new block (if any) will be activated on the idle SM. One or more thread blocks to be processed by an SM are partitioned into warps with a size of normally 32 threads. These warps are then scheduled by the SM schedule unit for its execution. All the threads activated at a SM can be larger than the actual SPs for hiding the memory latency. Furthermore, threads in the same block can be synchronized for correctness.

Memory Hierarchy

In Table 2, the frequently-used memories and their properties are listed. Texture memory and local memory are rarely used under the circumstances of general purpose computation, so they are not included in the table. Actually, as long as we make good use of the memories listed in the table, the best optimization effect can be achieved.

Registers scattered on the SPs have the fastest memory access speed. But they can only be used by a single thread and cannot be used in inter-thread communication. For communication between threads, shared memory residing in the SM that is shared by the whole threads in one block can be used. Owing to its on-chip property, it also has a fast access speed. This makes it a key point for acceleration. It can be regarded as the cache, small but fast.

Besides the above, others are the off-chip device memories which often have a relative slow access speed. The global memory owns the biggest capacity which usually has several GBs, and its bandwidth is also higher than the bandwidth of CPU's memory. From compute capability 2.0 and above, the on-chip cache makes it more efficient. Furthermore, it's shared by all the blocks. Constant memory, as its name suggests, is used to save the constant variables. As it has on-chip cache, its access speed is very fast. From [START_REF] Yang | Fixing Performance Bugs: An Empirical Study of Open-Source GPGPU Programs[END_REF], the constant cache often has higher bandwidth than shared memory in spite of the off-chip property. This conclusion inspires us to use constant memory to save read-only variables. After the introduction of GPU, conclusion can be drawn that GPU is very suitable for problems which have a two-level data parallelism-one level for different blocks, and another level for different threads. If a problem doesn't show a two-level parallelism, it can also be applied on GPU as long as it has at least one-level parallelism but maybe need more tricky works. For example, the classic matrix multiplication problem must be fragmented to small square matrices to create the two-level parallelism factitiously. For another, merging two sorted arrays [START_REF] Green | GPU merge path: a GPU merging algorithm[END_REF] must use more difficult techniques to transplant the problem into GPU's platform. Fortunately, in the area of audio signal processing, the problems often show a great parallelism. For an independent music, it's usually processed by frame. The inter-frame is one level parallelism at which the weak scaling [START_REF][END_REF] can be applied, and if the problem can also be processed parallel inner-frame, this is the other level parallelism at which the strong scaling will be applicable. If a frame cannot be processed parallel, we can still get one level parallelism from different music. In a nutshell, MIR is a perfect field that GPU can show its extraordinary computing capability.

Parallel Implementation

From Table 1, the parts which most need acceleration can be easily seen. Fortunately, these parts are also the easiest parts to be parallelized. At the first stage-Multi-pitch Extraction, the main job is to transform the music data and to calculate the pitch salience. In sinusoid extraction part, CUFFT [17] is a good choice to achieve the FFT on CUDA, which is an efficient official library. In pitch salience computation part, the key operation is to calculate the salience of different frequencies at every frame. Because the salience computation of a specific frequency has no contact with the others, the computation reveals a perfectly strong scaling which can be mapped to different threads in a one-to-one mapping. Besides these two time-consuming parts, there will be two filter parts. They occupy only a small proportion of the runtime, so they can be put on CPU. But considering their stay between the sinusoid extraction and the pitch salience computation, parallelizing them on GPU will reduce the communication time between CPU and GPU. Hence we will parallelize all parts of the first stage on GPU.

In Melody Identification stage, the main job is to create a series of pitch contours and to calculate the contour characteristics. These operations have a strong data dependency. What's more, this stage occupies only a very little proportion of the whole runtime as illustrated in Table 1. So this stage will be finished on CPU.

In the last post processing stage, the operation deals with the octave error and the pitch outlier using "melody pitch mean". The calculation of melody pitch mean is finished with the help of smoothing filter which needs a lot of computation. This filter has less data dependency and it can be parallelized on GPU even the proportion of this stage is also small. Different blocks will smooth different positions, and the threads in the same block will calculate the value of the same position using reduction.

The complete heterogeneous algorithm is shown below (Note that the first stage is included for completeness). In the following section, we will focus on specific parts of the algorithm as some stages execute on CPU.

Sinusoid Extraction

Before transforming, the music data should be transferred from host memory space to device memory space, and then rearranged to fit for FFT and consequent operations, as depicted in Fig. 2. The reason for rearrangement is to make the data suitable for CUFFT functions. The rearrangement is executed after the transfer because less data will be transferred and the rearrangement can be implemented in parallel.

Afterwards, the rearranged frames will be multiplied by Hann window. Then cufftPlan1d and cufftExecR2C are used to transform the frames. The results of FFT will be complex. The real part and the imaginary part of complex are used to get the module of the transform result. At the same time, as transform is symmetric, only half of the result at one frame will be useful. In this part, many operations just transform the frame data from one form to another, so the pattern concurs with each other-transform the data one by one in every frame. Therefore, we adopt the same strategy to parallelize them: all the threads in one block will deal with an individual frame and different blocks will process different frames. This is a good example of strong scaling and weak scaling.

Algorithm

Spectral Peaks filter

The spectral peaks of transform encompass the possible pitches and are extracted for further dispose. The key of the filter is getting all the peaks first. For traditional serial algorithm, it's a very easy goal to achieve as it just traverses the array forward and returns the found peaks. When transplanting to GPU platform, the complexity appears.

Although different frames can be processed in parallel, the parallelism in a frame cannot be utilized well. The difficulty is not how to find the peaks but how to save them. If merely counting the numbers, the different threads in a block can count the peaks on shared memory and count them all using efficient reduction method. But if we need to save the peaks at the same time, the story will be different. There is a possibility that some threads in a warp find the peaks simultaneously and the conflict happens when saving them. The problem here is that counting variable will be visited simultaneously by different threads in a block. One solution relies on the use of atomic operation which serializes the visit of memory. Unfortunately, this solution has a severe drawback-reducing the access efficiency, so does the speed.

In order to balance the efficiency and the correctness, different strategies are adopted to specific length of frames. For the original transform result, the frame length is half of the FFT length plus one, i.e. 4097. The order of magnitude of the peaks will stabilize at 10 2 levels. For the sake of high efficiency, the "space for time" strategy is used. More specifically, the equal sized array with the original frame is allocated for holding the peaks, as depicted in Fig. 3. The found peaks will be put in position which has the same index with the peaks in the original frame. In this way, the threads can save the peaks without the access conflict. The reason we waste the space to save the peaks is that the peaks will be further disposed, and none of all will be saved. The consequent work is to find the max peak which is used to carry out the filter. Finding the max value of an array is a classic problem on CUDA which can be finished using the CUBLAS [18] library or achieved by ourselves using reduction. Once getting the max peak, we can perform the filter of peaks. The problems here are the same as finding peaks, namely how to save the remaining peaks. But unlike the previous finding peaks, the order of magnitude of the peaks has reduced to 10 1 levels from experiment, so this means the remaining peaks in the array is rather sparse (length of thousands of array has dozens of peaks). From the above result, the probability of conflict when the threads in a warp save the peaks simultaneously will be low. In this case, we can use atomic operation to finish the filter.

Although the atomic operation assures the correctness of counting, the result of atomic operation cannot be used directly as it may be added by another thread before saving the peaks. The solution is using the return value of the atomic function of CUDA. Because it returns the old value and incrementing the old value will exactly be the right position for the remaining peak. In addition, owing to the fast access speed and shared by the whole threads in a block, the shared memory is a good place to put the counting variable. The part of salience peaks filter is the same as the spectral peaks filter, so it will not be described in detail.

Pitch Salience

After filtering the spectral peaks, they will be used to calculate the pitch salience. The calculation reveals a perfect parallelization as the calculation of each bin has no data dependency with other bins' calculation. So it can be easily transplanted to GPU platform, the strategy is different blocks process different frames, and threads in a block dispose a frame's calculation of 480 different pitch saliences. On one hand, due to the frequent access of spectral peaks and their small size, putting them in shared memory is a good idea to accelerate the access speed. On the other hand, because the calculation of pitch salience has so many floating point operations, the requirement for registers are huge and this will limit the threads number activated at the same time. Our solution is to use the constant memory to hold the constant variables needed in calculation. This is a good approach to optimize. The other used optimizations comprise the extraction of common sub-expressions and loop unrolling, etc. Although no advanced techniques are adopted, we can still achieve a high speed-up from the following evaluation.

Allocated frames

Original frames

Evaluation

Various experiments are presented to clarify the effectiveness of GPU in our approach. At first, the overall performance of the system and the parallelized parts on each hardware platform are presented. After this, the performance along with the music length is evaluated. Furthermore, the efficiency of our peaks filter is demonstrated. At last, the influence to the ultimate accuracy is tested.

We have implemented our algorithm on two different commodity GPUs, whose hardware specifications are listed in Table 3. The CPU version runs on an Intel(R) Core(TM) i5 CPU 750 platform.

The corpus used to evaluate the efficiency of the system is from MIR-1k which is a publicly available dataset proposed in [START_REF] Hsu | On the improvement of singing voice separation for monaural recordings using the MIR-1K dataset[END_REF]. It has 1000 song clips with a duration ranging from 3 to 12 seconds, and the total length is up to 133 minutes.

Overall Performance

The same polyphonic music clip as in performance analysis of serial approach is used to measure the overall performance of the system and the parallelized parts. The runtime of the parallelized parts and the whole system on different platforms is demonstrated in Fig. 4. From Fig. 4, we can see all accelerated parts achieve a positive speed up, especially two most time-consuming parts. In the pitch salience computation part, the runtime is reduced from seconds to lower than 10ms. This is a huge acceleration. The high acceleration can attribute to the perfect parallelism and the use of shared memory. Inspired by these two parts, the whole system also gets a considerable accelerationfor a length of 21 s audio clip, the extraction time is reduced to less than 40 ms on GTX 480. The speedup is nearly to 100 times. This implementation is sufficient for real-time applications, such as query by humming. What's more, our system can be applied to massive datasets which are often used to verify the effectiveness of algorithms. The approach can reduce the time to tune the variables or to develop new methods dramatically. So researchers can pay more attention to the algorithm itself rather than the performance.

Except the two parts mentioned above, other parts get a relative small acceleration. Even though they occupy a small proportion of execution time, it's meaningful to find out a reasonable explanation. In post processing stage, a 2-times acceleration is achieved approximately on GTX 285 platform. The final cause is shown in two aspects: at first, the proportion which can be accelerated is only 70 %. This means the speed-up upper bound is only 3.3 times according to the Amdahl's law. Secondly, the algorithm complexity of smoother filter is only O(n), so the ratio of operations to elements transferred is O [START_REF] Poliner | Melody transcription from music audio: Approaches and Evaluation[END_REF]. Performance benefits can only be more readily achieved when this ratio is higher. Furthermore, although a high parallelization exists in the smooth filter, the time is wasted on the space allocation on GPU and transfer between the host and the device. So we can reach such a conclusion that a high proportion of parallelizable part and computation complexity are the prerequisites for high speed-up. The reason pitch salience computation can also achieve a very high speed-up even its computation complexity is O(n) is that it has no transfer between the host and the device. In addition, the previous hidden coefficient is large. The peaks filter part will be better explained in detail later.

Influence of music length

As the melody extraction is processed mostly frame by frame, the extraction time will increase gradually with the growing of audio length, as illustrated in Fig. 5. Although the entire trend of the line growth is incremental, there still exist some points which descend with the growing of audio length. This is because of the influence of specific audios. That means the extraction procedure has different levels in difficulty. It is obvious that polyphonic music with lower energy of background will be much easier to extract the melody than the music with a strong energy of background intuitively as the former has less candidate pitches, and this will reduce the computation time of pitch salience.

Efficiency of parallelized finding peaks

From Fig. 4, we can also see that the two peak filter parts have a relative low speedup. This phenomenon arises from the massive space allocation on device memory under the situation of no transfer between the device and the host, and the low algorithm complexity. For a problem with a low computation complexity, the execution time will be dominated by the space allocation on GPU and transfer between the host and the device. Under such a circumstance, reusing GPU device memory will be important and it's also possible to reuse previous allocated space to put the found peaks. For example, the space for rearrangement of original audio is large enough to hold the peaks. In this way, the time of release and allocation space can be reduced once.

The peaks filter can be divided into more fine-grained steps: firstly, find all the peaks, and then the max peak. Thirdly, execute the first level peaks filter. At last, do the second level peaks filter. If just ignoring the space allocation to compare the computation time, we can still see the positive effect of our parallelization, as illustrated in Table 4. Owing to the "space for time" strategy, the operation of finding peaks achieves an obvious speed-up. The classic problem of finding max doesn't demonstrate a speedup, and the too few data are to blame. The filter using atomic operation shows a small speed-up. It is understandable that few data and access conflict both exist in the extraction. The second filter works better than the first one as the possibility of conflict is smaller than that in the first filter. From Table 4, we can also see that the allocation and release time occupy a big proportion if adding up all the time and comparing it with Fig. 4.

Influence to the accuracy

The operation of double-precision floating point on GPU is time-consuming, so the single-precision floating point is adopted in our system. Although the precision is reduced, the accuracy doesn't drop. The overall accuracy on our parallel system is 73.7 %, the same as the serial approach on dataset MIR-1k. But the extraction time is reduced from nearly one hour to less than one minute. This promotes the efficiency of development tremendously. So we can spend more time on the improvement of the accuracy of the algorithm itself, and do not need to concern about the performance.

Conclusions and Future Work

Melody extraction from polyphonic music is a valuable problem because the melody can be used in many valuable applications. However, the relative long extraction time and the low accuracy limit its extensions. The extraction can be accelerated by GPU due to its high computation capability. In this paper a fast extraction approach based on GPU is presented. The results show that GPUs are well suited for audio signal processing problems as its characteristic is that the frame is processed one by one. Our parallel implementation reduces computation times by nearly two orders of magnitude on GTX 480. The acceleration is so tremendous that our system can be applied to some real-time applications, such as query by humming. Moreover, another benefit from acceleration is that it can greatly reduce the development time, so we can take more time to improve the accuracy of melody extraction.

In the future, we will implement a real query by humming application using our parallel extraction approach. At the same time, the tune of parameters and new methods about melody extraction will be verified on GPU platform deeply.

Fig. 1 .

 1 Fig. 1. Block diagram of our approach, it includes stages: Multi-pitch Extraction, Melody Identification and Post Processing. The output of first stage is salience peaks which are used to construct the contours in second stage. At the last stage, the melody is generated.

Fig. 3 .

 3 Fig. 3. Space for time strategy for finding the peaks (All elements in allocated frames are zero except the peaks)

Fig. 4 .

 4 Fig. 4. Performance comparison of different parts and the whole system (Time for CPU and GTX 480 is marked for clarity)

Fig. 5 .

 5 Fig. 5. Extraction time along with the growing of audio length on different platforms

Table 1 .

 1 Execution time of every part (ms).

		Sinusoid	Spectral	Pitch sali-	Salience	Melody	Post	Total
		extraction	peaks filter	ence	peaks filter	identification	processing	time
	Running time	307.8	9.7	2611.8		4.1	0.2	7.7	2941.9
	Proportion 10.5%	0.3%	88.8%		0.1%	0	0.3%	1.0

Table 2 .

 2 Properties of the frequently used memories.

	Name	Action Scope Speed	Cached
	Register	Thread	Fastest	N
	Shared Memory	Block	Fast	N
	Global Memory	Grid	Slow	Compute Capability 2.0 and above
	Constant Memory	Grid	Slow	Y

  Hybrid Melody Extraction from Polyphonic Music

	Input: Polyphonic Music pm
	Output: Melody m
	1: Stage I: Pre-Processing CPU ::
	2:	CPU :: read the polyphonic music to CPU memory
	3:	CPU → GPU :: transfer the music to GPU memory
	4: Stage II: Multi-pitch Extraction GPU ::
	5:	GPU :: rearrange the music data
	6:	for every frame f parallel do
	7:	GPU :: sinusoid extraction
	8:	GPU :: spectral peaks filter
	9:	GPU :: pitch salience computation
	10:	GPU :: salience peaks filter
	11:	endfor
	12: Stage III: Melody identification CPU ::
	13:	GPU → CPU :: read salience peaks to CPU memory
	14:	CPU :: pitch contour
	15:	for every contour c do
	16:	CPU :: contour characteristics
	17:	CPU :: singing voice detection
	18:	endfor
	19: Stage IV: Post Processing mainly on CPU ::
	20:	for i=1:3 do
	21:	CPU → GPU → CPU :: melody pitch mean
	22:	CPU :: remove octave error
	23:	CPU :: remove pitch outlier
	24:	endfor
	25:	CPU: return the melody

Fig. 2. Rearrangement of the original audio (For music with 44.1 kHz sample rate, frame length is 2048, hop size is 441, FFT length is 8192, zero padding is 6144)

Table 3 .

 3 GPU Specifications.

	Properties	GTX 285	GTX 480
	Number of cores	240	480
	SMs	30	15
	Cores per SM	8	32
	Global Memory(GB)	1	1.5
	Memory Bandwidth (GB/s)	159	177
	Shared Memory (KB)	16	48
	Register (KB)	16	32
	GPU core clock rate (MHz)	648	700
	Compute Capacity	1.3	2.0

Table 4 .

 4 Execution Time of EachStep at the Spectral Peaks Filter (ms).

	Hardware	Find peaks	Max peak	First filter	Second filter
	CPU	8.54	0.29	0.38	0.44
	GTX 285	1.26	0.32	0.29	0.07
	GTX 480	0.24	0.38	0.11	0.03

Acknowledgments

. This research was funded by the Hi-tech Research and Development Program of China (863 Program) under Grant No.2011AA01A205, the National Natural Science Foundation of China under Grant No.61232009, the Doctoral Fund of Ministry of Education of China under Grant No.20101102110018, Beijing Natural Science Foundation under Grant No.4122042, the fund of the State Key Laboratory of Software Development Environment under Grant No.SKLSDE-2012ZX-07 and Shanghai Science and Technology Innovation Action Plan under Grant No.11511500400.