WO2008020839A2 - Cache friendly method for performing inverse discrete wavelet transform - Google Patents

Cache friendly method for performing inverse discrete wavelet transform Download PDF

Info

Publication number
WO2008020839A2
WO2008020839A2 PCT/US2006/031878 US2006031878W WO2008020839A2 WO 2008020839 A2 WO2008020839 A2 WO 2008020839A2 US 2006031878 W US2006031878 W US 2006031878W WO 2008020839 A2 WO2008020839 A2 WO 2008020839A2
Authority
WO
Grant status
Application
Patent type
Prior art keywords
data
input data
output data
method
memory
Prior art date
Application number
PCT/US2006/031878
Other languages
French (fr)
Other versions
WO2008020839A3 (en )
Inventor
Jagadeesh Sankaran
Original Assignee
Texas Instruments Incorporated
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING; COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/14Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/148Wavelet transforms
    • HELECTRICITY
    • H04ELECTRIC COMMUNICATION TECHNIQUE
    • H04NPICTORIAL COMMUNICATION, e.g. TELEVISION
    • H04N19/00Methods or arrangements for coding, decoding, compressing or decompressing digital video signals
    • H04N19/42Methods or arrangements for coding, decoding, compressing or decompressing digital video signals characterised by implementation details or hardware specially adapted for video compression or decompression, e.g. dedicated software implementation
    • H04N19/423Methods or arrangements for coding, decoding, compressing or decompressing digital video signals characterised by implementation details or hardware specially adapted for video compression or decompression, e.g. dedicated software implementation characterised by memory arrangements
    • HELECTRICITY
    • H04ELECTRIC COMMUNICATION TECHNIQUE
    • H04NPICTORIAL COMMUNICATION, e.g. TELEVISION
    • H04N19/00Methods or arrangements for coding, decoding, compressing or decompressing digital video signals
    • H04N19/60Methods or arrangements for coding, decoding, compressing or decompressing digital video signals using transform coding
    • H04N19/63Methods or arrangements for coding, decoding, compressing or decompressing digital video signals using transform coding using sub-band based transform, e.g. wavelets

Abstract

This invention is a method for inverse Wavelet transform using a breadth- first output data calculation which uses input data to calculate at least one output data for each iteration of a software loop even if the same input data is used in a later iteration for calculating other output data. This reduces data movement between memory and the data processor core thus reducing the possibility of cache misses and memory stalls due to access conflicts. The input data and computed output data are preferably stored as subwords packed within data words in memory. In inverse Wavelet transformation, this method performs vertical spatial frequency expansion and horizontal spatial frequency expansion for each level of Wavelet encoding. This invention arranges data flow providing a more efficient use of memory bandwidth and cache space than other known methods. A method for inverse Wavelet transform comprises the steps of storing input data as subwords packed within data words in a memory; calculating output data in a software loop in a breadth-first fashion by recalling input data (902) and filter coefficients (903, 905) for computation of at least one output data for each iteration of the software loop even if the same input data is used in a later iteration for calculating other output data (904, 906); and storing calculated output data as subwords packed within data words in the memory.

Description

CACHE FRIENDLY METHOD FOR PERFORMING

INVERSE DISCRETE WAVELET TRANSFORM

The technical field of this invention is digital data encoding and, more particularly, inverse discrete wavelet transform coding. BACKGROUND

Wavelet encoding of image data transforms the image from a pixel spatial domain into a mixed frequency and spatial domain. In the case of image data the wavelet transformation includes two dimensional coefficients of frequency and scale. FIGS. 1 to 6 illustrate the basic technique of wavelet image transformation. The two dimensional array of pixels is analyzed X and Y directions and a set for transformed data that can be plotted in respective X and Y frequency. FIG. 1 illustrates transformed data 100 with the upper left corner the origin of the X and Y frequency coordinates. This transformed data is divided into four quadrant subbands. Quadrant 101 includes low frequency X data and low frequency Y data denoted as LL. Quadrant 102 includes low frequency X data and high frequency Y data denoted LH. Quadrant 103 includes high frequency X data and low frequency Y data denoted HL. Quadrant 104 includes high frequency X data and high frequency Y data denoted HH.

Organizing the image data in this fashion with a wavelet transform permits exploitation of the image characteristics for data compression. It is found that most of the energy of the data is located in the low frequency bands. The image energy spectrum generally decays with increasing frequency. The high frequency data contributes primarily to image sharpness. When describing the contribution of the low frequency components the frequency specification is most important. When describing the contribution of the high frequency components the time or spatial location is most important. The energy distribution of the image data may be further exploited by dividing quadrant 101 into smaller bands. FIG. 2 illustrates this division. Quadrant 101 is divided into subquadrant 111 denoted LLLL, subquadrant 112 denoted LLLH, subquadrant 113 denoted LLHL and subquadrant 114 denoted LLHH. As before, most of the energy of quadrant 101 is found in subquadrant 111. FIG. 3 illustrates a third level division of subquadrant 111 into subquadrant 121 denoted LLLLLL, subquadrant 122 denoted LLLLLH, subquadrant 123 denoted LLLLHL and subquadrant 124 denoted LLLLHH. FIG. 4 illustrates a fourth level division of subquadrant 121 into subquadrants 131 denoted LLLLLLLL, subquadrant 132 denoted LLLLLLLH, subquadrant 133 denoted LLLLLLHL and subquadrant 134 denoted LLLLLLHH.

For an n-level decomposition of the image, the lower levels of decomposition correspond to higher frequency subbands. Level one represents the finest level of resolution. The n-th level decomposition represents the coarsest resolution. Moving from higher levels of decomposition to lower levels corresponding to moving from lower resolution to higher resolution, the energy content generally decreases. If the energy content of level of decomposition is low, then the energy content of lower levels of decomposition for corresponding spatial areas will generally be smaller. There are spatial similarities across subbands. A direct approach to use this feature of the wavelet coefficients is to transmit wavelet coefficients in decreasing magnitude order. This would also require transmission of the position of each transmitted wavelet coefficient to permit reconstruction of the wavelet table at the decoder. A better approach compares each wavelet coefficient with a threshold and transmits whether the wavelet value is larger or smaller than the threshold. Transmission of the threshold to the detector permits reconstruction of the original wavelet table. Following a first pass, the threshold is lowered and the comparison repeated. This comparison process is repeated with decreasing thresholds until the threshold is smaller than the smallest wavelet coefficient to be transmitted. Additional improvements are achieved by scanning the wavelet table in a known order, with a known series of thresholds. Using decreasing powers of two seems natural for the threshold values.

FIG. 5 illustrates scale dependency of wavelet coefficients. Each wavelet coefficient has a set of four analogs in a next lower level. In FIG. 19, wavelet coefficients B, C and D are shown with corresponding quads Bl, B2, B3, B4, Cl, C2, C3, C4, Dl, D2, D3 AND D4. Each of these wavelet coefficients has a corresponding quad in the next lower level. As shown in FIG. 5: wavelet coefficients Bl, B2, B3 and B4 each have corresponding quads Bl, B2, B3 and B4 in the next lower level; wavelet coefficient C2 has a corresponding quad C2; and wavelet coefficient D2 has a corresponding quad D2.

FIG. 6 illustrates an example Morton scanning order and a corresponding set of wavelet coefficients used in a coding example. The Morton scanning retains the same order for each 4 by 4 block in increasing scale. In the wavelet coefficient example H indicates the coefficient is greater than the threshold and L indicates the coefficient is less than the threshold. The encoded data in this example is "HHLH HLLH LLLL HLHL."

Use of the wavelet compressed data requires a reversal of the encoding process called an inverse discrete wavelet transform. SUMMARY

This invention is a method for inverse Wavelet transform using a breadth-first output data calculation. This breadth-first computation results in fewer intermediate variables that need to be stored. Such a breadth-first approach used input data to calculate at least one output data for each iteration of a software loop even if the same input data is used in a later iteration for calculating other output data. This reduces data movement between memory and the data processor core thus reducing the possibility of cache misses and memory stalls due to access conflicts. The input data and computed output data are preferably stored as subwords packed within data words in memory. In inverse Wavelet transformation this method performs vertical spatial frequency expansion and horizontal spatial frequency expansion for each level of Wavelet encoding. The software loop schedules multiplications between input data and corresponding non-zero filter coefficients and does not schedule multiplications between input data and corresponding zero filter coefficients. The technique calculates output data from input data and corresponding even filter coefficients, and then calculates output data from input data and corresponding odd filter coefficients.

This method preferably employs a data processor having data cache of a predetermined memory capacity. The image is processes as sub-band stripes having a size selected corresponding to the memory capacity of the data cache. In many cases, the filter coefficients include at least one pair of equal filter coefficients. For these equal pair of filter coefficients, the method adds the corresponding pairs of input data and multiplies the by said corresponding equal filter coefficient, hi the preferred embodiment this addition employs an instruction that simultaneously and separately adds differing portions of packed input data generating separate packed sums. This invention proposes a more efficient multiplier utilization involving no additional buffering in performing inverse discrete wavelet transforms. This invention arranges data flow providing a more efficient use of memory bandwidth and cache space than other known methods. BRIEF DESCRIPTION OF THE DRAWINGS

These and other aspects of this invention are illustrated in the drawings, in which:

FIG. 1 illustrates transformed wavelet data divided into four quadrant subbands (Prior Art);

FIG. 2 illustrates further division of quadrant 201 into smaller bands (Prior Art); FIG. 3 illustrates a third level division of subquadrant 211 into yet smaller bands

(Prior Art);

FIG. 4 illustrates a fourth level division of subquadrant 221 into still smaller bands (Prior Art);

FIG. 5 illustrates scale dependency of wavelet coefficients (Prior Art); FIG. 6 illustrates an example Morton scanning order and a corresponding set of wavelet coefficients used in a coding example (Prior Art);

FIG. 7 illustrates a vertical stage of the inverse wavelet transform (Prior Art);

FIG. 8 illustrates a horizontal stage of the inverse wavelet transform (Prior Art);

FIG. 9 illustrates the algorithm of this invention in flow chart form; FIG. 10 illustrates the operation of an advanced multiply-sum instruction that is useful in this invention (Prior Art);

FIG. 11 illustrates the operation of a dual addition instruction that is useful in this invention (Prior Art);

FIG. 12 illustrates the equal symmetrical pairs of coefficients in a first part of a Vertical Stage of the Antonini Inverse 9-7 JPEG2000 Wavelet Transform; and

FIG. 13 illustrates the equal symmetrical pairs of coefficients in a second part of a Vertical Stage of the Antonini Inverse 9-7 JPEG2000 Wavelet Transform. DETAILED DESCRIPTION OF THE EMBODIMENTS

FIGS. 7 and 8 illustrate the steps in computing the traditional inverse discrete wavelet transform. FIG. 7 illustrates a vertical stage of the inverse wavelet transform. The wavelet encoded data 200 is supplied to re-sampling filters 211 and 212. Their outputs are filtered by respective high pass filter 221 and low pass filter 222. The result from summer 230 is the partially reconstituted picture 240. FIG. 8 illustrates a horizontal stage of the inverse wavelet transform. The partially decoded data 240 is supplied to re-sampling filters 251 and 252. Their outputs are filtered by respective high pass filter 261 and low pass filter 262. The result from summer 270 is the reconstituted picture 280.

FIGS. 7 and 8 show that the inverse wavelet transform can be viewed as the summation of the outputs from filtering the input low pass and high pass subbands by a factor of 2. This type of interpolation often occurs in conjunction with re-sampling filters. The inverse wavelet transform constitutes an example of increasing the sampling rate by a factor of 2 as compared to the sampling rate of the existing low pass and high pass subbands.

Several factors would facilitate computation of the inverse discrete wavelet transform. It would desirable to avoid explicit interpolate of the data by a factor of 2. The process of interpolation by a factor of M in general involves inserting M-I zeros between input samples. It would also desirable to implement the inverse wavelet transform to minimize additional control or storage required due to not implementing explicit interpolation. The basic steps are examined for a data sequence and four-tap filter to illustrate the conventional approach and the previous approaches shown in the literature.

Let X represent the data for the approximation coefficients and Y represent the data for the detail coefficients, with X and Y having the same size. Then:

X = {xo, xi, x2,...xN-i}, Y = {yo, yi, y2, y3,...yN-i}

The interpolated data can be viewed as follows:

Xint = (X0, O5 xi, O, X2, O, X3, 0,...XN-I } Yint = {yo, O, yi, O, y2, O, y3, 0,...yN-1}

The process of interpolation by a Daubechies_4 wavelet can now be viewed as shown in Listing 1 :

Figure imgf000008_0001
where: hi is the corresponding low pass filter coefficient; and gi is the corresponding wavelet filter coefficient.

Kethamawaz and Grisworld (Efficient Wavelet Based Processing, EE Dept. Technical Report 1999) reformulated this computation to be more efficient by formulating the following computation that seeks to maximize the number of multiplies that can be performed by looking across multiple output samples. This is shown in Listing 2.

Figure imgf000010_0001

This reformulation takes an input sample and multiplies it by all the filter taps and stores the partial products. The process is repeated for the next input sample except that its intermediate product is added with the previously computed partial outputs two stages down from the previous iteration as shown in Listing 2. For example X]ho is added with the partial product X0Ii2, which is the partial output 2 stages away from the previous iteration. This algorithm is a depth first search (DFS) solution to the problem in which every input is used to perform computations across the depth of the wavelet. This algorithm requires access to only one input (horizontal) and one line of input (vertical) case.

This algorithm has several implementation implications when used on a real-time digital signal processor (DSP) with cache. The computed intermediate products are the product of two 16-bit numbers and thus need retained in held in memory at 32-bits. Implementing this algorithm in a vectorized manner requires an output buffer to hold two complete interpolated output lines of the size of the width of the image at 32-bits. This requirement greatly increases the needed storage. Since the algorithm calculates results in a staggered manner after several iterations of accumulation, a separate loop may be required to shift the accumulated outputs from 32-bits back to 16-bits depending on the length of the wavelet. This amount of data needing storage and the amount of read and write traffic degrades performance significantly for most cache architectures.

This invention includes data specific transformations to make the algorithm more cache-friendly. This invention uses an alternative constituting a breadth first search algorithm. This invention requires storage of half the number of inputs (horizontal) and half the number of input lines (vertical) at any given time. The input lines do not need to be the same width as the width of the subband. This feature allows processing the subbands as stripes and trades data cache degradation by adjusting the amount of the working data set. The following example is shown for the Daubechies wavelet (D-4), the wavelet of this example, but is applicable to any wavelet family.

This algorithm has several implementation implications when used on a real-time digital signal processor (DSP) with cache. The computed intermediate products are the product of two 16-bit numbers and thus need retained in held in memory at 32-bits. Implementing this algorithm in a vectorized manner requires an output buffer to hold two complete interpolated output lines of the size of the width of the image at 32-bits. This requirement greatly increases the needed storage. Since the algorithm calculates results in a staggered manner after several iterations of accumulation, a separate loop maybe required to shift the accumulated outputs from 32-bits back to 16-bits depending on the length of the wavelet. This amount of data needing storage and the amount of read and write traffic degrades performance significantly for most cache architectures.

This invention includes data specific transformations to make the algorithm more cache-friendly. This invention uses an alternative constituting a breadth first search algorithm. The invention changes the order of computation of the output results. When using dedicated hardware, the prior art trend for increased through put is use an algorithm with a maximum number of multiplies. A special hardware circuit is constructed to simultaneously calculate these multiplications. This results in a process that is efficient from the standpoint of computation hardware. The prior art generally employs such a depth-first computation even when implemented on a programmable digital signal processor. There are disadvantages to this approach. This approach requires storage of large amounts of intermediate data. Often this intermediate data is multiplication products in extra bit form. It is well known that multiplication of two n-bit factors results in a 2n-bit product. The number of such intermediate products are typically greater than can be accommodated in the register set requiring storage in memory or cache. The greater data width of such intermediate products means they cannot be packed into data words as typical for the starting data and filter coefficients. Thus the number and size of intermediate products requires much data transfer traffic between the DSP core and cache increasing the likelihood of cache missing and resort to slower main memory.

The breadth first technique of this invention focuses on immediate calculation of output results. This means that fewer results are computed in parallel but that the latency of each final result is smaller. This smaller latency produces less intermediate data that needs to be stored. Often, as shown below, this smaller amount of immediate data can be completely stored in the register set. This permits rounding and packing the output data in the same manner as the input data is packed. This reduces memory traffic by requiring storage of only output data in packed form. Accordingly, this invention is less likely to generate cache misses that require time consuming memory accesses.

This invention requires storage of half the number of inputs (horizontal) and half the number of input lines (vertical) at any given time. The input lines do not need to be the same width as the width of the subband. This feature allows processing the subbands as stripes and trades data cache degradation by adjusting the amount of the working data set. The following example is shown for the Daubechies wavelet (D-4), the wavelet of this example, but is applicable to any wavelet family.

This invention partitions the existing low pass filter H and the wavelet filter G into two sub-filters each as shown:

H = (H0, H1, H2, H3) ≡ Hev = (H0, H2}, Hod = {Hi, H3) G = {Go, Gi, G2, G3) ≡ Gev = {Go, G2), God = {Gh G3)

Thus this invention uses a poly-phase scheme to implement the transformation. The first output sample is produced by taking half the number of inputs from each subband and convolving with the respective even filters. The second output sample uses the same inputs, but convolves with the odd filter. This algorithm produces outputs required for steady state only, so only output samples starting at output sample 2 are valid. This algorithm can be extended to produce the filtered output samples at startup by pre-appending the sequence with zeroes. This modification is not required since the inverse transform only deals with steady state outputs. This is shown in Listing 3.

The above output computation can now be simplified by elimination of the zero data terms as shown in Listing 4.

Figure imgf000014_0001

This can be re-written using the "*" symbol to denote the convolution operator and the poly-phase filters Hev, Hod, Gev and God as shown in Listing 5 to represent the underlying math.

Output Sample 2: {x0, xj } *Hev + {y0, y, } *Gev

Output Sample 3 : {x0, X1 } *Hod + {y0, yi } *God

Output Sample 4: {xj, x2} *Hev + Jy1, y2} *Gev

Output Sample 5: {xls x2} *Hod + {yu y2} *God

Listing 5

An analysis of these equations shows that the number of multiplies required to produce an output sample are half the length of the associated scaling and wavelet filter. The elimination of the useless multiplies with the inserted zero data produces this result. Thus, the current method achieves efficient utilization of the inputs with the concurrent evaluation of the outputs in a non-staggered (one-shot) manner. In addition the input advances by one input (horizontal) or one input line (vertical) at a time. This permits some data re-use. This also removes the need for temporary buffering of the accumulated 32-bit precision results. The non-staggered method of computation produces complete results at the end of every iteration. The accumulated results can be held in the register file prior to shifting them down to 16-bits. This eliminates the need for additional storage of 32-bit products, resulting in a substantial improvement in cache performance.

FIG. 9 illustrates this algorithm in flow chart form. The process initiates in block 901. The next data is recalled in block 902. The even coefficients in h and g are recalled in block 903. The convolution and sum as shown in Listings 4 and 5 is performed in block 904. As shown in Listing 4 this convolution-sum calculation requires four multiplies and three sums. The manner of performing this operation will be described further below. The odd coefficients in h and g are recalled in block 905. Block 906 performs the convolution-sum using the odd coefficients. Decision block 907 determines if the computation has reached the end of the image. If so (Yes at decision block 907), then the inverse discrete wavelet transform is complete and the process ends at end block 908. If not (No at decision block 907), process control passes back to block 902 to recall the next x and y data points. The manner of computing the convolution-sums of blocks 904 and 906 depends upon the instruction set of the digital signal processor doing the computation. Several examples will be described below. In these examples assume that the x and y data is 16-bits with adjacent x and adjacent y data points packed into a 32-bit data word. The h and g coefficients are also assumed to be 16-bits each packed into 32-bit data words. Assume that coefficients ho and h2, and g0 and g2 are packed into corresponding even data words and coefficients h1 and h3, and g1 and g3 are packed into corresponding odd data words.

In the simplest instruction set case, no complex multiply-sum instructions are available. The Texas Instruments TMS320C6000 family of digital signal processors includes four 16-bit by 16-bit multiply instructions. These four instructions employ a selected 16-bit half of each of two input operands and form a 32-bit product. The four instructions include: 1) least signification 16-bits of the first operand times the least significant 16-bits of the second operand; 2) least signification 16-bits of the first operand times the most significant 16-bits of the second operand; 3) most signification 16-bits of the first operand times the least significant 16-bits of the second operand; and 4) most signification 16-bits of the first operand times the most significant 16-bits of the second operand. The four multiplications are performed using appropriate ones of these four multiply instructions. Each of these products is stored in a 32-bit general purpose register. The four products are summed in three add operations. Care must be taken to avoid or recover from any overflow during these sums. The final sum is right then shifted 16 bits to retain the most significant bits. This final sum is packed into a 32-bit data register with another final sum for output as two packed 16-bit output data points. Some digital signal processors include pack instructions which can do this final packing operation in one instruction cycle. If no such instruction is available, then the packing can take place as follows. One of the 16-bit data words is left shifted 16-bits into the 16 most significant bits of a data word. This left shifted data word and the other data word (with the 16 most significant bits set as zeros) are logically ANDed. This effectively packs the two output data points in the most significant and least significant halves of the 32- bit data word.

This method employs 8 or 9 instructions not counting dealing with sum overflow. These are four multiplies, three adds and a single pack operation or a shift operation followed by a logical AND operation. The TMS320C6000 includes two separate data paths, each with a multiply unit. Thus the x computation could be done by the A data path and the y computation could be done by the B data path. Since two multiplies and four arithmetic/logical operations can be dispatched each instruction cycle, a pipelined software loop may be able to compute one output data point for each two instruction cycles. The number of general purpose data registers required by an inner loop should be well within the 16 general purpose data registers for each data path.

Other digital signal processors have advanced instructions permitting faster operation. The Texas Instruments TMS320C6400 digital signal processor includes an instruction called Dot Product With Shift and Round, Signed by Unsigned Packed 16-Bit (DOTPRSU2). The operation of this instruction is illustrated in FIG. 10. The most significant and least significant halves of the first source operand (Source 1) are multiplied by the corresponding most significant halves of the second source operand (Source X). The two products are added and rounded by addition of a "I1 bit in the 15 bit position (hex8000). Then sum is right shifted by 16 bits. The most significant bits are sign extended, that is filled with l's if the most significant bit is 1 and filled with O's if the most significant bit is 0. This instruction performs half the convolution-sum computation required. This instruction requires the first operand to be a signed number. The x and y data would typically be pixels and therefore unsigned data. However, the h and g coefficients could easily be signed numbers. The multiply-sum quantities for the x data and the y data would be added and then right shifted and packed into 16-bits.

The DOTPRSU2 is a multi-cycle instruction requiring for instruction cycles to complete. However, it can be fully pipelined with a new instruction started every cycle. Because each multiply unit in the TMS320C6400 can complete one convolution each cycle, this data signal processor can be software pipelined to produce one output data point per instruction cycle.

Analysis reveals that this technique will work for odd-tap filter wavelets like the JPEG2000 9-7 wavelet. The first output sample is the one where the input is convolved with the odd poly-phase filter, and the second output sample is the one convolved with the even poly phase filter. This ordering of coefficients is exactly opposite to the order for even tap wavelets (D-4) where the first output sample is the one where the input is convolved with the even poly-phase filter and the second output sample is the one where the input is convolved with the odd poly-phase filter. The use of this technique can be shown to be general and applicable for other interpolation factors that often arise in the context of re-sampling.

Listings 6 and 7 are example loop kernels incorporating this invention in scheduled assembly code for the Texas Instruments TMS320C62x digital signal processor for the Vertical Stage of the Antonini Inverse 9-7 JPEG2000 Wavelet Transform. Listing 6 is for the 7 data term case and listing 7 is for the 9 data term case. Note that each listing computes two separate output samples per iteration.

Figure imgf000018_0001

Figure imgf000019_0001

Figure imgf000019_0002
Figure imgf000020_0001
Figure imgf000021_0001

LISTING 7

The TMS320C62X is an 8-way very long instruction word (VLIW) processor employing two separate datapaths/register files with limited cross path capability. The eight execution units include data units Dl and D2, logic units Ll and L2, multiply units Ml and Ml and arithmetic units Sl and S2. The "1" units use the A register file and the "2" units use the B register file with limited cross connection. The instruction mnemonics are interpreted as follows. The initial symbol " | " indicates this instruction executes in parallel with the immediately prior instruction. AU instructions can be made conditional. This listing only notes those instructions that are conditional. The next symbol in square brackets "[]" indicates the predicate register and the "!" symbol indicated inverted logic. Next is the instruction type mnemonic, then the execution unit (D, L, M and S). An "X" suffix on the execution unit designator indicates the use of a crosspath transferring an operand from a register file opposite that of the execution unit. The operand list is last. This operand list is source 1, source2, destination except for loads and stores and immediate operands.

The process is initialized via instructions not shown in these loop kernels. Only one of the loops of Listings 6 and 7 is used dependent upon whether this loop is operating on 7 input samples (Listing 6) or nine input samples (Listing 7) in the 9-7 scheme. This initialization includes storing the following constants. In Listing 6: hi in register A5; h3 in register BO; gl in register B2; g3 in register A6; the initial data addresses in registers AO, AlO, Al l, A12, B7, B8 and B9. In Listing 7: h0 in register Bl 1 ; h2 in register A4; g0 in register A6; g2 in register B12; g4 in register AlO; the initial data addresses in registers Al l, A12, A13, A14, A15, Bl, B2, B3 and BlO. Listing 6 has five execute packets (A, B, C, D and E) having five to eight numbered instructions. Listing 7 has six execute packets (F, G, H, I, J and I) having five to seven numbered instructions. Following initialization the proper loop is selected. The loop kernels include load 32-bit word instructions for loading the data. This data is 16-bit data with two data points packed in each 32-bit word. Listing 6 includes 7 LDW instructions at A3, A4, B6, B7, C5, C6 and E6. Listing 7 includes 9 LDW instructions at F4, F5, G2, G4, H3, H4, 12, J5 and 36. These load instructions use an indirect address stored in respective source registers (B9, AO, AlO, B7, Al l, B8 and A12 for Listing 6 and Al l, BlO, A12, B3, A15, Bl, Al 3, B2 and A14 for Listing 7). The "*R++" address mode post increments the address register. The address stored in the address register increments to point to the next 32-bit data word following the data fetch. The sets the address registers to point to and fetch sequential 32-bit words in memory. The corresponding destination registers thus store two packed 16- bit data words.

Each loop kernel adds some data using the ADD2 instruction (Listing 6 at instructions A2, B4 and C2 and Listing 7 at instructions F6, G3, H2 and K4.) Operation of the ADD2 instruction is illustrated in FIG. 11. The ADD2 instruction separately adds the 16-bit upper and lower parts of the source operands. If the lower sum B + D overflows and generates a carry, that carry does not propagate to the upper sum A + C. This uses the distributive property of multiplication over addition and a property of the coefficients in the Antonini Inverse 9-7 JPEG2000 Wavelet Transform. These coefficients are equal in symmetrical pairs. For example, in Listing 6 the coefficients h0 and h2 are equal and the coefficients g0 and g3 are equal and the coefficients gl and g2 are equal. The filter calculation includes:

prodl = xoh0 + x2h2 Let h0 = h2 = h, so: prodl = xoh + x2h prodl = h(xo + x2)

The original calculation requires two multiplications and one addition. The algebraically transformed calculation requires one addition and one multiplication. Performing the data addition before the multiplication permits the elimination of one multiplication. Since multiplication hardware is at a premium on the TMS320C62x digital signal processor relative to addition hardware, this serves to reduce number of instructions for each iteration of the loop.

Next comes the filter multiplications. These multiplications are at instructions B3, B5, C3, C4, D3, D5, E2 and E5 in Listing 6 and instructions F3, G5, G6, H5, H6, 13, 15, J2, J4 and K3 in listing 7. The algebraic transformation eliminated 6 of 14 multiplies in Listing 6 and 8 of 18 multiplies in Listing 7. These multiplications include a mixture of MPY, MPYH, MPYLH and MPYHL as required to select the proper pair of 16-bit operands in the packed 32-bit words. The 32-bit product results are stored in respective data registers as 32- bit quantities.

The next process is the addition for accumulation. These additions take place separately for the two output samples, hi Listing 6 the accumulations are at instructions AO, Al, BO, Bl, B2, CO and E4. In Listing 7 the accumulations are at instructions Fl, F2, GO, Gl, HO, 14, Jl, J3, Kl and K4.

Each listing then rounds the accumulations by right shifting 15 bits. These instructions are instructions Cl and DO in Listing 6 and instructions Hl and 10 in Listing 7.

The results are stored via STH instructions Dl and El in Listing 6 and instructions Il and KO in listing 7. The STH instruction stores the 16 least significant bits of the source register into the memory address specified in the address register. Instructions Dl and Il use the addressing form "*R++(4)" which specifies auto increment of the address register after the store operation by 4 bytes (32 bits). Instructions El and KO use the addressing form "*- R(2)" which provides a negative offset of 2 bytes before the store operation. This effectively packs two 16-bit results into a single 32-bit data word.

Listing 6 includes a prolog collapse technique. The two store instructions STH Dl and El are made conditional upon register A2 being zero. Register A2 is decremented in subtract instruction EO whenever the register is non-zero. The main steady-state loop can serve the prolog initialization function of filling the intermediate results before steady state is reached. During these prolog iterations of the main loop the results are invalid. Making the store instructions conditional prevents these initial, invalid results from being stored. The register A2 is initiated with an integer equal to the number of such prolog iterations of the loop. In this example register A2 is initialized as "2." The subtract instruction EO decrements register A2 for this number of loops. During the prolog loops the conditional store instructions Dl and El are not executed. When register A2 is decremented to zero, subtract instruction EO is disabled and the store instructions Dl and El are enabled.

Listing 6 and 7 include loop control instructions. The loop control instructions in Listing 6 are instructions D2, EO and E5. Instruction D2 effectively decrements the loop variable stored in register Al. Instruction E5 branches to the beginning of the loop (LOOP) when register Al is non-zero. When register Al has decremented to zero, instruction E5 exits the loop to the next instruction, which is outside the loop kernel of Listing 6. Instruction EO decrements the contents of register A2. Making the STH instructions Dl and El conditional on A2 being zero ([!A2]). The loop control instructions in Listing 7 are instructions JO and FO. Instruction JO effectively decrements the loop variable stored in register BO. Instruction FO branches to the beginning of the loop (LOOPH) when register Bo is non-zero. When register BO has decremented to zero, instruction FO exits the loop to the next instruction, which is outside the loop kernel of Listing 7.

Table 1 lists a summary of the loop kernels of Listing 6 and 7.

Figure imgf000024_0001
Listings 6 and 7 are merely examples of implementation of this algorithm on a particular data processor architecture. One skilled in the art would realize that the method of this invention is readily practiced on other data processor architectures. One essential feature is storing only packed data in the memory only. The memory stores the input data as packed data. The filter results are likewise stored as packed data into the memory. This feature reduces the memory size needed and reduces the number of required memory accesses.

This invention is also applicable to other filtering processes than the inverse Wavelet Transform. For example, conversion of an image to a different size could use this technique. In up-sampling (conversion to a larger image size) the input coefficients in the filter function will often be zero, hi accordance with this invention, these multiplications are omitted from performance. The breadth-first computation technique is also useful in this application. Other applications include audio time- warping which changes the run-time of an audio file without changing its pitch.

Claims

CLAIMSWhat we claim is:
1. A method for inverse Wavelet transform comprising the steps of: storing input data as subwords packed within data words in a memory; calculating output data in a software loop in a breadth-first fashion by recalling input data and filter coefficients for computation of at least one output data for each iteration of the software loop even if the same input data is used in a later iteration for calculating other output data; storing calculated output data as subwords packed within data words in the memory.
2. The method of Claim 1, wherein: said step of calculating output data in a software loop includes scheduling multiplications between input data and corresponding non-zero filter coefficients and not scheduling multiplications between input data and corresponding zero filter coefficients.
3. The method of Claim 1, wherein: said step of calculating output data in a software loop includes calculating output data from input data and corresponding even filter coefficients, and calculating output data from input data and corresponding odd filter coefficients.
4. The method of Claim 1 , wherein: said step of calculating output data in a software loop includes for each level of wavelet encoding calculating an expansion of vertical spatial frequency, and thereafter calculating an expansion of horizontal spatial frequency.
5. The method of Claim 1, wherein: said method is preformed upon a data processor having data cache of a predetermined memory capacity; said method further comprising the step of processing an image as sub-band stripes having a size selected corresponding to said predetermined memory capacity of said data cache.
6. The method of Claim 1, wherein: said filter coefficients include at least one pair of equal filter coefficients; said method further comprising the steps of: adding pairs of input data to be multiplied by a corresponding pair of equal filter coefficients; and multiplying a sum of said pairs of input data by said corresponding equal filter coefficient.
7. The method of Claim 6, wherein: said step of adding pairs of input data includes simultaneously separately adding differing portions of input data packed into a single memory word generating separate sums packed into a single memory word.
PCT/US2006/031878 2006-08-15 2006-08-15 Cache friendly method for performing inverse discrete wavelet transform WO2008020839A3 (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/US2006/031878 WO2008020839A3 (en) 2006-08-15 2006-08-15 Cache friendly method for performing inverse discrete wavelet transform

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/US2006/031878 WO2008020839A3 (en) 2006-08-15 2006-08-15 Cache friendly method for performing inverse discrete wavelet transform

Publications (2)

Publication Number Publication Date
WO2008020839A2 true true WO2008020839A2 (en) 2008-02-21
WO2008020839A3 true WO2008020839A3 (en) 2008-06-26

Family

ID=39082476

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/US2006/031878 WO2008020839A3 (en) 2006-08-15 2006-08-15 Cache friendly method for performing inverse discrete wavelet transform

Country Status (1)

Country Link
WO (1) WO2008020839A3 (en)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6567081B1 (en) * 2000-01-21 2003-05-20 Microsoft Corporation Methods and arrangements for compressing image-based rendering (IBR) data using alignment and 3D wavelet transform techniques
US6996287B1 (en) * 2001-04-20 2006-02-07 Adobe Systems, Inc. Method and apparatus for texture cloning
US7149362B2 (en) * 2001-09-21 2006-12-12 Interuniversitair Microelektronica Centrum (Imec) Vzw 2D FIFO device and method for use in block based coding applications
US7236637B2 (en) * 1999-11-24 2007-06-26 Ge Medical Systems Information Technologies, Inc. Method and apparatus for transmission and display of a compressed digitized image
US7283684B1 (en) * 2003-05-20 2007-10-16 Sandia Corporation Spectral compression algorithms for the analysis of very large multivariate images

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7236637B2 (en) * 1999-11-24 2007-06-26 Ge Medical Systems Information Technologies, Inc. Method and apparatus for transmission and display of a compressed digitized image
US6567081B1 (en) * 2000-01-21 2003-05-20 Microsoft Corporation Methods and arrangements for compressing image-based rendering (IBR) data using alignment and 3D wavelet transform techniques
US6996287B1 (en) * 2001-04-20 2006-02-07 Adobe Systems, Inc. Method and apparatus for texture cloning
US7149362B2 (en) * 2001-09-21 2006-12-12 Interuniversitair Microelektronica Centrum (Imec) Vzw 2D FIFO device and method for use in block based coding applications
US7283684B1 (en) * 2003-05-20 2007-10-16 Sandia Corporation Spectral compression algorithms for the analysis of very large multivariate images

Also Published As

Publication number Publication date Type
WO2008020839A3 (en) 2008-06-26 application

Similar Documents

Publication Publication Date Title
Moreland et al. The FFT on a GPU
Nuzman et al. Auto-vectorization of interleaved data for SIMD
US5960012A (en) Checksum determination using parallel computations on multiple packed data elements
Andra et al. A VLSI architecture for lifting-based forward and inverse wavelet transform
Daubechies et al. Factoring wavelet transforms into lifting steps
US6871208B1 (en) Parallel adder-based DCT/IDCT design using cyclic convolution
US4849922A (en) Circuit for computing the quantized coefficient discrete cosine transform of digital signal samples
Lee Subword parallelism with MAX-2
Mulders et al. On lattice reduction for polynomial matrices
US6219458B1 (en) Overlapped reversible transforms for unified lossless/lossy compression
US5880979A (en) System for providing the absolute difference of unsigned values
Wu et al. A high-performance and memory-efficient pipeline architecture for the 5/3 and 9/7 discrete wavelet transform of JPEG2000 codec
US5493524A (en) Three input arithmetic logic unit employing carry propagate logic
US5596763A (en) Three input arithmetic logic unit forming mixed arithmetic and boolean combinations
US5590350A (en) Three input arithmetic logic unit with mask generator
US7430578B2 (en) Method and apparatus for performing multiply-add operations on packed byte data
US5936872A (en) Method and apparatus for storing complex numbers to allow for efficient complex multiplication operations and performing such complex multiplication operations
US6067613A (en) Rotation register for orthogonal data transformation
US5420809A (en) Method of operating a data processing apparatus to compute correlation
Aksoy et al. Search algorithms for the multiple constant multiplications problem: Exact and approximate
US5983253A (en) Computer system for performing complex digital filters
US5634065A (en) Three input arithmetic logic unit with controllable shifter and mask generator
US6098163A (en) Three input arithmetic logic unit with shifter
US6237016B1 (en) Method and apparatus for multiplying and accumulating data samples and complex coefficients
US5838377A (en) Video compressed circuit using recursive wavelet filtering

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 06801555

Country of ref document: EP

Kind code of ref document: A2

NENP Non-entry into the national phase in:

Ref country code: DE

NENP Non-entry into the national phase in:

Ref country code: RU

122 Ep: pct application non-entry in european phase

Ref document number: 06801555

Country of ref document: EP

Kind code of ref document: A2