Reducing Memory Usage by the Lifting-based Discrete
Wavelet Transform with a Unified Buffer on a GPU
Takuya Ikuzawaa,b , Fumihiko Inoa,∗, Kenichi HagiharaaaGraduate School of Information Science and Technology, Osaka University,
1-5 Yamada-oka, Suita, Osaka 565-0871, Japanb
Industrial ICT Solutions Company, Toshiba Corporation,
72-34 Horikawa-cho, Saiwai-ku, Kawasaki, Kanagawa 212-0013, Japan

Abstract
In this study, to improve the speed of the lifting-based discrete wavelet transform (DWT) for large-scale data, we propose a parallel method that achieves low
memory usage and highly efficient memory access on a graphics processing unit
(GPU). The proposed method reduces the memory usage by unifying the input
buffer and output buffer but at the cost of a working memory region that is smaller
than the data size n. The method partitions the input data into small chunks, which
are then rearranged into groups so different groups of chunks can be processed inparallel. This data rearrangement scheme classifies chunks in terms of data dependency but it also facilitates transformation via simultaneous access to contiguous
memory regions, which can be handled efficiently by the GPU. In addition, thisdata rearrangement is interpreted as a product of circular permutations such thata sequence of seeds, which is an order of magnitude shorter than input data, allows the GPU threads to compute the complicated memory indexes needed for∗

Corresponding author. Tel.: +81 6 6879 4353; fax: +81 6 6879 4354.
Email address: ino@ist.osaka-u.ac.jp (Fumihiko Ino)

Preprint submitted to Journal of Parallel and Distributed Computing

March 24, 2016

parallel rearrangement. Because the DWT is usually part of a processing pipelinein an application, we believe that the proposed method is useful for retaining theamount of memory for use by other pipeline stages.
Keywords: Discrete wavelet transform, lifting scheme, memory-saving
computation, in-place algorithm, GPU
1. Introduction
The discrete wavelet transform (DWT) [1] is a well-known method for timefrequency analysis, which analyzes a signal in both the time and frequency domains using various time-frequency representations. The DWT employs wavelets
as these representations, or scaled and translated versions of a finite-length or fastdecaying oscillating waveform called a mother wavelet. The DWT is useful forclassifying the input signal into high-frequency components and low-frequency
components in the same manner as the traditional fast Fourier transform (FFT).
Given a signal with n samples (we assume that n is a power of 2), the DWT canbe processed within O(n) time, whereas the FFT requires O(n log n) time; consequently, the DWT is more computationally efficient than the FFT. Similar to
the FFT, the DWT is a separable filter. For example, a two-dimensional (2-D)
DWT filter can be realized by row-wise and column-wise (horizontal and vertical) 1-D DWT filters. The DWT is useful in various fields, such as data compression for images [2, 3, 4], artifact removal [5], data decoding for earth observation
satellites [6] and clustering analysis for geometrical data [7]. These typical applications usually employ multiresolution analysis [1], where the DWT is applied2

recursively to low-frequency components at each resolution level of the hierarchy. Therefore, a fast DWT implementation is required to process DWT-based
applications in real-time.
Acceleration methods for the DWT can be classified into two types: algorithmbased and accelerator-based approaches. An example of the former is the lifting
scheme [8], which halves the computational requirements of the traditional filterbank scheme [9]. Furthermore, the lifting scheme is an in-place algorithm, whichis capable of overwriting the algorithm’s output onto the input buffer; thus, thelifting scheme reduces the amount of memory usage by half compared with thefilter bank scheme. In this study, an in-place algorithm is defined as an algorithm
that requires O(1) memory space, excluding the input buffer of size n.
The DWT has also been implemented in many accelerators, such as graphics
processing units (GPUs) [10], field programmable gate arrays [11], and cell broadband engines [12, 13]. The GPU is one of the most popular commodity devices
among these hardware types, where it provides a much higher memory bandwidththan the CPU. For example, the latest Maxwell architecture [14] provides 336
GB/s of memory bandwidth, which is one magnitude higher than that of the double data rate fourth generation (DDR4) memory. Consequently, many memorybound applications have been successfully accelerated on GPUs [15, 16, 17]. The
DWT is a memory-bound operation rather than a compute-bound operation, so we
consider that the GPU is a promising accelerator for the DWT.
Laan et al. [10] accelerated the lifting scheme on a GPU that is compatible
with the compute unified device architecture (CUDA) [18]. They implemented3

row-wise and column-wise 1-D DWTs as kernel functions, which are then offloaded from the CPU to the GPU. Their kernels facilitate efficient memory access by using a data layout that stores high-frequency components in a contiguous
memory region. After the high-frequency components, low-frequency components are also stored in a contiguous memory region. This data layout is useful for maximizing the effective memory bandwidth on a GPU because memorytransactions are performed per warp, i.e., a series of GPU threads that executes thesame instruction every clock cycle [18], where the threads in a warp are allowed
to access a contiguous memory region simultaneously via a 128-byte memorytransaction.
However, this layout distinguishes the output buffer from the input buffer,
thereby eliminating the advantage of in-place transformation. For signal data
larger than half of the GPU memory capacity, the data must be swapped from theGPU memory to the CPU memory. This data transfer between the CPU and GPU
degrades the overall transformation performance because the peak bandwidth of
the PCI Express (PCIe) bus is 15.8 GB/s, which is one magnitude lower than thatof the GPU memory. In addition, the DWT is usually part of a processing pipelinein an application, so it is better to retain the amount of memory for use by otherpipeline stages.
In this study, in order to obtain a fast lifting-based DWT for large-scale data,
we propose a parallel method that reduces the amount of memory usage andthat facilitates highly efficient memory access on a CUDA-compatible GPU. The
proposed method achieves memory-saving transformation by unifying the input
4buffer and the output buffer but with the cost of a tiny working memory space,which is less than the input size n. Naive unification avoids parallelizing thelifting scheme because the unified buffer causes many flow dependencies duringtransformation. To address this issue, the proposed method partitions the input
data into small chunks, which are then rearranged in an appropriate manner so theparallel DWT can be performed over the unified buffer. Thus, this data rearrangement can be regarded as a preprocessing phase that classifies the input data into
groups so different groups are independent of each other, where data-dependent
tasks within the same group are processed sequentially, but different groups canbe processed in parallel. Furthermore, chunk operations can be processed in asingle-instruction, multiple-data (SIMD) manner because each chunk comprises
multiple data elements.
In addition to the buffer unification mentioned above, the proposed method
performs the DWT via simultaneous access to contiguous regions. Thus, after the
chunks have been generated using this GPU-friendly access pattern, the generated
chunks can also be placed in appropriate positions by this favorable access pattern.
Allowing the GPU threads to compute the complicated memory addresses where
their responsible chunks must be placed is a challenging issue. In particular, irregular addresses must be computed from very small amounts of information with adata size less than n to maintain the advantage of the in-place transformation. To
tackle this issue, the proposed method represents a sequence of memory addresses
as a circular permutation. This interpretation allows us to precompute a small sequence of numbers, which can then be used as seeds to produce the sequence of5

memory addresses for use in data rearrangement.
The remainder of this paper is organized as follows. Section 2 introduces
related research on GPU-accelerated DWT implementations. Section 3 summarizes the lifting scheme, which forms the basis of the proposed method. Section
4 describes the proposed method and Section 5 presents the experimental results.
Finally, Section 6 gives the conclusions.

2. Related Work
Table 1 compares related research into the acceleration of the DWT on GPUs.
Each study employed different hardware so the peak performance is presented
with the transformation throughput to allow standardization. To the best of ourknowledge, no previous GPU-based implementation has achieved a compact data
representation and data arrangement that enables GPU-friendly memory accesses
for the lifting-based DWT.
Wong et al. [19] implemented the filter bank scheme on a GeForce 7800 GTX
GPU, where their five-level DWT implementation was based on OpenGL [20], ora graphics application programming interface (API). They selected the filter bankscheme rather than the lifting scheme because the data dependencies in the lifting
scheme are more complex than those in the filter bank scheme. These complex
data dependencies cannot be handled efficiently by OpenGL, which is specialized
for graphics applications, because OpenGL has many restrictions on its programming capability compared with CUDA. For example, OpenGL prohibits the explicit management of an on-chip cache called shared memory [18]. Despite these6restrictions, Tenllado et al. [21] implemented the lifting scheme using OpenGL.
However, they found that the filter bank scheme was faster than the lifting schemedue to the complex data dependencies mentioned above.
Franco et al. [22] implemented the filter bank scheme on a Tesla C870 GPU.
They adopted a row-column approach, which reuses a horizontal DWT kernel asa vertical DWT kernel by transposition. Their implementation transformed an
8192 × 8192-pixel image within 71 ms and the transformation throughput reached
0.9 Gpixel/s. However, the overall performance was limited by data transfer between the CPU and GPU. Transferring the data between the CPU and GPU required 163 ms, so the transformation throughput dropped to 0.3 Gpixel/s after
considering these overheads. A similar row-column approach was presented by
Galiano et al. [23].
To the best of our knowledge, Matela et al. [3] first implemented the filterbank scheme on a CUDA-compatible GPU. Their implementation was based ona block-based approach, which allows the horizontal and vertical DWT kernels
to be merged into a single kernel. Their implementation copies 2-D blocks into
the shared memory, which can be commonly accessed from threads in the same
thread block [18]. This on-chip memory is useful for saving the GPU memory
bandwidth, so that the ratio R of the memory throughput over the peak memory
bandwidth was more than 1.0. That is, on-chip data was efficiently shared among
threads. In their implementation, they distinguished the output region from the
input region.
Laan et al. [10] proposed a hybrid of the row-column and block-based ap7

proaches. Their horizontal process is same as that in the row-column approach.
By contrast, their vertical process avoids transposition by applying the vertical
DWT to vertically-long blocks. Furthermore, a sliding window mechanism wasdeveloped to reuse on-cache data fetched from the GPU memory. Similar to theblock-based approach [3], this hybrid approach increased the ratio R to more than
1.0. Using a GeForce 8800 GTX card, their implementation achieved a transformation throughput of 1.5 Gpixel/s for a 4096 × 4096-pixel image. This throughput was roughly ×3 higher than those achieved by a row-column scheme and an
OpenGL-based implementation [21]. As presented theoretically in [8], Laan etal. demonstrated that the lifting scheme was two times faster than the filter bankscheme on the CUDA-compatible GPU. Kucis et al. [24] presented an extensionof Laan’s method that reduces synchronization by allowing redundant computation between threads. This extended method achieved a transformation throughput
of 8.3 Gpixel/s on a GeForce GTX 580 card, which was 30% faster than the basic
method. However, after considering the data transfer between the CPU and GPU,the transformation throughput decreased to 0.4 Gpixel/s.
Sharma et al. [25] used OpenCL [26] to implement the lifting scheme on a
GeForce GTX 285 GPU. OpenCL is a standard for heterogeneous computing platforms, so their implementation ran on CUDA-incompatible GPUs and multicore
CPUs. Their OpenCL-based implementation achieved a transformation throughput of 0.4 Gpixel/s for a 4096 × 4096-pixel image, but it was slower than that of
their CUDA-based implementation, which was 0.7 Gpixel/s.

8

Table 1: Comparison of related research. P denotes the peak memory bandwidth for the GPU
employed. ρ and ρ+ denote the transformation throughputs without and with data transfer time
between the CPU and∑
GPU, respectively. For the L-level DWT, the number N of transformed pixL
els is given by N = k=1 N 2 /4k−1 because the level k transforms an N/2k−1 × N/2k−1 -pixel
region. R = N A/T /P is the ratio of the memory throughput over the peak memory bandwidth
P , where A and T are the amount of memory access per pixel and the execution time, respectively.
The proposed method assumes that a chunk comprises b = 512 elements.
DWT scheme
API
GPU
MemoryData size Memory
PR
Throughput (Gpixel/s)(level)capacity (MB) (pixel)
usage(GB/s)
ρρ+Wong [19]Filter bank (5) OpenGL 7800 GTX
512
2K×2K
2n
38
—
—
0.003
Tenllado [21] Lifting (5) OpenGL 7800 GTX
512
2K×2K
2n
38 0.33–0.73 0.3–0.4—
Tenllado [21] Filter bank (5) OpenGL 7800 GTX
512
2K×2K
2n
38 0.57–0.66 0.3–0.6—
Franco [22] Filter bank (1) CUDA Tesla C870
1536
8K×8K
2n
77
0.49
0.9
0.3
Laan [10]Lifting (3)CUDA 8800 GTX
768
4K×4K
2n
86
1.22
1.5
—
Matela [3]Lifting (1)CUDA
GTX 295
896
2K×2K
2n
112
1.67
2.6
—
Galiano [23]Lifting (1)CUDA
GTX 280
1024
4K×4K
2n
142
0.97
2.5
—
Sharma [25]Lifting (3) OpenCL GTX 285
1024
4K×4K
2n
159
0.11
0.4
—
Sharma [25]Lifting (3)CUDA
GTX 285
1024
4K×4K
2n
159
0.18
0.7
—
Kucis [24]Lifting (1)CUDA
GTX 580
1536
4K×4K
2n
192
3.11
8.3
0.4
This paper
Lifting (1)CUDA
GTX 5801536
16K×16K n + ⌈n/2b⌉ 192
1.27
2.8
0.7
This paper
Lifting (1)CUDA GTX 980 Ti
6192
32K×32K n + ⌈n/2b⌉ 336
1.38
5.3
2.0
Study

3. Lifting based DWTLet x0 , x1 , . . . , xn−1 be the input signal with n samples. We assume that n
is a power of two. Given the input signal, the DWT outputs its low-frequency
components c0 , c1 , . . . , cn/2−1 and high-frequency components d0 , d1 , . . . , dn/2−1 .
In the lifting scheme [8], the t-th element ct of low-frequency components and dt
of high-frequency components, where 0 ≤ t < n/2, are given by

ct = x2t +

k∑
0 −1

ui dt ,

(1)

i=−k0
k1
∑

dt = x2t+1 −pj x2(t+j) ,(2)

j=−k1 +1

where ui (−k0 ≤ i < k0 ) and pj (−k1 < j ≤ k1 ) are coefficients determined by
the mother wavelet.

9The data dependencies in Eqs. (1) and (2) indicate that in-place transformation can be achieved by overwriting dt to x2t+1 , then ct to x2t . This access pattern implies a cyclic layout, which stores the low-frequency element ct and highfrequency element dt alternately: c0 , d0 , c1 , d1 , . . . , cn/2−1 , dn/2−1 . Using this layout, the DWT can be implemented with two kernels, a high-pass filter kernel forcomputing high-frequency components and a low-pass filter kernel for computing
low-frequency components.
One drawback of the cyclic layout is that it requires data packing to obtain thefiltered result, where the high- and low-frequency components are stored in odd
and even elements, respectively (note that the indexes start from zero). Furthermore, the cyclic layout cannot maximize the memory throughput during transformation because both kernel functions cause a memory stride of two, where
the high- and low-pass kernels output odd and even elements, respectively. As afurther issue, the memory stride doubles at every hierarchical level of the multiresolution analysis, where the DWT is applied recursively to low-frequency components.
3.1. Laan’s Previous Method
Laan et al. [10] presented a solution to the data layout issue. Their method increases the efficiency of memory access by using a block layout, which stores lowfrequency components and high-frequency components in contiguous regions:c0 , c1 , . . . , cn/2−1 , d0 , d1 , . . . , dn/2−1 . However, this method separates the output
buffer from the input buffer due to the data dependencies inherent in Eqs. (1)10and (2). Consequently, it eliminates the advantage of in-place transformation andconsumes double the memory space, as mentioned in Section 1.
Consider here a block layout that naively unifies the input buffer and the output
buffer. As shown in Fig. 1(a), the high-pass filter kernel outputs high-frequency
element dt to xt+n/2 , where 0 ≤ t < n/2. However, this execution eliminates partof the original elements, xn/2 , xn/2+2 , . . . , xn−2 , which are needed by the second
kernel to compute cn/4 , cn/4+1 , . . . , cn/2−1 , as indicated in Eq. (1).
Two alternative approaches can be employed to avoid this elimination: a serialization approach that sequentially processes data-dependent tasks or a duplication approach that copies all of the elements, x0 , x1 , . . . , xn−1 , to another memory
region before outputting the computational results. The former approach fails totake advantage of massively parallel computing on the GPU. The latter approachcorresponds to the previous method [10], which separates the output buffer fromthe input buffer.

4. Proposed Method
Similar to Laan’s previous method [10], the proposed method employs a blocklayout to exploit the high peak memory bandwidth of the GPU. Furthermore, theproposed method rearranges the data elements so the odd elements follow the evenelements, as shown in Fig. 1(b). This rearrangement is performed on a per-chunkbasis, where each chunk contains b elements. Next, the method computes Eqs. (1)and (2) to complete the lifting operation [10]. For simplicity, we assume 1-D datain the following discussion.
11

d0

d2

d1

d n / 2− 2 d n / 2−1

Write

L , xn / 2− 2 , xn / 2−1 xn / 2 , xn / 2+1 , xn / 2+ 2 , L

x0 , x1 , x2 , L
Read

c0

cn / 4Lc1

cn / 4+1, xn − 2 , xn −1
L cn / 2−1(a)d 0 d1 d 2

d n / 2− 2 d n / 2−1

Writex0 , x2 , x4 , L
c0 c1 c2

L , xn − 4 , xn − 2 x1 , x3 , x5 , L

Read

L

L , xn − 2 , xn −1

cn / 2− 2 cn / 2−1(b)Figure 1: Memory access pattern in lifting-based DWT. (a) When writing high-frequency components, a naive unification scheme eliminates the original elements needed for the latter half of thelow-frequency components. (b) By contrast, the proposed scheme maintains the original elements
after writing the high-frequency components.

4.1. Data Rearrangement
The proposed method rearranges the input data according to the following two
steps.
1. Chunk generation: Given an input signal x0 , x1 , . . . , xn−1 with n samples,the method constructs an even chunk and an odd chunk within every segment of 2b elements (see Fig. 2). An even chunk contains even elements
whereas an odd chunk contains odd elements. In the following discussion,
let m (= ⌈n/2b⌉) be the number of even chunks. For simplicity, the number of odd chunks is assumed to be the same as that of even chunks; thus,
the total number of chunks is given by 2m. To facilitate highly efficient
memory access, the chunk size is assumed to be a multiple of 128 bytes be-

12

2b

2b

2b
x0 , x1 , x2 , x3 , x4 , x5 , L

L

L , xn − 4 , xn −3 , xn − 2 , xn −1

x0 , x2 , x4 , L x1 , x3 , x5 , L x64 , x66 , L x65 , x67 , L

L

L , xn − 4 , xn − 2 L , xn −3 , xn −1

C0

C2

C1

C3

C2 m − 2

C2 m −1

Figure 2: Chunk generation. Even and odd chunks are generated within every segment of 2belements. After generation, even and odd chunks appear alternately.

Before

C0

C1

C2

C3

C4

C5

C6

C7

After

C0

C2

C4

C6

C1

C3

C5

C7

Figure 3: Chunk rearrangement (m = 4). Odd chunks follow even chunks after rearrangement.
Note that C0 is an even chunk because the indexes start from zero.

cause the current CUDA accesses the global memory via 128-byte memorytransactions [18].
2. Chunk rearrangement: 2m chunks are rearranged such that m odd-chunks
appear after m even-chunks, as shown in Fig. 3. More formally, the 2i-th
chunk and the (2i + 1)-th chunk, where 0 ≤ i < m, are placed in the i-th
chunk and the (i + m)-th chunk after rearrangement.
Algorithm 1 shows the proposed chunk generation algorithm, which describes
the kernel function that needs to be executed by massive GPU threads in parallel.
This algorithm assumes that a thread block [18] containing 2b threads is responsible for generating two chunks within a segment of 2b elements. Each thread loads
its responsible element from the global memory to a shared array s of size 3 × b

13

(line 3), before performing local synchronization within the thread block. Notethat the last column of this array is a dummy region padded for avoiding bankconflicts [18] when writing the loaded element back to the appropriate location in
the global memory (line 10). This row-wise algorithm can easily be extended toa column-wise algorithm without using the shared memory. As described in [10],column-wise operations can be achieved by assigning a row of 32 neighboringcolumns to a warp, where every warp is allowed to access a contiguous memory
region.
Algorithm 1 Chunk generation kernel.
Input: Input signal x0 , x1 , . . . , xn−1 with n samples and b elements in a chunk
Output: Per-chunk rearranged elements x0 , x2 , . . . , xn−1 containing alternate
odd and even chunks
1: i ← thread index;
2: j ← thread block index;
3: si%2,⌊i/2⌋ ← x2bj+i ;▷ Load responsible element into a shared array of size

3×b

▷ Local synchronization▷ Even chunk generation

4: synchthreads();
5: if i < b then6:k ← 2i;7: else8:k ← 2(i − b) + 1;
9: end if
10: x2bj+i ← sk%2,⌊k/2⌋ ;▷ Odd chunk generation

Chunk rearrangement can be represented as a product of l circular permutations: σ0 ◦ σ1 ◦ . . . ◦ σl−1 , where σj (0 ≤ j < l) represents the j-th circularpermutation. Let (i j) be a transposition representing an exchange between thei-th chunk and the j-th junk. The chunk rearrangement for m = 4 shown in Fig. 3

14can then be formulated as a product of two circular permutations, (1 4 2) ◦ (3 5 6).
According to this interpretation, a circular permutation (1 4 2) represents a datadependent task that circulates chunks C1 , C4 and C2 . Thus, a task is a series ofchunk movements that must be processed sequentially. However, different circularpermutations are independent of each other, so they can be processed in parallel.
Therefore, the proposed method assigns a circular permutation to a thread block,which then sequentially processes a series of chunk movements that corresponds
to the assigned circular permutation. Note that a chunk contains multiple data
elements; thus, the data elements within a chunk can be processed simultaneously
by b threads from the same thread block in a SIMT manner.
A critical issue is how to formulate circular permutations according to thenumber n of samples. To address this issue, we introduce the following lemma.
Lemma 1. Let s be the size of an arbitrary circular permutation σj (0 ≤ j < l),which represents a series of chunk movements in step 2. Then, s ≤ ⌈log 2m⌉.
Proof . Let (a0 a1 a2 . . . as−1 ) be an arbitrary circular permutation, which represents a series of chunk movements in step 2, where this permutation indicatesthat chunks Ca0 , Ca1 , . . . , Cas−1 must be circulated. Then, for all 0 ≤ k < s,
0 ≤ ak < 2m because there are 2m chunks in total. According to step 2, we have
Eq. (3).
{
ak+1 /2,
if ak+1 is even,ak =(3)⌊ak+1 /2⌋ + m, otherwise.
Note that for all 0 ≤ k < s, ak can be represented as a binary number with⌈log 2m⌉ bits. In this binary representation, applying a circular right shift to ak+1
gives ak , as indicated by Eq. (3). For example, the permutation (1 4 2) has threeelements, 1, 4, and 2, which correspond to binary representations of 001, 100,
and 010, respectively, and the last two elements can be obtained by applying acircular right shift to the first element 001 (see Fig. 4). Therefore, for an arbitrary
ak , applying a circular right shift at most ⌈log 2m⌉ times to ak produces itself;thus, s ≤ ⌈log 2m⌉.
□
15

0

(1 4

2)

(001 100

010 )

0

1

Figure 4: Interpretation of chunk rearrangement with a circular permutation. A circular permutation of length M can be mapped to a necklace with M beads and two colors, where rotations are
considered equivalent. In this example, M = 3.

According to Eq. (3), we have Eq. (4):ak+1 = 2ak mod (2m − 1),(4)

where 0 ≤ k < s. Thus, the representative element a0 of a circular permutation
produces all the elements of that circular permutation. This means that the final
memory location for an arbitrary chunk can be computed from a0 , i.e., a seed.
Algorithm 2 describes the proposed method for rearranging chunks according
to step 2. This algorithm assumes that the process uses l thread blocks of sizeb and three inputs: (1) 2m chunks that need to be rearranged, (2) the number
l of circular permutations, and (3) their representative elements p0 , p1 , . . . , pl−1 ,where pj is the representative element for circular permutation σj (0 ≤ j < l).
Given these inputs, the thread block with index j loads a chunk that corresponds
to the representative element pj in a register (line 2). Next, another chunk that
corresponds to the next element in circular permutation σj is moved into the empty
memory location where the copied chunk was present (line 5). Similar chunkmovements are iterated until the representative element is reached (lines 4–7).
Finally, the algorithm moves the representative element pj from the register into
the last empty location (line 8). Note that this rearrangement is optimal in terms
16

of the amount of memory access because all of the chunks that must move to adifferent location are loaded once and stored once.
Algorithm 2 Chunk rearrangement kernel based on a product of circular permutations.
Input: Chunks C0 , C1 , . . . , C2m−1 that need to be rearranged, number 2m ofchunks, and representative elements p0 , p1 , . . . , pl−1 of circular permutations
Output: Rearranged chunks C0 , C1 , . . . , C2m−1
1: j ← thread block index;
2: W ← Cpj ;3: i ← pj ;4: while 2i mod (2m − 1) ̸= pj do
5:
Ci ← C2i mod (2m−1) ;6:i ← 2i mod (2m − 1);
7: end while
8: Ci ← W ;▷ b elements are copied in parallel▷ b elements are copied in parallel▷ b elements are copied in parallel

4.2. Precomputation of Representative Elements
In order to perform the in-place DWT, the GPU threads must compute thenumber l of circular permutations and all the representative elements p0 , p1 , . . . , pl−1 ,for an arbitrary problem size n. In particular, the representative elements must be
assigned to thread blocks without duplication. However, it is not easy to parallelize this assignment to GPU threads, mainly due to the lack of a global synchronization mechanism.
Our solution to this issue is to precompute the number l of circular permutations and representative elements p0 , p1 , . . . , pl−1 on the CPU. Algorithm 3 describes our sequential algorithm for computing all the representative elements.
This algorithm enumerates all of the elements in a circular permutation σj and17stores the minimum value as pj . The representative elements are computed per
circular permutation, so an array X is used to avoid duplicated enumeration.
Algorithm 3 Precomputation of representative elements for circular permutations.
Input: Number 2m of chunks
Output: Number l of circular permutations and representative elements
p0 , p1 , . . . , pl−1 of circular permutations
1: for i ← 0 to 2m − 1 do
▷ Initialization
2:
Xi ← 0;3: end for
4: l ← 0;
5: for i ← 0 to 2m − 1 do
▷ Compute pl6:
if Xi = 0 then▷ Not registered yet
7:pl ← i;8:
Xi ← 1;9:j ← 2i mod (2m − 1); ▷ Generate number j such that j and i belong

to the same circular permutation
10:
while j ̸= i do
11:
Xj ← 1;12:j ← 2j mod (2m − 1);13:
end while14:l ← l + 1;
15:
end if
16: end for
The representative elements p0 , p1 , . . . , pl−1 are then transferred to constantmemory [18] (i.e., GPU memory) before generating chunks on the GPU. Note thatthe capacity of constant memory is 64 KB in the current CUDA. This capacity is
sufficiently large to store representative elements for the largest problem size nthat can be processed on a GPU. A detailed analysis is provided in Section 5.1.

18

4.3. Space Complexity Analysis
The proposed method requires two different working spaces to perform the
DWT: (1) constant memory used for representative elements, and (2) shared memory and registers used for chunk generation and circulation.
The constant memory usage depends on the number l of representative elements. The following lemma [27] is introduced to formulate l for arbitrary n.
Lemma 2. According to Burnside’s lemma [28], the number of necklaces with
M beads and K (> 1) colors (see Fig. 4), where rotations are considered to beequivalent, is
( )
M
1 ∑
M
1 ∑ gcd(i,M )
µ
Ki =
K
,
(5)
M
i
M i=1
i|M

where µ is the classical Möbius function and gcd(i, M ) is the greatest common
divisor of i and M .
Thus, the number l+2 is equivalent to that of necklaces with M = ⌈log 2m⌉ beads
and K = 2 colors, where rotations are considered to be equivalent. Note that “+2”
considers two exceptions, permutations (0 0 . . . 0) and (1 1 . . . 1), which nevermove chunks during the arrangement process. These exceptions correspond to thehead and tail chunks. Lemma 2 gives the following theorem.
Theorem 1. The number l of representative elements is smaller than m = ⌈n/2b⌉
if m ≥ 128.
Proof . Eq. (5) can be simplified further as follows.
1 ∑ gcd(i,M )
1 ∑ i
l+2=
K
≤
K
M i=1
M i=1
M

M

K(K M − 1)
M (K − 1)
K M +1<M (K − 1)
=

19

(6)
(7)
(8)

M = ⌈log 2m⌉, K = 2 and x ≤ ⌈x⌉ ≤ x + 1 for an arbitrary real number x, so
we have8m
log 2m≤ m, if m ≥ 128.

l <

(9)
(10)
□

Eq. (10) indicates that 2m chunk indexes are reduced to at least half because acircular permutation comprises at least two elements (i.e., two chunk indexes).
Therefore, the number l of representative elements is smaller than m = ⌈n/2b⌉.
Asymptotically, l is bounded by O(n/ log n) from Eq. (9).
In addition, the proposed method generates O(n/ log n) thread blocks during
chunk circulation because the method uses l thread blocks in total. The usageof registers per thread block (i.e., per circular permutation) is equivalent to the
chunk size, where the chunk size is 4b in bytes if the signal data is stored in a
4-byte data type. Therefore, the proposed algorithm consumes O(n/ log n) registers if all of the thread blocks are processed simultaneously. However, due tothe limited amount of computational resources, such as the shared memory andregisters, CUDA restricts the maximum number of active thread blocks that canbe processed simultaneously [18]. Therefore, our CUDA-based implementation
runs with O(1) registers. Note that the proposed method achieves this restriction
by data arrangement, which eliminates the data dependencies between different
thread blocks. Such independent thread blocks are necessary to reduce the number of active thread blocks.
In summary, the proposed method consumes ⌈n/2b⌉ additional space locations
20

to perform the lifting-based DWT. Currently, we use b = 512, so this additional
memory usage is an order of magnitude lower than the n needed by existing methods that distinguish the output buffer from the input buffer.
4.4. Tiling-based Transformation
Tiling-based approaches [29, 30] are useful for dealing with large-scale datathat exceeds the GPU memory capacity. These approaches split the data into
tiles that are transformed separately. However, Eqs. (1) and (2) indicate that eachelement needs its neighbors to update its value. Consequently, separated transformation creates blocking artifacts at the tile boundaries, thereby excessive tiles candegrade the transformation quality.
This degradation can be avoided by having halo regions for each tile (see Fig.
5); the input data is decomposed into small overlapping segments (with halo regions) to exactly transform all elements inside the tiles. Both the proposed andLaan’s previous methods can be enhanced with a tiling-based approach by specifying a segment as the input. However, the proposed memory-saving method accepts roughly twice larger segments than those accepted by the previous method.
As for tile creation, a 1-D block decomposition scheme is deployed to minimize the overhead of data decomposition. An alternative is a 2-D block decomposition scheme, but it complicates data transfer between the CPU and GPU because
elements of a 2-D tile is usually scattered throughout memory. In contrast, elements of a 1-D tile is stored in a contiguous memory region, which facilitates datatransfer between the CPU and GPU.

21

…

Halo

Tile

Segment…
Figure 5: A tiling-based approach with halo regions.
Table 2: Specifictions of hardware devices deployed for experiments. Effective bandwidths were
measured using the bandwidthTest program [18]. DMA stands for direct memory access.

GPU card
980 Ti
Compute capability
5.2
Architecture
Maxwell
GPU memory capacity (GB)6
Peak memory bandwidth (GB/s)
336.5
Effective memory bandwidth (GB/s)
249.1
Peak PCIe bandwidth (GB/s)
15.8
Effective PCIe bandwidth (GB/s)
11.5
Number of DMA engines
2780 Ti
680
580
3.5
3.0
2.0
Kepler Kepler Fermi
32
2
336.0 192.2 192.4
236.3 152.5 169.2
15.8
15.8
8.0
11.8
12.2
6.3
1
1
1

Note that tiling-based approaches can be efficiently processed in a pipeline;
the data segments are streamed through the pipeline to overlap CPU-GPU datatransfer with kernel execution. Such a pipeline can be implemented with CUDA
streams [18].

5. Experimental Results
To evaluate the proposed method, we performed comparisons with Laan’s previous method [10], which enabled GPU-friendly memory accesses with a blocklayout scheme (see Section 3.1). In the experiments, we investigated the mem22

ory usage and execution time for the Deslauriers-Dubuc (13,7) wavelet [31] andN × N -pixel images, where 2K ≤ N ≤ 32K. Each pixel comprised 4-byte data.
Table 2 summarizes the specifications for our experimental machine, which had
an Intel Core i7-3770K CPU and either a GeForce GTX 980 Ti, 780 Ti, 680 or
580. The machine ran with Windows 8.1 (64 bit), CUDA 7.5 [18] and the driver
version 355.82. All of the implementations were compiled using the Visual Studio2013 with the -O2 option.
5.1. Memory Consumption
Figure 6 shows the amount of memory usage. Compared with the previous
method, the proposed method halved the amount of memory usage by buffer unification. Consequently, on a GeForce GTX 980 Ti GPU, the proposed method
processed an N × N -pixel image of up to N = 32K pixels without data decomposition. By contrast, the previous method processed a maximum size of N = 16K
pixels without data decomposition. Note that the entire GPU memory cannot be
used for the input/output buffer because part of the GPU memory is reserved forthe frame buffer.
Next, we investigated the additional memory size α (= 4l) required to store allof the representative elements (Fig. 7). As shown in this figure, the size α reached3.09 KB for a large image of 32K × 32K pixels, which consumed 4 GB of the
GPU memory. Consequently, we consider that the proposed method is useful forminimizing additional memory usage with a practical data size, although it does
not achieve in-place transformation.

23

8,192

Memory usage (MB)10,000Previous
Proposed

1,000

2,048
1,024

128
100

4,096

512
256

32
64

10

16

1
2K

4K

8K
Image size N

16K

32K

Figure 6: Comparison of memory usage with different image sizes N . An image comprises N ×N
pixels.

Working memory size (KB)

3.5

3.09

3
2.5
1.75

2
1.5
0.97
1
0.5

0.56
0.30

0
2K

4K

8K
Image size N

16K

32K

Figure 7: Working memory size required for representative elements. An image comprises N × N
pixels.

24

5.2. DWT Performance
The DWT is part of a processing pipeline in a practical application. Therefore, many GPU-accelerated applications typically implement the DWT as wellas other pipeline stages on the GPU to accelerate the overall pipeline. Consequently, we measured the execution time for four configurations where the input
data and output data were stored in either the CPU memory or the GPU memory. In the following, the notation x → y denotes a configuration, where x andy represent the locations of the input data and the output data, respectively, andx, y ∈ {C, G}. The notations C and G correspond to the CPU memory and the
GPU memory, respectively.
Figure 8 shows the execution times using the proposed method (without tiling)
and the previous method (with tiling) for large data. Note that the previous method
was extended to transform this large image successfully, where we decomposed
the image into two segments to transform them in order, as presented in Section
4.4. However, the pipelining capability was switched off, so that data segments
were processed in a synchronous manner.
The proposed method minimized the amount of memory usage, and it avoided
data transfer between the CPU and GPU. Thus, the proposed method was 2.4–
3.9 times faster than the tiled previous method when both the input and output
data were stored in the GPU memory (G → G). However, the proposed method
requires a preprocessing phase, which increased the kernel execution time by 2.0–
2.6 times. Similarly, the proposed method achieved speedups of 1.2–1.8 times
compared with the tiled previous method when either the input data or output
25

600

Kernel execution

Execution time (ms)

Data transfer

Data transfer

G→C

400
300
200
100

C→G

G→G

C→C

C→G

C→C

G→G

(c)

G→C

Proposed

Prev. (tiled)Prev. (tiled)
C→G

Proposed

Prev. (tiled)

Proposed

Proposed

Prev. (tiled)Prev. (tiled)

Proposed

Prev. (tiled)

Proposed

Prev. (tiled)

Proposed0

Prev. (tiled)50

Proposed

100

G→C

C→C

Kernel execution

Prev. (tiled)150

Data transfer

Proposed

200

500
450
400
350
300
250
200
150
100
50
0

Prev. (tiled)

Kernel execution

250

G→G

C→G

(b)
Execution time (ms)

Data transfer

Proposed

G→C

(a)
300

Prev. (tiled)

Proposed

Prev. (tiled)

Proposed

Proposed

Prev. (tiled)Prev. (tiled)

Proposed

Proposed

Prev. (tiled)Prev. (tiled)
G→G

Execution time (ms)

Kernel execution

500

0
Proposed

Execution time (ms)

1000
900
800
700
600
500
400
300
200
100
0

C→C

(d)

Figure 8: Execution time with different input/output configurations for large data. Results for (a)GeForce GTX 980 Ti with N = 32K, (b) 780 Ti with N = 24K, (c) 680 with N = 16K, and (d)580 with N = 16K. The data sizes were 4.0 GB, 2.3 GB and 1.0 GB for N = 32K, N = 24K
and N = 16K, respectively.

26

data were stored in the GPU memory (G → C and C → G). These configurations required data transfer between the CPU and GPU regardless of the data
size. Moreover, the tiled previous method swapped half of the input data before
the second transformation because its input and output buffers exhausted the GPU
memory. This swapping operation increased the data transfer volume, and thus thetiled previous method was slower than the proposed method. Finally, the proposed
method was 12–18% slower than the tiled previous method when both the input
and output data were stored in the CPU memory (C → C). However, CPU–GPU
data transfer determines the transformation throughput for this configuration, sothat these overheads can be partially hidden by pipelined execution, as describedin Section 5.4.
Figure 9 shows the execution times for small data, which could be transformed
without tiling. Compared with the previous method, the proposed method required
1.1–2.7 more time to transform these small data. In this case, the previous methodcompleted the DWT in almost the same data transfer time as the proposed method.
Consequently, the preprocessing overheads increased the execution time for the
proposed method. As mentioned for large data, these overheads can be partiallyaddressed by pipelined execution.
5.3. Kernel Execution Efficiency
The DWT is a memory-bound operation, so we evaluated the execution efficiency of our kernel functions in terms of the effective memory throughput. Figure
10 shows the effective memory throughput E = M/T for the proposed method,

27

250

250
Kernel execution

Data transfer

200

Execution time (ms)

Execution time (ms)

Data transfer

150
100
50

150
100
50

Proposed
Prev.
Prev. (tiled)

Proposed
Prev.
Prev. (tiled)

Proposed
Prev.
Prev. (tiled)

Proposed
Prev.
Prev. (tiled)

Proposed
Prev.
Prev. (tiled)

Proposed
Prev.
Prev. (tiled)

Proposed
Prev.
Prev. (tiled)

0

Proposed
Prev.
Prev. (tiled)

0

G→G

G→C

C→G

C→C

G→G

G→C

C→G

C→C

(a)

(b)

70

120
Kernel execution

Data transfer

20

60
40

G→G

G→C

C→G

Prev. (tiled)

Prev.

Proposed

Prev. (tiled)

Prev.

Proposed

Prev. (tiled)

Prev.

Proposed0

Prev. (tiled)

0

Prev.

20

Proposed

10

C→C

(c)

Proposed
Prev.
Prev. (tiled)3080

Proposed
Prev.
Prev. (tiled)40

Proposed
Prev.
Prev. (tiled)

Execution time (ms)50

Kernel execution

100

Proposed
Prev.
Prev. (tiled)

Data transfer
60

Execution time (ms)

Kernel execution

200

G→G

G→C

C→G

C→C

(d)

Figure 9: Execution time with different input/output configurations for small data. Results for (a)GeForce GTX 980 Ti with N = 16K, (b) 780 Ti with N = 16K, (c) 680 with N = 8K, and(d) 580 with N = 8K. The data sizes were 1 GB and 256 MB for N = 16K and N = 8K,
respectively.

28

Table 3: Number of data elements accessed. We used r = 256 and h = 6 for the DeslauriersDubuc (13,7) wavelet employed in this study.

Breakdown
Chunk generation
Chunk rearrangement
Lifting DWTAmount of memory access
4N 2
4N (N − 2b)4(1 + h/r)N 2

where M and T are the amount of memory accesses and the kernel execution time,
respectively.
The amount M of memory accesses is determined as follows. First, the number of elements accessed in the chunk generation step is given by 2N 2 because
every chunk is read and written once during this step. This step is repeated twice,each for the horizontal 1-D DWT and the vertical 1-D DWT, so that 4N 2 memory
accesses are required for an N × N -pixel image. Similarly, the chunk rearrangement step incurs 2N (N − 2b) memory accesses, where 2b corresponds to the head
and tail chunks, which never move during this step, and thereby 4N (N − 2b) locations are accessed for an N × N -pixel image. Finally, the thread blocks in thelifting step access halo regions because computation for an element requires its
neighboring elements. For simplicity, we considered 1-D data with a size of n.
For 1-D data, n/r thread blocks referred to r + 2h elements to output n elements,
where r and h = k0 + k1 are the numbers of responsible elements per thread blockand neighboring elements in the halo region, respectively. Consequently, the lifting step for 1-D data accessed (r + 2h)n/r + n = 2(1 + h/r)n memory locations.
By extending this result to 2-D data, which requires a row-wise 1-D DWT and acolumn-wise 1-D DWT, we obtained 4(1 + h/r)N 2 .

29

Memory throughput (GB/s)

300
250
200
150
100
980 Ti50

780 Ti

680

580

0
2K

4K

8K
Image size N

16K

32K

Figure 10: Effective memory throughputs obtained without tiling. Owing to memory exhaustion,
kernels on 780 Ti, 680 and 580 cards failed to transform a 32K × 32K-pixel image.

As shown in Fig. 10, the effective memory throughputs were bounded by thepeak GPU memory bandwidth of the deployed card (Table 2). The proposed
method achieved high memory throughputs, although it changes the order ofdata elements in the preprocessing phase. This highly efficient arrangement is
attributable to regular memory access in contiguous regions. However, Table 3
shows the main drawback of the proposed method, which is the increased amountof memory accesses. Compared with the previous method, the proposed method
required three times more memory accesses to perform the lifting-based DWTwith a unified buffer. This theoretical result explains why the proposed method
required 2.0–2.6 times longer kernel time than the previous method (see Section
5.2).

30

5.4. Pipelined DWT Throughput
We next evaluated tiling-based pipelined versions presented in Section 4.4.
We measured the transformation throughput for configuration C → C, where weassumed that an image of size N × N was stored in the CPU memory. The image
size N ranged from 2K to 64K, so that the largest image could not be stored
entirely in the GPU memory.
Figure 11 shows the throughputs measured for different image sizes N , where
the number of tiles was set to N/b. Except for the 980 Ti, the throughputs
achieved by the pipelined proposed method were close to those achieved using thepipelined previous method. By contrast, the non-pipelined proposed version was
12–25% slower than the non-pipelined previous version due to the preprocessing phase. Thus, the proposed method degraded the transformation throughput
for on-memory data, but this performance degradation can be partially hidden byoverlapping, particularly for large-scale data that cannot be stored entirely in the
GPU memory.
With respect to the 980 Ti card, this card is equipped with two direct memory
access (DMA) engines, which enable simultaneous bidirectional transfer between
the CPU and GPU. Consequently, as compared with the non-pipelined version, thepipelined version reduced the data transfer time. However, this bidirectional transfer caused memory access contention, so there was a gap between the pipelined
and non-pipelined throughputs, as shown in Fig. 11(a).
Figure 12 investigates this contention issue in more detail. This figure explains(1) why the pipelined proposed version was slower than the pipelined previous
31

Transformation throughput (Gpixel/s)

Transformation throughput (Gpixel/s)

2.5
2.0
1.5
1.0
0.5

Previous (no-pipe)

Proposed (no-pipe)

Previous (pipe)

Proposed (pipe)

0.0
2K

4K

8K
16K
Image size N

32K

1.6
1.4
1.2
1.0
0.8
0.6
0.4

Previous (no-pipe)

Proposed (no-pipe)

0.2

Previous (pipe)

Proposed (pipe)

4K

32K

0.0

64K

2K

2.0
1.5
1.0
0.5

Previous (no-pipe)

Proposed (no-pipe)

Previous (pipe)

Proposed (pipe)

0.0
2K

4K

64K

(b)
Transformation throughput (Gpixel/s)

Transformation throughput (Gpixel/s)

(a)

8K
16K
Image size N

8K
16K
Image size N

32K

1.0
0.8
0.6
0.4
0.2

(c)

Proposed (no-pipe)

Previous (pipe)

Proposed (pipe)

4K

32K

0.0
2K

64K

Previous (no-pipe)

8K
16K
Image size N

64K

(d)

Figure 11: Comparison of transformation throughputs with different image sizes N . Results for(a) GeForce GTX 980 Ti, (b) 780 Ti, (c) 680, and (d) 580. The number of segments was fixed to
N/b.

32

(a)

(b)

(c)

(d)
Figure 12: Timeline views of tiling-based pipelined implementations. In each subfigure, the upper
four lines show the activities of CUDA streams. The lower three lines summarize the activities interms of kernel execution (gray), CPU→GPU data transfer (green) and GPU→CPU data transfer(purple). GeForce GTX 980 Ti has two DMA engines, so that (a) the kernel execution overheadof the proposed method increased data transfer time compared to (b) the previous method. In
contrast, GTX 680 has a single DMA engine, so that there was no significant gap between (c) the
proposed method and (d) the previous method.

33

version on the GTX 980 Ti card and (2) why the former was as fast as the latter on the GTX 680 card. With multiple DMA engines, the GPU memory canbe simultaneously accessed by three CUDA operations: (1) kernel execution, (2)CPU→GPU data transfer and (3) GPU→CPU data transfer, as shown in Figs.
12(a) and 12(b). This simultaneous access slightly affected the performance ofthe pipelined proposed version, which took longer kernel execution time than thepipelined previous version. In fact, the proposed version declined the data transfer
rate by 10% compared to the previous version. Therefore, the declined rate determined the transformation time though kernel execution was fully overlapped with
CPU-GPU data transfer. In contrast, such a decline was not observed on the GTX
680 card, which exclusively transfers data between the CPU and GPU owing to asingle DMA engine; at most two CUDA operations can access the GPU memory
simultaneously: (1) kernel execution and (2) either CPU→GPU data transfer or
GPU→CPU data transfer, as shown in Figs. 12(c) and 12(d).
Thus, the pipelined DWT throughput for configuration C → C was mainly
determined by the peak PCIe bandwidth. Because the 580 card has a PCIe generation 2 interface, its DWT throughput was half of that achieved with the 780 ti or
680 cards, which has a PCIe generation 3 interface. Multiple DMA engines further increased the throughput with a full overlap but declined the data transfer rate
due to memory contention. The pipelined scheme is useful for partially hiding theoverhead of the proposed method.
As for the practicality of the proposed method, the method is useful for processing larger data at a time. For example, a decoding system [6] for satellite34

images adopted a series of decoder accelerated on a GPU. Another example isa clustering system [7] that groups similar objects into disjoint classes. Because
a cluster algorithm was applied to tiles separately, the proposed memory-saving
method allows more objects to be grouped from an approximately twice larger
tile.

6. Conclusion
In this study, we proposed a parallel DWT method that achieves both a compact data representation and GPU-friendly memory accesses for the lifting scheme.
The proposed data representation unifies the output buffer with the input buffer.
The key idea of unification is chunk-based data rearrangement, which separates
data-independent tasks from data-dependent tasks, and thus massively parallel
threads are allowed to move data elements according to circular permutations. The
proposed method also achieves highly efficient memory access by allowing GPU
threads to access contiguous memory regions. However, the proposed method requires at most ⌈n/2B⌉ numbers as seeds for a series of memory addresses that areneeded for parallel rearrangement.
The results of experiments showed that the proposed method was a maximumof 3.9 times faster than a previous method that separates the input buffer fromthe output buffer. In particular, the proposed method was useful when either the
input buffer or output buffer existed in the GPU memory. We also found that the
proposed method achieved a high memory throughput that was close to the peakmemory bandwidth. However, our kernel was affected by the increased amount of35

memory accesses, i.e., three times that of the previous kernel. This disadvantage
can be partially resolved by overlapping kernel execution with CPU-GPU datatransfer. This pipelined implementation achieved a transformation throughput thatwas comparable to the previous method.

Acknowledgments
This study was supported in part by the Japan Society for the Promotion of
Science KAKENHI Grant Numbers 15K12008 and 15H01687, and the Japan
Science and Technology Agency CREST program, “An Evolutionary Approachto Construction of a Software Development Environment for Massively-Parallel
Computing Systems.” We are also grateful to the anonymous reviewers for theirvaluable comments.

References[1] S. Mallat, A Wavelet Tour of Signal Processing: The Sparse Way, 3rd Edition, Academic Press, 2008.
[2] I. V. Šimek, R. R. ASN, GPU acceleration of 2D-DWT image compressionin MATLAB with CUDA, in: Proc. 2nd UKSIM European Symp. Computer
Modeling and Simulation (EMS’08), 2008, pp. 274–277.
[3] J. Matela, GPU-based DWT acceleration for JPEG2000, in: Proc. Doctoral
Workshop Mathematical and Engineering Methods in Computer Science(MEMICS’09), 2009, pp. 136–143.
36

[4] J. Ao, S. Mitra, B. Nutter, Fast and efficient lossless image compression
based on CUDA parallel wavelet tree encoding, in: Proc. Southwest Symp.
Image Analysis and Interpretation (SSIAI’14), 2014, pp. 21–24.
[5] L. Zhu, Y. Zhou, D. Zhang, D. Wang, H. Wang, X. Wang, Parallel multilevel 2D-DWT on CUDA GPUs and its applications in ring artifact removal,
Concurrency and Computation: Practice and Experience.
[6] C. Song, Y. Li, B. Huang, A GPU-accelerated wavelet decompression system with SPIHT and Reed-Solomon decoding for satellite images, IEEE J.
Selected Topics in Applied Earth Observations and Remote Sensing 4 (3)(2011) 683–690.
[7] A. A. Yıldırım, C. Özdoğan, Parallel wavecluster: A linear scaling parallelclustering algorithm implementation with application to very large datasets,
J. Parallel and Distributed Computing 71 (7) (2011) 955–962.
[8] W. Sweldens, The lifting scheme: A construction of second generation
wavelets, SIAM J. Mathematical Analysis 29 (2) (1998) 511–546.
[9] S. G. Mallat, A theory for multiresolution signal decomposition: The
wavelet representation, IEEE Trans. Pattern Analysis and Machine Intelligence 11 (7) (1989) 674–693.
[10] W. J. van der Laan, A. C. Jalba, J. B. Roerdink, Accelerating wavelet liftingon graphics hardware using CUDA, IEEE Trans. Parallel and Distributed
Systems 22 (1) (2011) 132–146.
37

[11] M. E. Angelopoulou, P. Y. K. Cheung, K. Masselos, Y. Andreopoulos, Implementation and comparison of the 5/3 lifting 2D discrete wavelet transform computation schedules on FPGAs, J. Signal Processing Systems 51 (1)(2008) 3–21.
[12] D. A. Bader, V. Agarwal, S. Kang, Computing discrete transforms on the
Cell Broadband Engine, Parallel Computing 35 (3) (2008) 119–137.
[13] M. Błażewicz, M. Ciżnicki, P. Kopta, K. Kurowski, P. Lichocki, Twodimensional discrete wavelet transform on large images for hybrid computing architectures: GPU and CELL, in: Proc. 17th European Conf. Parallel
Computing Workshops (Euro-Par’11 Workshops), Part I, 2011, pp. 481–490.
[14] NVIDIA Corporation, NVIDIA GeForce GTX 980 (Nov. 2014).
URL

http://international.download.nvidia.com/

geforce-com/international/pdfs/GeForce_GTX_980_
Whitepaper_FINAL.PDF[15] J. D. Owens, M. Houston, D. Luebke, S. Green, J. E. Stone, J. C. Phillips,
GPU computing, Proceedings of the IEEE 96 (5) (2008) 879–899.
[16] F. Ino, Y. Munekawa, K. Hagihara, Sequence homology search using finegrained cycle sharing of idle GPUs, IEEE Trans. Parallel and Distributed
Systems 23 (4) (2012) 751–759.
[17] Y. Okitsu, F. Ino, K. Hagihara, High-performance cone beam reconstruction

38

using CUDA compatible GPUs, Parallel Computing 36 (2/3) (2010) 129–141.
[18] NVIDIA Corporation, CUDA C Programming Guide Version 7.5 (Sep.
2015).
URL

http://docs.nvidia.com/cuda/pdf/CUDA_C_

Programming_Guide.pdf[19] T.-T. Wong, C.-S. Leung, P.-A. Heng, J. Wang, Discrete wavelet transformon consumer-level graphics hardware, IEEE Trans. Multimedia 9 (3) (2007)
668–673.
[20] D. Shreiner, M. Woo, J. Neider, T. Davis, OpenGL Programming Guide, 5th
Edition, Addison-Wesley, Reading, MA, 2005.
[21] C. Tenllado, J. Setoain, M. Prieto, L. Piñuel, F. Tirado, Parallel implementation of the 2D discrete wavelet transform on graphics processing units: Filter bank versus lifting, IEEE Trans. Parallel and Distributed Systems 19 (3)(2008) 299–310.
[22] J. Franco, G. Bernabé, J. Fernández, M. Ujaldón, The 2D wavelet transformon emerging architectures: GPUs and multicores, J. Real-Time Image Processing 7 (3) (2012) 145–152.
[23] V. Galiano, O. López, M. P. Malumbres, H. Migallón, Parallel strategies
for 2D discrete wavelet transform in shared memory systems and GPUs, J.
Supercomputing 64 (1) (2013) 4–16.
39

[24] M. Kucis, D. Barina, M. Kula, P. Zemcik, 2-D discrete wavelet transformusing GPU, in: Proc. 26th Int’l Symp. Computer Architecture and High Performance Computing Workshop (SBAC-PADW’14), 2014, pp. 1–6.
[25] B. Sharma, N. Vydyanathan, Parallel discrete wavelet transform using
the open computing language: a performance and portability study, in:Proc. 24th IEEE Int’l Parallel and Distributed Processing Symp. Workshops(IPDPSW’10), 2010.
[26] Khronos OpenCL Working Group, The OpenCL specification version 1.1,
http://www.khronos.org/registry/cl/ (Jun. 2011).
[27] N. Metropolis, G.-C. Rota, Witt vectors and the algebra of necklaces, Advances in Mathematics 50 (2) (1983) 95–125.
[28] W. Burnside, Theory of Groups of Finite Order, 1897.
[29] T. Acharya, P.-S. Tsai, JPEG2000 Standard for Image Compression: Concepts, Algorithms and VLSI Architectures, Wiley-Interscience, 2004.
[30] J.-W. Kim, J. Song, S. Lee, I.-C. Park, Tiled interleaving for multi-level 2D discrete wavelet transform, in: Proc. 40th IEEE Int’l Symp. Circuits andSystems (ISCAS’07), 2007, pp. 3984–3987.
[31] G. Deslauriers, S. Dubuc, Symmetric iterative interpolation processes, Constructive Approximation 5 (1) (1989) 49–68.

40