Hybrid architecture for multiple transforms for signal processing applications

Abstract

A hybrid architecture for transforms such as N-point Discrete Fourier Transform (DFT), Discrete Cosine Transform (DCT), Discrete Sine Transform (DST) and Discrete Wavelet Transform (DWT) has been proposed and implemented using triple matrix product method. There is limited work reported on efficient single architecture that can perform multiple transforms simultaneously or serially depending on the application requirement. The Hybrid architecture implemented, can compute various transforms efficiently. A controller is designed which can perform different transforms using the hybrid architecture based on the input provided. The implemented systolic array can be used for computing the diagonal elements of triple-matrix product. The designed architecture produces the output of transform sequence in order, which avoids reordering at output. The implemented architecture can be used to handle large sized transforms by repeatedly using fixed size architecture for a large number of points without increasing the number of Processing Elements (PEs). The proposed architecture has been validated with a watermarking algorithm that uses DCT and DWT transforms and its performance analyzed. The proposed hybrid architecture is implemented on Spartan-7 xc7s100fgga676-1. The simulation results are given and analyzed against standalone architecture.

Keywords

Hybrid architecture transforms implementation triple matrix product method

1 Introduction

Signal Processing Transforms are used to transform a signal from one domain to another domain. Some of the widely used transforms in signal processing are Discrete Fourier Transform (DFT), Discrete Sine Transform (DST), Discrete Cosine Transform (DCT), Discrete Wavelet Transform (DWT), Fractional Fourier Transform (Fr FT), Discrete Hartley Transform (DHT), Hilbert Transform, Hadamard Transform, Hough Transform, Discrete Walsh Transform etc...

Systolic array is a form of pipe-lining, sometimes in more than one dimension. Systolic computing is a set of simple processing elements with regular and local connections which takes external inputs and processes them in a predetermined manner in a pipe-lined fashion. The first systolic structure was developed by Kung [1] for DFT. After that, many designs using systolic architectures were proposed. An early paper by Aravena explains the parallel representation of linear algorithms, which converts the matrix-vector operations to sum of small size triple-matrix products [2]. In the paper [4], a novel hybrid architecture to compute 8-point 1-D transforms such as DCT, DFT and HWT is proposed. The hardware is implemented based on Matrix decomposition algorithm i.e, Matrix Decomposition of the 8x8 DCT, DFT and HWT. The transforms are proposed for 1-D, 8x8 and forward. In the paper [3], a unified VLSI implementation of the DCT /DST /IDCT/IDST is proposed. The hardware implementation is done based on cyclic convolution structures of transforms. The transforms are proposed only for DCT /DST /IDCT/IDST.

Though there are unifies architectures reported for few multiple transforms, there is no efficient single architecture that can perform multiple transforms simultaneously or serially depending on the application requirement. The objective is to design and implement a scalable architecture for multiple transforms like DFT, DCT, DST, DHT and DWT.

In [4], 2D systolic implementation of DFT is proposed. In paper [5], the DFT algorithm proposed in [4] is rewritten to form a Triple-Matrix-Product form to map it to 2D systolic array. In this paper, the same is extended to multiple transforms such as DCT, DST and DWT and applied to a watermarking application.

Rest of this report is organized as follows. Section 2 describes the matrix-based representation of the algorithms. Hybrid multi-transform architecture that can compute different transforms is described in Section 3. Section 4 contains Simulation results and Conclusions are given in Section 5.

2 Matrix based representation of transforms

The conversion of DFT algorithm into a triple-matrix product representation and mapping on to a 2D systolic array structure is explained in [5]. In this work, the concept is extended to DCT, DST and DWT.

2.1 Triple matrix product form of 1-D DCT

Discrete Cosine Transform (DCT) of an N point sequence x(n) is given by,

$\begin{matrix} X (k) = C (k) \sum_{n = 0}^{N - 1} x (n) cos [\frac{(2 n + 1) k π}{2 N}] \\ for k = 0, 1, 2 . . . N - 1 \end{matrix}$ (1) $C (k) = {\begin{matrix} \frac{1}{\sqrt{2}} & for k = 0 \\ 1 & for 1 \leq k \leq N - 1 \end{matrix}$ Decomposing N = N₁N₂ and re-indexing n = N₁n₁ + n₂. Then, $\begin{matrix} X (k) = \sum_{n_{1} = 0}^{N_{2} - 1} \sum_{n_{2} = 0}^{N_{1} - 1} x (N_{1} n_{1} + n_{2}) * \\ cos [\frac{(2 (N_{1} n_{1} + n_{2}) + 1) k π}{2 N}] \end{matrix}$ (2) $\begin{matrix} X (k) = \sum_{n_{1} = 0}^{N_{2} - 1} \sum_{n_{2} = 0}^{N_{1} - 1} x (N_{1} n_{1} + n_{2}) * \\ cos [\frac{2 N_{1} n_{1} k π}{2 N} + \frac{(2 n_{2} + 1) k π}{2 N}] \end{matrix}$ (3) Equation (3) can be written in the matrix form as $X (k) = A_{{kn}_{1}} * X_{n_{1} n_{1}} * A_{n_{2} k} - B_{{kn}_{1}} * X_{n_{1} n_{1}} * B_{n_{2} k}$ (4) where, $X_{n_{1} n_{2}} = x (N_{1} n_{1} + n_{2}) isa N_{2} * N_{1} matrix,$ $A_{{kn}_{1}} = cos (\frac{2 N_{1} n_{1} k π}{2 N}) isaN * N_{2} matrix,$ $A_{n_{2} k} = cos (\frac{(2 n_{2} + 1) k π}{2 N}) isa N_{1} * Nmatrix,$ $B_{{kn}_{1}} = sin (\frac{2 N_{1} n_{1} k π}{2 N}) isaN * N_{2} matrix,$ $B_{n_{2} k} = sin (\frac{(2 n_{2} + 1) k π}{2 N}) isa N_{1} * Nmatrix .$ DCT given by (4) is the sum of two matrices diagonal element of the triple matrix product. Each of the triple matrices are computed and the outputs of the two triple matrix products are subtracted in order to obtain the DCT output.

2.2 Triple matrix product form of 1-D DST

Discrete Sine Transform (DST) of an N point sequence x(n) is given by,

$\begin{matrix} X (k) = \frac{2}{N} e (k) \sum_{n = 0}^{N - 1} x (n) sin [\frac{(2 n + 1) k π}{2 N}] \\ for k = 0, 1, 2 . . . N - 1 \end{matrix}$ (5) $e (k) = {\begin{matrix} \frac{1}{\sqrt{2}} & for k = 0 \\ 1 & for 1 \leq k \leq N - 1 \end{matrix}$ Decomposing N = N₁N₂ and re-indexing n = N₁n₁ + n₂. Then, $\begin{matrix} X (k) = \sum_{n_{1} = 0}^{N_{2} - 1} \sum_{n_{2} = 0}^{N_{1} - 1} x (N_{1} n_{1} + n_{2}) * \\ sin [\frac{(2 (N_{1} n_{1} + n_{2}) + 1) k π}{2 N}] \end{matrix}$ (6) $\begin{matrix} X (k) = \sum_{n_{1} = 0}^{N_{2} - 1} \sum_{n_{2} = 0}^{N_{1} - 1} x (N_{1} n_{1} + n_{2}) * \\ sin [\frac{2 N_{1} n_{1} k π}{2 N} + \frac{(2 n_{2} + 1) k π}{2 N}] \end{matrix}$ (7) Equation (7) can be written in the matrix form as $X (k) = A_{{kn}_{1}} * X_{n_{1} n_{1}} * A_{n_{2} k} + B_{{kn}_{1}} * X_{n_{1} n_{1}} * B_{n_{2} k}$ (8) where, $X_{n_{1} n_{2}} = x (N_{1} n_{1} + n_{2}) isa N_{2} * N_{1} matrix,$ $A_{{kn}_{1}} = sin (\frac{2 N_{1} n_{1} k π}{2 N}) isaN * N_{2} matrix,$ $A_{n_{2} k} = cos (\frac{(2 n_{2} + 1) k π}{2 N}) isa N_{1} * Nmatrix,$ $B_{{kn}_{1}} = cos (\frac{2 N_{1} n_{1} k π}{2 N}) isaN * N_{2} matrix,$ $B_{n_{2} k} = sin (\frac{(2 n_{2} + 1) k π}{2 N}) isa N_{1} * Nmatrix .$ DST given by (8) is the sum of two matrices diagonal element of the triple matrix product. Each of the triple matrices are computed and the outputs of the two triple matrix products are added in order to obtain the DST output.

2.3 Matrix-matrix product form of 1-D DWT

Based on the size of the input sequence, it is possible to represent the wavelet coefficients in the form of a NxN matrix. In such cases, the output of the wavelet transform will be the matrix multiplication of a NxN matrix with a Nx1 column vector resulting in a Nx1 column vector. The matrix form can be expressed as below: $[\begin{matrix} Y_{0} \\ Y_{1} \\ Y_{2} \\ Y_{3} \end{matrix}] = [\begin{matrix} h_{0} & h_{1} & h_{2} & h_{3} \\ h_{3} & - h_{2} & h_{1} & - h_{0} \\ h_{2} & h_{3} & h_{0} & h_{1} \\ h_{1} & - h_{0} & h_{3} & - h_{2} \end{matrix}] [\begin{matrix} x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \end{matrix}]$ (9) If Haar wavelet is used, then the 4x4 Haar matrix is represented as $[\begin{matrix} Y_{0} \\ Y_{1} \\ Y_{2} \\ Y_{3} \end{matrix}] = [\begin{matrix} h_{0} & h_{1} & 0 & 0 \\ 0 & 0 & h_{0} & h_{1} \\ h_{1} & - h_{0} & 0 & 0 \\ 0 & 0 & h_{1} & - h_{0} \end{matrix}] [\begin{matrix} x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \end{matrix}]$ (10)

2.4 Triple matrix product form of 2-D DWT

The 2-D DWT performs the transformation of a 2-D input signal of size NxN into a transformed domain. This is carried out in two stages, first, by multiplying NxN input matrix with a NxN DWT coefficient matrix and in the second stage the result is multiplied with the transpose of the coefficient matrix. Essentially, it gives a triple matrix product form as shown below: $Y_{NXN} = C_{NXN} * X_{NXN} * C_{NXN}^{T}$ where C_NXN is the DWT coefficient matrix. The 2-D DWT using DB4 wavelet in triple matrix form is as below:

$\begin{matrix} [\begin{matrix} Y_{0} \\ Y_{1} \\ Y_{2} \\ Y_{3} \end{matrix}] = [\begin{matrix} h_{0} & h_{1} & h_{2} & h_{3} \\ h_{3} & - h_{2} & h_{1} & - h_{0} \\ h_{2} & h_{3} & h_{0} & h_{1} \\ h_{1} & - h_{0} & h_{3} & - h_{2} \end{matrix}] * \\ [\begin{matrix} x_{11} & x_{12} & x_{13} & x_{14} \\ x_{21} & x_{22} & x_{23} & x_{24} \\ x_{31} & x_{32} & x_{33} & x_{34} \\ x_{41} & x_{42} & x_{43} & x_{44} \end{matrix}] [\begin{matrix} h_{0} & h_{3} & h_{2} & h_{1} \\ h_{1} & - h_{2} & h_{3} & - h_{0} \\ h_{2} & h_{1} & h_{0} & h_{3} \\ h_{3} & - h_{0} & h_{1} & - h_{2} \end{matrix}] \end{matrix}$ (11)

3 Proposed hybrid multi-transform architecture

The proposed hybrid architecture has two types of PEs A and B as in [5]. The systolic arrays A1X and A2X has PEs A and the systolic array A1X_B1, A1X_B2, A2X_B1 and A2X_B2 has PEs B as shown in Fig. 1. PEs A are used to compute double matrix product. The results of PEs A are multiplied with PEs B to obtain triple matrix product.

Fig. 1

Proposed hybrid architecture.

Let the triple-matrix-product be C=A.X.B. The matrices A and B are co-efficient matrices and X is input matrix. The input data matrix X has real elements, A and B matrix will have both real and complex values as elements. Hence, $DxT = A * X * B$ (12) where, $A = A_{real} + j * A_{imag}$ $A = A_{real} + j * A_{imag}$ $\begin{matrix} DxT = A * X * B = (A_{real} + j * A_{imag}) * X * \\ (B_{real} + j * B_{imag}) \end{matrix}$ (13) $\begin{matrix} = (A_{real} * X * B_{real} - A_{imag} * X * B_{imag}) + \\ j * (A_{real} * X * B_{imag} + A_{imag} * X * B_{real}) \end{matrix}$ (14)

All the inputs are sent in a sequential manner. A_real, A_imag, B_real and B_imag are sent to arrays A1_COFF, A2_COFF, B1_COFF and B2_COFF respectively as shown in Fig. 1. Output of the systolic arrays A1X_B1, A1X_B2, A2X_B1 and A2X_B2 are A_real * X * B_real, A_real * X * B_imag, A_imag * X * B_real and A_imag * X * B_imag respectively. Outputs of systolic arrays A1X_B2 and A2X_B1 are added to get either imaginary part of DFT or DST output based on the input given to the controller. Similarly, outputs of systolic arrays A1X_B1 and A2X_B2 are subtracted to get either real part of DFT or DCT output based on the input given to the controller. The DWT outputs can be directly obtained at the output of systolic arrays A1X_B2 and A2X_B1.

3.1 Mapping of 1-D DCT

As per (4), the DCT computation requires two sets of A.X.B arrays and the result needs to be subtracted. Hence it is clear that the real part of the DxT gives the DCT output if the DCT coefficients are fed. $DCT = (A_{real} * X * B_{real} - A_{imag} * X * B_{imag})$ (15)

3.2 Mapping of 1-D DST

As per (8), the DST computation also requires two sets of A.X.B arrays and adding the outputs of individual arrays, it may be observed that the output is similar to that of the imaginary part of DxT. $DST = (A_{real} * X * B_{imag} + A_{imag} * X * B_{real})$ (16)

3.3 Mapping of 1-D DWT

As per (4), the 1-D DWT computation requires only one set of A.X. array for computation by feeding. Hence the other set of AX array available in Unified architecture can be used for computing another set of 1-D DWT or it can also be used to compute higher order DWT.

3.4 Mapping of 2-D DWT

The 2-D DWT requires the complete computation of triple matrix product unlike the computation of diagonal elements as in case of sinusoidal transforms. Hence the output of the first array is fed as the input to the second A.X array.

A 2x1 multiplexer chooses the inputs either from the external channel (select=0, and is the case for DFT, DCT, DST, 1-D DWT) or from the output of the first array (select=1 for 2-D DWT) based on the choice of the transform.

4 Simulation results

The architecture is implemented for N=16, which is square of the integer 4. The systolic arrays A1X and A2X has 4*4 array of PEs A and the systolic arrays A1X_B1, A1X_B2, A2X_B1 and A2X_B2 has 4*1 array of PEs B.

4.1 Hardware implementation

The hybrid structure and standalone structures are implemented on FPGA device using vivado. They are shown in Figs. 1, 2, and 3 respectively. The standalone architecture can calculate DCT and DWT simultaneously. Table 1 shows the hardware resource utilization report generated through vivado during Run Implementation stage. The target device used is Spartan-7 xc7s100fgga676-1. The application requires calculation of DCT and DWT, which is done serially on hybrid architecture. The hybrid architecture active usage is equivalent to Isolated DCT architecture, because it doesn’t utilize the systolic arrays A1X_B2 and A2X_B2. Therefore, hybrid active used resources are less than combined standalone structures, which can be inferred from Table 2.

Fig. 2

Isolated DCT architecture.

Fig. 3

Isolated DWT architecture.

Table 1

Resource utilization

Resource	Isolated DCT	Isolated DWT	Hybrid
LUT	626	178	834
FF	1211	245	1512
DSP	40	4	48
IO	344	104	362

Table 2

Resource utilization

Resource	Isolated DCT+DWT	Hybrid active usage	% Area reduction
LUT	804	626	22.14
FF	1456	1211	16.83
DSP	44	40	9.09
IO	448	344	23.21

4.2 Time analysis

The watermarking algorithm proposed in [6] is implemented for N=64 i.e, each frame will have 64 samples of data. Two level DWT is performed on the frame in which 16 values remain. These values of two frames are tabulated in Table 5. Mean of the 16 DCT output values for two adjacent frames is calculated which is used to obtain watermark key by comparing with adjacent frame values. The algorithm is implemented on proposed Hybrid architecture and a Standalone architecture which can compute DWT and DCT at the same time. The results obtained are used to extrapolate for N=1600.

Hybrid architecture has structural dependency, i.e, it cannot execute DWT and DCT at the same time, even if there is no data dependency. Architecture with standalone DWT and standalone DCT doesn’t have structural dependency. i.e, it can run DWT and DCT at the same time, if there is no data dependency between currently running DCT and DWT. In the watermarking algorithm, DCT is computed after two-level DWT computation, i.e, there is data dependency between DCT and DWT.

The DFT implementation requires 2N + 4N₂ multipliers, 2N + 4N₂ + 2 adders. For DCT and DST implementation, it requires 2N + 2N₂ + 1 adders, 2N + 2N₂ multipliers. For DWT implementation, it requires M² adders, M² multipliers, where M= $\sqrt{N}$ . For N=16, DCT computation requires 48 multipliers, 50 adders and DCT computation requires 16 multipliers, 16 adders. Let one clock cycle=T_A. The time taken to run one set (8 values) of DWT is 9*T_A time units and on set (16 values) of DCT is 27*T_A time units. In the proposed architecture, T_A=10.

Once the architecture computes first set of DWT values, the standalone architecture can start computing first set of DCT and second set of DWT at the same time, whereas the hybrid architecture cannot. It can either compute 1st set of DCT or second set of DWT values. In this case, the second set of DWT is computed first and then DCT is computed as shown in Table 3. It can be observed that, time saved for two frames by standalone architecture is 270 time units which is the time required to run one DCT. Time saved for n frames is (n-1)*270. For N=1600, no. of DWT runs to be computed is 200 and no. of DCT runs to be computed is 25.

Table 3
Comparison between hybrid and standalone

Hybrid architecture Standalone architecture

1st set of DWT (0-720) 1st set of DWT (0-720)

2nd set of DWT (721-1440) 2nd set of DWT (721-1440)

1st set of DCT (1441-1710) 1st set of DCT (721-990)

T2nd set of DCT (1711-1980) 2nd set of DCT (1441-1710)

Hybrid architecture	Standalone architecture
1st set of DWT (0-720)	1st set of DWT (0-720)
2nd set of DWT (721-1440)	2nd set of DWT (721-1440)
1st set of DCT (1441-1710)	1st set of DCT (721-990)
T2nd set of DCT (1711-1980)	2nd set of DCT (1441-1710)

The data is extrapolated for 1025 frames with each frame having 1600 samples. As per Table 4, Hybrid architecture needs 35.47% more time. As DCT and DWT computations are done serially on Hybrid Architecture, more time is required when compared to standalone architecture which can compute DCT and DWT parallelly.

Table 4

Time comparison for N=1600

Simulation	Hybrid	Standalone
DWT for 1-Frame	18,000(90x200)	18,000(90x200)
DCT for 1-Frame	6,750(270x25)	270(270x1)
Total time for 1-Frame	24,750	18,270
Total time for 1025-Frames	2,53,68,750	1,87,26,750

4.3 Error between Matlab and verilog results

The algorithm is simulated for two frames containing 64 samples each as shown in Fig. 4 and the level-2 DWT and mean values of DCT are captured and compared. The average error for outputs of 2 level DWT outputs shown in Table 5 is 0.012785% and for DCT outputs shown in Table 6 is 0.0395%.

Fig. 4

Algorithm.

Table 5

DWT outputs (round off to 6 digits)

Matlab output	Verilog output	Error %
1.412525	1.412415	–0.008
1.252898	1.252808	–0.007
0.743510	0.743408	–0.014
0.849776	0.849670	–0.012
1.094192	1.094116	–0.007
1.183180	1.183105	–0.006
0.836847	0.836731	–0.014
0.820559	0.820435	–0.015
1.225195	1.225037	–0.013
1.069751	1.069641	–0.010
0.547769	0.547668	–0.018
0.724341	0.724182	–0.022
1.025594	1.025513	–0.008
0.575747	0.575562	–0.032
0.824274	0.824158	–0.014
1.357872	1.357788	–0.006
1.192364	1.192261	–0.009
1.163867	1.163757	–0.009
1.079697	1.079590	–0.010
1.480276	1.480164	–0.008
0.992737	0.992615	–0.012
1.424531	1.424438	–0.007
1.146984	1.146851	–0.012
0.628918	0.628784	–0.021
1.263610	1.263550	–0.005
1.302949	1.302795	–0.012
0.586852	0.586731	–0.021
0.867590	0.867432	–0.018
0.836658	0.836548	–0.013
0.798535	0.798401	–0.017
1.002538	1.002319	–0.022
1.351594	1.351501	–0.007

Table 6

DCT outputs (round off to 6 digits)

Matlab output	Verilog output	Error %
1.174599	1.174377	-0.019
1.081090	1.080444	-0.060

5 Conclusions

A hybrid architecture for multiple transforms has been proposed. Implemented Hybrid architecture is verified against the values obtained in Matlab. A watermarking application which uses DWT and DCT is identified. The algorithm sequence is implemented in verilog.

The proposed hybrid architecture occupies less area and has lesser resource utilization. The designed architecture produces the output of transform sequence in order, which avoids reordering the output.The implemented architecture can be used to handle large sized transforms by repeatedly using fixed size architecture for a large number of points without increasing the number of Processing Elements (PEs). Hybrid architecture needs 35.47% more time. There is 22.14%, 16.83%, 9.09% and 23.21% reduction in LUT’s, FF’s, DSP units and IOB’s respectively in the hybrid architecture.

As future work, we want to modify the architecture so that an application can make use of maximum resource in parallel to have an higher yield.

Footnotes

Acknowledgment

We would like to express our deepest gratitude to our collaborator Dr. Mamatha.I from Amrita University, whose constructive suggestions helped to speed up the implementation.

References

Kung

H.-T.

, Why systolic architectures?, IEEE computer15(1) (1982), 37–46.

Aravena

J.L.

, Triple matrix product architectures for fast signal processing, IEEE transactions on circuits and systems35(1) (1988), 119–122.

Chiper

D.-F.

, Swamy

M.N.S.

, Ahmad

M.O.

and Stouraitis

, Systolic algorithms and a memory-based design approach for a unified architecture for the computation of DCT/DST/IDCT/IDST, IEEE Transactions on Circuits and Systems I: Regular Papers52(6) (2005), 1125–1137.

Wahid

K.A.

, Islam

Md.A.

, Shimu

S.S.

, Lee

M.H.

and Ko

S.-B.

, Hybrid Architecture and VLSI Implementation of the Cosine–Fourier–Haar Transforms, Circuits, Systems and Signal Processing29(6) (2010), 1193–1205.

Mamatha

, Sudarshan

T.S.B.

, Tripathi

and Bhattar

, Triple-Matrix Product-Based 2D Systolic Implementation of Discrete Fourier Transform, Circuits, Systems, and Signal Processing34(10) (2015), 3221–3239.

, Min

, Xiaoming

, Zhiguo

and Cheng

, Novel zero-watermarking scheme based on DWT-DCT, China Communications13(7) (2016), 122–126.