WO2001018743A9 - Fast and efficient computation of cubic-spline interpolation for data compression - Google Patents

Abstract

A fast and efficient method and system for computation of cubic-spline interpolation for data compression is described (II). In one aspect, the present invention is a method and system for obtaining a transform of the compressed signal and computing an inverse transform of the transform of the compressed signal to obtain the compressed signal (I).

Description

FAST AND EFFICIENT COMPUTATION OF CUBIC-SPLINE INTERPOLATION

FOR DATA COMPRESSION

FIELD OF THE INVENTION

This invention relates to data compression. More specifically, the invention relates to a new cubic-spline interpolation (CSI) for both 1-D and 2-D signals to sub-sample signal and image compression data.

BACKGROUND OF THE INVENTION

In most multimedia systems the amount of image data is so large that the use of image data compression is almost mandatory. Image data compression allows the image to be transmitted over the Internet in real time. Also it reduces the requirements for image storage. Presently both spatial and temporal data reduction techniques are available and continue to improve the performance of image data compression. The fundamental problem of image data compression is to increase the compression ratio and to reduce the computational complexity within an acceptable fidelity.

Interpolation is one of the more important functions that can be used in the process of estimating the intermediate values of a set of discrete sampling points. Interpolation is used extensively in image data compression to magnify or reduce images and to correct spatial distortions. For example, see R. G. Keys, "Cubic Convolution Interpolation for Digital Image Processing," IEEE Trans, on Acoustics, Speech, and Signal Processing, vol. ASSP-29, no.6, pp.l 153-1160, Dec. 1981, [1], the contents of which are hereby expressly incorporated by reference. In general, the process of decreasing the data rate is called decimation and the process of increasing data samples is called interpolation as described in H. S. Hou, and H. C. Andrews, "Cubic Splines for Image Interpolation and Digital Filtering," IEEE Trans, on Acoustics, Speech, and Signal Processing, vol. ASSP-26, no.6, pp.508-517, Dec. 1978, [2], the contents of which are hereby expressly incorporated by reference. It is well known that several interpolation functions such as linear interpolation (see, . K. Pratt, Digital Image Processing, second edition, John Wiley & Sons, Inc., New York, 1991, [3], the contents of which are hereby expressly incorporated by reference.) cubic-convolution interpolation (see [1], and [3]), cubic B-spline interpolation (described in C. de Boor, A Practical Guide to Splines. New York: Springer-Verlag, 1978, [4]; M. Unser, A. Aldroubi, and M. Eden, "B-Spline Signal Processing: Part II-Efficient Design and Applications," IEEE Trans, on Signal Processing, vol.41, pp.834-848, Feb. 1993, [5]; M. Unser, A. Aldroubi, and M. Eden, "Enlargement or Reduction of Digital Images with Minimum Loss of Information," IEEE Trans, on Image Processing, vol.4, pp.247-258, March 1995, [6]; and [2]) can be used in the image data compression process.

The disadvantage of these interpolation schemes is that in general they are not designed to minimize the error between the original image and its reconstructed image. In 1981 Reed (I. S. Reed, Notes on Image Data Compression Using Linear Spline

Interpolation, Department of Electrical Engineering, University of Southern California, Los Angeles, California, 90089-2565, U.S.A., Nov. 1981 [7], the contents of which are hereby incorporated by reference) and in 1998 Reed and Yu (I. S. Reed and A. Yu, Optimal Spline Interpolation for Image Compression, United States Patent, No. 5822456, Oct. 13, 1998 [8], the contents of which are hereby incorporated by reference) developed a linear spline interpolation scheme for re-sampling the image data. This linear spline interpolation is based on the least-squares method with the linear interpolation function.

Using an extension of the ideas of Reed in [7,8], a modified linear spline interpolation algorithm, called the cubic-spline interpolation (CSI) algorithm, is developed in this invention for the sub-sampling of image data. (The linear spline interpolation explained in [8] and used by America On Line™ (AOL) will be called the "AOL algorithm" in this document from hereon.)

It follows from [1] that the cubic-convolution interpolation, which is different from the B-spline interpolation, can be performed much more efficiently than that of the cubic B-spline interpolation method. In this invention, the new CSI scheme combines the least- squares method with a cubic-spline function developed by Keys [1] for the decimation process. Also the cubic-spline reconstruction is used in the interpolation process. Therefore, the CSI constitutes a new scheme that is quite different from both cubic B- spline interpolation [2,4-6] and cubic-convolution interpolation [1,3]. The concept of the CSI for both 1-D and 2-D signals is describes and demonstrated in the following sections. In addition, it is shown by computer simulation that the CSI scheme obtains a better subjective quality for the reconstructed image than linear interpolation, cubic-convolution interpolation, cubic B-spline interpolation and linear spline interpolation. An important advantage of this new CSI scheme is that it can be computed by a use of the FFT technique. The complexity of the calculation of the CSI scheme is substantially less than other conventional means.

W. B. Pennebaker and J. L. Mitchell, JPEG Still Image Data Compression Standard, Van Nostrand Reinhold, New York, 1993, [9], the contents of which are hereby incorporated by reference, describes the JPEG still image data compression standard. It is well known that the JPEG (see [9]) algorithm is the international compression standard for still-images. The disadvantage of the conventional JPEG algorithm is that it causes visually disturbing blocking effects when high quantization parameter is used to obtain a high compression ratio. One embodiment of this invention includes a simpler and modified JPEG encoder-decoder to improve the JPEG standard with a high compression ratio and still maintain a good quality reconstructed image.

Recently, the authors in T. K. Truong, L. J. Wang, I. S. Reed, W. S. Hsieh, and T.C. Cheng "Image data compression using cubic convolution spline interpolation," accepted or publication in IEEE Transactions on Image Processing [10], the contents of which are hereby incorporated by reference, proposed the modified JPEG encoder-decoder for τ= 2 that utilizes the CSI scheme with a compression ratio of 4 to 1 as a pre-processing stage of the JPEG encoder and the cubic-spline reconstruction with a ratio of 1 to 4 as a post- processing stage of the inverse JPEG decoder to achieve a high compression ratio.

In such a modified JPEG encoder the CSI scheme is the pre-processing stage of the JPEG encoder. It can be implemented by the use of the FFT algorithm. In addition, the output of the modified JPEG encoder represents the compressed data to be transmitted. It can be pre-computed and stored. In such a modified JPEG decoder, the cubic-spline reconstruction constitutes the post-processing stage of the JPEG decoder. This postprocessing stage is different from the conventional post-processing algorithms that were proposed to reduce the blocking effects of block-based coding in B. Ramamurthi and A. Gersho, "Nonlinear space variant post-processing of block coded images," IEEE Trans, on Acoustics, Speech, Signal Processing, vol. ASSP-34, pp.1258-1267, 1986, [11], Y. Yang, N. Galatsanos, and A. Katsaggelos, "Projection-based spatially adaptive reconstruction of block-transform compressed images," IEEE Trans, on Image Processing, vol.4, pp.896-908, July 1995 [12], the contents of which are hereby incorporated by reference.

The proposed post-processing stage is an interpolation process that uses the cubic- convolution interpolation. In [10], a computer simulation shows that the modified JPEG encoder-decoder for τ= 2 obtains a better subjective quality and an objective PSNR of the reconstructed image than the JPEG algorithm described in T. Lane, Independent JPEG Group's free JPEG software, 1998, [13], the contents of which are hereby incorporated by reference; and [9]. Furthermore, the modified inverse JPEG decoder requires less computational time than the conventional JPEG decoder. But, the disadvantage of the modified JPEG encoder-decoder for τ— 2 is that the computational time required for the modified JPEG encoder is greater than the conventional JPEG encoder.

Thus, in one aspect, the present invention describes a fast method to compute the modified JPEG encoder. It is shown in this aspect of the invention that the speed of the new method for computing the modified JPEG encoder is approximately two times faster than that of the conventional JPEG encoder with still a good quality of reconstructed image. SUMMARY OF THE INVENTION

The present invention describes a new cubic-spline interpolation (CSI) for both 1-D and 2-D signals to sub-sample signal and image compression data. This new inteφolation scheme which is based on the least-squares method with a cubic-spline function can be implemented by the fast Fourier transform (FFT), and/or by a Winograd discrete Fourier transform (WDFT). The result is a simpler and faster inteφolation design than can be obtained by conventional means. It is shown by computer simulation that such a new CSI yields the most accurate algorithms for smoothing. Linear inteφolation, linear spline inteφolation, cubic-convolution inteφolation and cubic B-spline inteφolation tend to be inferior in performance. In addition it is shown in this invention that the CSI scheme can be performed by a fast and efficient computation. The proposed method uses a simpler technique in the decimation process. It requires substantially fewer additions and multiplications than the original CSI algorithm. In one aspect, the present invention is a method and system for defining a cubic- spline filter; correlating the filter with the signal to obtain a correlated signal; autocorrelating the filter to obtain autocorrelated filter coefficients; computing a transform of the correlated signal and the autocorrelated filter coefficients; dividing the transform of the correlated signal by the transform of the autocorrelated filter coefficients to obtain a transform of a compressed signal; and computing an inverse transform of the transform of the compressed signal to obtain the compressed signal. The signal, the filter, and the transforms may be one dimensional or two dimensional. Further, the transforms may be a fast Fourier transform (FFT) or a Winograd discrete Fourier transform (WDFT) with an overlap-save scheme. Also, a zonal filter may be defined to simplify the steps of correlating and autocorrelating.

Furthermore, a new type of overlap-save scheme can be utilized to solve the boundary-condition problems that occur between two neighboring sub-images in the actual image for higher compression ratios. It is also shown in this invention that a very efficient 9-point Winograd discrete Fourier transform (WDFT) can be used to replace the FFT needed to implement the CSI scheme image for higher compression ratio of 9 to 1. Finally, a fast new CSI algorithm is used along with the Joint Photographic Experts Group (JPEG) standard to design a modified JPEG encoder-decoder for image data compression. As a consequence, for the higher compression ratios the proposed modified JPEG encoder-decoder obtains a better quality of reconstructed image and also requires less computational time than both the conventional JPEG method and the America on Line (AOL) algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

The objects, advantages and features of this invention will become more apparent from a consideration of the following detailed description and the drawings, in which:

FIG. 1 is an exemplary periodic function with a period of τ= 5;

FIG. 2 is an exemplary 1-D cubic-spline function;

FIG. 3 is an exemplary shifted cubic-spline function for n=6;

FIG. 4 is an exemplary side view of a 2-D cubic-spline function; FIG. 5 is an exemplary reconstructed function between sampling periods;

FIGs. 6 (a) - (d) are exemplary zonal masks of the 2-D cubic-spline function used in the computation of Y_{Jt h} in (25) and A_JιJ2tkιk2 in (26) for the CSI scheme;

FIG. 7 is an exemplary overlap map grid points for calculating α for τ=3; FIG. 8 is an exemplary overlap map grid points for calculating β for τ=3; FIG. 9 is an exemplary overlap map grid points for calculating γ for τ=3;

FIG. 10 is an exemplary overlap map grid points for calculating δ for τ=3; FIG. I l is a l9^χ19 sub-images of size 9><9 in an exemplary image of size 171 x 171; FIG. 12 is a reconstructed image with serious artifacts that is generated by using FCSI method implemented by the direct use of the 9 x 9 Winograd DFT for compression; FIG. 13 is an illustrative example of the FCSI algorithm implemented by a 5x5

Winograd DFT and the overlap-save sub-image method;

FIG. 14 is a reconstructed image with no apparent artifacts that is generated by using a FCSI implemented by the 9x9 Winograd DFT and the overlap-save method for compression; FIG. 15 is an exemplary modified JPEG encoder for τ= 2, 3

FIG. 16 is an exemplary modified JPEG decoder for τ= 2, 3;

FIGs. 17 (a) - (d) illustrate some reconstructed images with a compression ratio of 100:1; and

FIGs. 18 (a) - (d) illustrate some reconstructed images with a compression ratio of 200:1.

DETAILED DESCRIPTION

The present invention describes a new cubic-spline inteφolation (CSI) for both 1-D and 2-D signals to sub-sample signal and image compression data. This new inteφolation scheme which is based on the least-squares method with a cubic-spline function can be implemented by the fast Fourier transform (FFT). Furthermore, in one embodiment of the present invention, a new type of overlap-save scheme is utilized to solve the boundary- condition problems that occur between two neighboring sub-images in the actual image for higher compression ratios. In one embodiment of this invention, an efficient 9-point

Winograd discrete Fourier transform (WDFT) is used to replace the FFT needed to implement the CSI scheme image for higher compression ratio of 9 to 1. Finally, a fast new CSI algorithm is used along with the Joint Photographic Experts Group (JPEG) standard to design a modified JPEG encoder-decoder for image data compression.

In order to speed up the modified JPEG encoder for τ= 2 with a good image quality, the compression ratio of the CSI scheme is extended to 9 to 1 (τ= 3). However, if the compression ratio for the CSI scheme is increased, there is need for additional computations that involve considerably more additions and multiplications. Therefore, in this invention a new faster and efficient algorithm for CSI is developed. The basic idea of this new algorithm, called the fast cubic-spline inteφolation (FCSI) scheme, is based on the CSI scheme, but it has a simpler form than that used for the original CSI scheme. The FCSI scheme substantially reduces the complexity of the additional computations that are required for the increased compression ratio. Moreover, the constants a,β,γ and δ , used to compute the terms AJJ _{k k} for the new FCSI scheme, are accurately calculated in detail in this invention.

For τ= 3, the actual image to be compressed, a novel overlap-save sub-image method is used to solve for the needed boundary conditions. Also a simplified and efficient algorithm that uses a 9-point Winograd DFT (WDFT) instead of the FFT is implemented to compress an actual image. Computer runs show that for some gray images of size 512 by 512, the computational time of the FCSI encoder that is implemented by a new type of overlap-save sub-image method and the 9-point Winograd DFT on a 400-MHz Intel Pentium II personal computer using C-code can be dramatically reduced. The FCSI encoder requires only around 0.15 sec when compared with around 0.57 sec for the original CSI encoder at a compression ratio of 9 to 1 (τ= 3). Also, the new FCSI scheme obtains a PSNR that is similar to all of the other more complicated zonal filters considered in this invention. Finally, this FCSI scheme with τ= 3 is combined with the JPEG standard to speed up the modified JPEG encoder for color image coding and still obtain a better quality of reconstructed image than the JPEG algorithm for higher compression ratios. In other words, the modified JPEG encoder for τ= 3 requires 0.71 sec, 0.38 sec and 0.67 sec less time than the modified JPEG encoder for τ= 2 [10], the conventional JPEG encoder [13] and the AOL algorithm [8], respectively.

This document is organized as follows: In Section II the encoding algorithm for this new inteφolation method is derived in detail. In addition it is shown now that the properties of the FFT and the convolution theorem can be used to compute the CSI scheme. The decoding algorithm is explained in Section III. In Section IN the FCSI computation is described. The constants needed for the FCSI are calculated in Section N. A novel FCSI algorithm that uses a new type of overlap-save sub-image technique and the Winograd DFT algorithm is developed in Section NL In Section Nil the modified JPEG encoder-decoder is presented. Finally, experimental results are presented in Section NIII. II. ENCODING ALGORITHM FOR THE NEW INTERPOLATION METHOD

Encoding with the CSI scheme, utilizes the decimation process needed to perform image data compression. The philosophy of the CSI scheme is to recalculate the sampled values of the signal or image data by means of the least-squares method using the cubic- spline function. It is shown in this section that this new proposed method applies to both 1-D and 2-D signals as follows: A. CSI for the 1-D Signal

Let τ be a fixed, positive integer. Also let the data function X(t) be periodic with period nτ , where n is an integer. An example of X(t) for n = 5 is illustrated in FIG. 1 From [1], the 1-D cubic- spline function R(t) , shown in FIG. 2, is defined by

(3/2)|t|³ - (5/2)|t|² + 1 , 0 ≤ |t| < l

R(t) = ] -(l/2)|t|³ +(5/2)|t|² -4|t| + 2 , l ≤ |t| < 2 (1)

0 , 2 < |t|

Next one needs to define the shift function of the cubic-spline function R(t) as follows: Ψ_t( = R(t -kτ) for 0 ≤ Jfc ≤ n-1 . (2)

An example of the shift functions Ψ_k(t) for n = 6 is shown in FIG. 3. The goal is to approximate X(t) by a π-point sum, given by

S(t) = ∑X_kΨ_k (t) = ∑X_kR(t -kτ) (3) i=0 k=0 in a least-squares fashion, where X_ϋ,---,X_n__x are the reconstructed values at the sampling points which represent the compressed data to be transmitted or stored. The function S(t) in (3) is cubic-spline reconstruction of the function X(t) using the weights ₀,---, „_, . It follows from [7] that the least-squares approximation of S(t) to X(t) is defined by

L(X₀,X_l,-X_lt__l) = ∑(X(t) -S(t))² = ∑ (X(t) ~ ∑X_kΨ_k(t))², (4) where the sum is taken over one period nτ plus 2τ of the data. FIG. 1 shows an exemplary periodic function with a period 5τ , FIG. 2 is an exemplary 1-D cubic-spline function, and FIG. 3 is an exemplary shifted cubic-spline function for n = 6.

Using the same procedure, described in [7,8], one can find the weights X_a,X_λ,---,X_n__x in (4) that minimize the function L(X₀,X_l ,---X_n__l) . To minimize (4) partial differentiation of L(X₀ , X_λ , ■ ■ ■ X_n_ ) with respect to _j for 0 ≤j ≤n -l yields the following set of equations:

or

where

and

Y_j = ∑X(t)Ψ_j(t) for 0≤j≤n-l. (7) t=-2τ

The term 7 in (7) can be reduced as follows: nτ 0^'+2)r

Y_j = ∑ X(t)R(t -jτ)= ∑ X(t)R(t -jτ). (8) t=-2τ t=(j-2)τ

Let t-jτ- m , then

Note that computing Y. in (9) involves n correlation coefficients of only 3τ-l points.

Now let the periodicity of the cubic-spline function be R(t) = R(t + nτ) , i.e., R(t) has a period of nτ. Then the matrix form of A_jk in (6) can be simplified to:

2r-l v._* = ∑ R(m + (j-k)τ)R(m) ra=-2r+l

(10)

It follows from [7] that A_J>k in (10) can be expressed in circulant form by letting

where (k-j)_n denotes the residue (k - j) modulo n, and

B_o=a,B =β,B₂=γ,B₃=δ,B₄=0, ^■■■, B„_₄ = 0 , B_n_₃ = δ , B„_₂ = γ , (12)

B„-ι=β-

Therefore, A_jk in (11) and (12) has the following symmetric, circulant representation:

A substitution of (13) into (5) yields the matrix equation,

A-X = Y, (14) where the matrix A is given in (13), X = (X₀ ,X_X,---, X_n__x )^τ and Y = (Y₀,Y_X,---, Y_n__x )^τ .

Since the nxn matrix on the left side of (14) is a circulant matrix, (14) reduces immediately to n-l

J = ∑ H). forO≤j≤ra-1. (15) k=0

But 5_;=5„_. for t = l,2,---,[«/2j, where |_ j denotes the greatest integer less than or equal to x . Thus (15) becomes r,=∑ * -». for0≤j<n-l. (16)

In (16) the FFT can be used to solve for the X_k . To see this, let the FFT of Y_j , X_k and B_j for 0≤j,k,m≤n-l be defined by Y_m , X_m and B_m , respectively. By using the convolution theorem describes in E. O. Brigham, The Fast Fourier Transform and its Application, Prentice-Hall International, Inc., Englewood Cliffs, New Jersey, 1988 [14]; and A. V. Oppenheim, and R. W. Schafer, Digital Signal Processing, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1975 [15], the contents of which are hereby expressly incoφorated be reference, it is easy to see that the solution of (16) in the frequency domain can be expressed as Ϋ_m = X_m ■ B_m or X_m = Y B_m , where B_m ≠ 0. Thus, using the inverse FFT of X_m , the X_k can be obtained for 0≤k≤n-l.

The encoding method for the 1-D signal is summarized as follows:

Choose a suitable value of integer r . The compression ratio is roughly τ . Apply (9) to find Y_} . Also apply (10), (11) and (12) to find B_j .

Find the FFT of Y_} and B_j to obtain Y_m and B_m, respectively. Also compute

Take the inverse FFT of X_m to get X_k which is the compressed data to be transmitted or stored. B. CSI for the 2-D Signal

Let X(t_x,t₂) be a doubly periodic signal (e.g., image) of periods nτ and n₂τ with respect to integer variables t_x and t₂ , where n and n₂ are also integers. A 2-D cubic- spline function, R(t_l,t₂), is defined by

R(t_x,t₂) = R(t)-R(t₂), (17) where R(t_x) and R(t₂) are 1-D cubic-spline functions, respectively. A 3-D plot of this cubic-spline function is shown in FIG.4. It is well known fact that 2-D inteφolation can be accomplished by the use of 1-D inteφolations with respect to each coordinate [1,3]. By analogy with the 1-D case, let _{k kl}(t_x,t₂)= R(t_x-k_xτ,t₂-k₂τ)

= R(t_x - k_x τ) ■ R(t₂ - k₂τ) for 0 < k_t < n_t - 1 and i = 1, 2. (18)

By a procedure similar to that used in the 1-D case in (3), one defines the 2-D CSI as follows:

where X_{k /c} are the reconstructed values at sampling points which represent the compressed image to be transmitted or stored. Again one wants to find the best weights X_{k k} such that

is a minimum. It follows from [7] that minimizing (20) yields

' "² = Σ Σ 2(X(t_x,t₂)-∑ ∑X_{kl,k2 k2}(t^₂))Ψ_ΛJ2(t,t₂) = 0, ^aAh,h 'ι=-^2r '2=-2«- *ι=0 *.=<> or

Σ Σ

,- = 1,2. (21) where n,r «₂r ΛΛ,*A = Σ ∑Ψ_Λ,Λ('I )Ψ .('I ), O≤Λ ≤^-1 and ^I = 1,2, (22) ,=-2 (=-2r and

X,,_Λ = O≤Ji ≤*ι - and i = l,2. (23)

The term Y_j __; in (23) reduces as follows:

^ Λ_,Λ_Λ= Σ £ X(t_χtt₂)R(f -j_xτ,t₂~j₂τ)

(,=-2r /,=-2τ

(24)

0^'ι+2)r ( ₂+2)r

= ∑ Σ X(t_x,t₂)R(t_x-j_xτ,t₂-j₂τ) h=Uι-2)τ t₂=J₂-2)τ

Let t_k ~ j_kτ = m_k for k = 1,2 , then 2τ-l 2τ-l

Λ l = Σ Σ X(m_x+ j_xτ,m₂+ j₂τ)R(m_x,m₂), (25) ,=-2r+l m₂=-2r+l where R(m_x,m₂) is the 2-D cubic-spline function shown in FIG.4. The terms AJJ _kk in (22) are treated in a similar fashion as follows:

2r-l 2r-l jJl ≠l = ∑ ∑ R(m ₊(j_x ^■ k )τ, m₂ + U2-k₂)^τ)^{R m m}2) (26) m,=-2r+l m₂=-2r+l

From Eq. (26) given above, the array A _kk can be expressed in 2-D circulant form as follows:

A

^A JhzAk ^~B(.h-λ⁾ _nι, fe ^■h)_r, (27) where (k_t - j_i)_n denotes the residue (k_: - j ) modulo «,• for z^' = l,2 and

where 0 < s₍. < n.- - 1 for z^' = 1,2. Note that if the array [ B ] is represented in matrix form, it is a block-circulant matrix. FIG.4 is the side view of the 2-D cubic-spline function. Since the matrix [B ] in (28) is a block-circulant matrix, (21) can be expressed by

«< -1 «₂ -1 h>h ^: Σ Σ "^ a B, C*I-Λ⁾«,, - 1 and i = 1,2 , (29) =04,=0 ⁽*--Λ⁾«, 0 ≤ j ≤ n, where (k_: -y^',-)_π/ denotes (k_t-j^mod ; for i = l,2. Using a procedure similar to that used in (16), then (29) becomes

H_|-lκ₂-l

^Yj,,J2 ⁼ ΣΣ t,,*! ^BUl -ftl)„, , (Λ-*2)»₂ 0 < y^',- < «,. - 1 and z = 1,2. (30) i-,=0/c₂=0

In (30) the 2-D FFT [14,15] can.be used to solve for the X_kh . Let the 2-D FFT of

7 Λ,J2 ^; X_{k k} zrA B for 0 < j, ,k_t ≤ «,• - 1 and O≤ i≤n,-!, i = l,2 be defined, respectively, by Y_mn , X_m>n and B_mn for 0 < m ≤ n_x -1 , 0 < n < n₂ -1. Then Eq. (30) in the frequency domain can be expressed as Ϋ_m>n= _mιn-B_mι„ or X,,_un-Y_mtn /£,„_,„, where B_mn ≠ 0. Finally, using the 2-D function of the inverse FFT of X_mn , the X_{k >k} are obtained for 0 < k, < n_i - 1 , i - 1,2.

The encoding method for the 2-D signal is summary in the following steps:

• Choose a suitable value of integer τ . The compression ratio is roughly τ¹.

• Apply (25) to find Y_ΛJι . Next, apply (26), (27) and (28) to find _3,_ιΛ .

• Take the 2-D FFT of 7 . and B to obtain Y_mn and B_mn, respectively. Also compute X„=^/^_B.

• Take the inverse FFT of X_mn to obtain X_{k >k} which is the compressed image to be transmitted or stored.

III. DECODING ALGORITHM

In the decoding process, using the reconstructed values at the sampling points (e.g., X_k and X_{k <k} ) obtained in Section II, the reconstructed points between the sampling points are obtained by means of the cubic-spline function. This decoding algorithm is called the cubic-spline reconstruction.

A. Decoding of the compressed 1-D Signal Since the n reconstruction values X₀,---,X_n__x are known, the reconstructed signal

S(t) can be obtained by the use of (3). In other words, the retrieved signal is the convolution of the cubic-spline function R(t) , defined in (1), and the sequence of n reconstructed values with sampling interval τ . The reconstructed function S(t_a ) between the two adjacent reconstructed values X_k and X_k+X is illustrated in FIG. 5 and given by the sum,

S(t = X_k__xΨ_k__x (t_a) + X_kΨ_k(t_a) + X_MΨ_k+x (t_a) + X_k+2Ψ_k÷2(t_a) , (31) where kτ < t_a < (k + l)r and Ψ^ (t) is defined in (2) and the boundary conditions as given in [1] are X_ = 3(X₀ -X_x) + X₂ and X_n = 3(X_n__x - X_n_₂ ) +X„_₃. FIG. 5 is a reconstructed function between the sampling points.

B. Decoding of the compressed 2-D Signal

Since the reconstruction values X_kl>kl for 0 ≤ / ,. ≤ n, -l, t = l,2 are known, the 2-D reconstructed image S(t_x,t₂) can be obtained by the use of (19). In other words, the retrieved image is the 2-D convolution of the 2-D cubic-spline function R(t_x,t₂) , given in

(17), and the 2-D sampled waveform X_{k k} . It is well known that the computationally simpler method described in [3], called bilinear inteφolation, can be used also to perform the 2-D inteφolation. The ideas of the bilinear inteφolation can be used to accomplish the 2-D cubic-spline reconstruction. In other words, the discrete data of each row can be inteφolated from the reconstructed values X_{k k} with similar inteφolation for the given discrete data of each column. It is easy to show that the reconstructed image S(t,,t ) between the four adjacent reconstructed values is given by

where k_xτ < t_t < (k + l)τ , k₂τ < t_j < (k₂ +l)τ and the boundary conditions as given in [1] are -U₂

~ k„, ) ⁺ k_t,2 > ^Xk_s,n

'

X_n -_\ - X_n-_\ -_\ ^_ ^«-2,-ι ) ⁺ ^«-₃,-i ' ^-ι,« ⁼ 3( ₀ „ - X_{X n} ) + X_{2 n} , X_n>n = ( „_ι „ -

IV. FAST COMPUTATION OF CSI

In fact, the CSI scheme needs a large number of pixels of the 2-D cubic-spline function R(t_x,t₂) in order to compute the Y_{j h} in (25) and the A_jjijkιk2 in (26). For example, computing each 7, • in (25) and JJ _tk]kι in (26) for τ= 2 (compression ratio 4: 1) and τ= 3 (compression ratio 9: 1) involve 81 and 169 pixels of the domain of function R(t_x ,t₂) , respectively. As a consequence, the complexity of the calculation of (25) and (26) needed in the CSI scheme is substantially increased when the compression ratio is extended from τ= 2 to τ= 3. To illustrate this, for τ= 3 some abbreviated solutions are used to overcome the above problem of extra computational complexity. First, one defines a zonal mask of 169 pixels of R(t ,t₂) shown in FIG. 6(a), called zonal filter 1. This means that to compute each Y_{j h} in (25) one needs to use a zonal mask of 169 pixels of R(t ,t₂) to correlate with a periodic image, X(t_x,t₂) . Also, to calculate each A_{] j k}^ in (26) one needs to use the 169 pixels for the function R(t_x,t₂) to autocorrelate with the other 169 pixels of R(t_x,t₂) .

FIGs. 6(a) -(d) are zonal masks of the 2-D cubic-spline function used in the computation of Y_{J Jι} in (25) and JJ _{k k} in (26) for the CSI scheme. FIG. 6(a) shows a 169 pixels in zonal filter 1, FIG. 6(b) shows a 133 pixels in zonal filter 2, FIG. 6(c) shows a 69 pixels in zonal filter 3, and FIG. 6(d) shows a 25 pixels in zonal filter 4. In order to reduce the computational complexity of this zonal filter 1, it is proposed next to use the zonal filter 2, shown in FIG. 6(b), the zonal filter 3, shown in FIG. 6(c), and the zonal filter 4, shown in FIG. 6(d). These zonal filters 2, 3 and 4 use zonal masks of 133, 69 and 25 pixels or grid points of R(t_x,t₂) , respectively, to compute each -₂ in (25) and each

^Af k₂ ⁱⁿ (²⁶)-

Next for τ= 3, the FCSI schemes described above are developed by the use of zonal filter 2, zonal filter 3 and zonal filter 4. Following the same procedure, used in the previous section, the FCSI scheme for τ= 3 is obtained readily. First the derivations of these FCSI schemes for τ ~ 3 are developed by a means similar to that of the original CSI scheme for τ= 2 [10]. Only the terms Y_{j J}, in (25) and A_JJιAki in (26) for the FCSI schemes are slightly different from the original CSI scheme. Also in the FCSI scheme, the complexities of the terms Y_{J Jι} and A_{jj k k} for zonal filter 2, zonal filter 3 and zonal filter 4 are summarized in Tables I and II and are described in the following: 1. By FIG. 6(a), the terms 7^ and A_{hh> kι} for zonal filter 1 are obtained completely in a manner similar to (25) and (26). It uses 169 pixels of R(t_x,t₂) to compute each Y_{jι h} and each A_{hh k kϊ} . This algorithm is very complex and is just the original CSI scheme.

2. By FIG. 6(b), the terms JJ _{k k} for zonal filter 2 are obtained, and the result is similar to (26), except that σ = 0. Each 7, • in (25) is computed for zonal filter 2 with 133 pixels of R(t_x,t₂) . This algorithm still involves extensive computations.

3. By FIG. 6(c), the terms A_{j htk]kι} for zonal filter 3 are obtained, and the result is similar to (26), but for this case μ = 0 and σ = 0 . Each 7₇ ■ in (25) is computed for zonal filter 3 with 69 pixels of R(t_x,t₂) . This algorithm is also complex.

4. By FIG. 6(d), the terms A_JJιAkι for zonal filter 4 are obtained, and the result is similar to (26), but in this case η = = p = μ = σ = 0 . Each Y_{j h} in (25) is computed for zonal filter 4 by the use of only 25 pixels of R(t_x,t₂) . This algorithm is more compact than any of the other computations given above in cases 1, 2 and 3. TABLE I COMPLEXITY OF A_hh USING FOUR ZONAL FILTERS.

TABLE II

COMPLEXITY OF COMPUTING

WITH A COMPRESSION RATIO OF 9 TO 1 FROM AN

IMAGE OF SIZE 512 BY 512 USING FOUR ZONAL FILTERS.

It is shown by computer runs in Section VIII that the zonal filter 4 obtains a PSNR that is similar to any of the other three zonal filters. Thus, this zonal filter 4 represents the most practical and simple zonal filter for the FCSI scheme. The primary advantage of the

FCSI scheme with zonal filter 4 over the original CSI scheme is that it substantially reduces the computational complexity.

V. CALCULATION OF CONSTANTS

In accord with (27), (28) and the previous section the terms _{j jltkιkι} needed for zonal filter 4 are given by the general formula: if y^' _j - k_x ≡ 0 mod ;z_j and y^' ₂ - k₂ ≡ 0 mod n₂ if y^' _j - k_x ≡ ±l mod n_x and y^' ₂ - k₂ ≡ 0 mod ?z₂ ory'_j - C_j ≡ OmodzZ_j and ^' ₂ -/c₂ ≡ ±l mod?z₂ if y^' _j —k_x ≡ ±2 mod H_j and y₂ - k₂ ≡ 0 mod «₂ ory^' _j —k_x ≡ 0mod »_j and y₂ - c₂ ≡ +2 mod«₂ if y^' _j -k ≡ +3 mod «i and ^' ₂ -k₂ ≡ 0 mod«₂ or y^' _j - k_x ≡ 0 mod n and y^' ₂ - k₂ ≡ ±3 mod n₂

otherwise or by the array,

The constants α, β, γ and δ are the autocorrelation coefficients between the 2-D spline function R(m ,m₂) . In the following analysis, τ and m are assumed to be integers and R(t) is the 1-D cubic-spline function, defined in (1). By the use of (33), one obtains constants a, β, γ and δ as follows:

1) Calculation of α: By FIG. 7 for τ = 3, the value of α is the sum of the overlapping values of the 2-D spline function R(m_x,m₂) , given by a = l + 4x ∑R²(tn/τ) (34) m=\

2) Calculation of β: By FIG. 8 for τ - 3, the value of β is the sum of the overlapping values of the 2-D spline function R(m_x,m₂) and the shifted 2-D spline function R(m + τ,m₂) , given by ? = ∑R(/n/r)R((z— m)/r)+ 2x ∑R((r + m)/r)R(/w/"r). (35) m=\ «ι=0

3) Calculation of γ: By FIG. 9 for τ = 3, the value of γ is the sum of the overlapping values of the 2-D spline function R(m_x,m₂) and the shifted 2-D spline function

R(m + 2τ,m₂) , given by r = ∑R(mlτ)R((2τ-m)lτ). (36) m=0

4) Calculation of δ: By FIG. 10 for τ = 3, the value of δ is the sum of the overlapping values of the 2-D spline function R(m_x,m₂) and the sliifted 2-D spline function R(m_x + 3τ,m₂) , given by δ = ∑R({τ + m)/τ)R((2τ-m)/τ). (37)

»ι=0

Specifically for τ = 3 the above equations yield the following parameters for a, β, γ

2x3 and δ: α = l + 4x ∑R²(rø/3) = 3.8916 ,

3-1 β = ∑R{m/3)R((3-m)/3)+ 2x ∑R((3 + w)/3)R(m/ 3) - 0.3786,

2x3 3 = ∑R(m/3)#((2x3 -m)/3) = -0.1070 , = ∑R((3 + m)/3)β((2x 3 -7?z)/3)= 0.0055 . m=0 m=0

FIG. 7 shows the overlap map grid points for calculating for τ=3, FIG. 8 shows the overlap map grid points for calculating β for τ=3, FIG. 9 shows the overlap map grid points for calculating γ for τ=3, and FIG. 10 shows the overlap map grid points for calculating δ for τ=3.

VI. A NOVEL FCSI ALGORITHM IMPLEMENTED BY WINOGRAD DFT AND

OVERLAP-SAVE METHOD

The Winograd DFT algorithm described in Dean P. Kolba and Thomas W. Parks, "A prime factor FFT algorithm using high-speed convolution," IEEE Trans, on Acoustic, Speech, and Signal Processing, vol. ASSP-25, No.4, pρ.281-294, Aug. 1977, [16]; S.

Winograd, "On computing the discrete Fourier transform," Mathematics of Computation, vol. 32, No.141, pp.175-199, Jan. 1978, [17]; and I. S. Reed, T. K. Truong, R. L. Miller, and B. Benjauthrit, "Further results on fast transforms for decoding Reed-Solomon codes over GE(2") for n = 4, 5, 6, 8," in the Deep Space Network Progress Report 42-50, Jet Propulsion Laboratory, Pasadena, CA, pp.132- 155, Jan. and Feb. 1979. [18], the contents of which are hereby incoφorated by reference, with a new type of overlap-save method is used in this section to implement the FCSI scheme for τ= 3.

Consider an image of size NxN = 512x512 pixels. If the FCSI scheme for τ= 3 is used to compress this original image, the size of the compressed image is reduced to [512 / 3^~|χ [^~512 / 3^"| = 171x171 pixels, where [x ] denotes the least integer greater than or equal to x. Since 171 is not a power of two, the 2-D FFT can not be used to solve for the X_kιJi2 in (30). To overcome this problem, one possible method is to expand the data Y_{j h} in (25) and B_{S >S} in (28), where O ≤y',-,5,. ≤ n, -I for i = i,2 from 171x171 pixels to

256x256 pixels, by appending zeros to the edges of the original reduced image. Then one would take the 2-D FFT of Y_{h h} and B_{S tS} to obtain 7_m>„ and B_{m n} , respectively. This would yield X_{m n} = 7„, „ IB_{m n} . Finally, one would take the inverse 2-D FFT of the X_{m n} to solve for the X_{k Aι} , where 0 < k_x,k₂ ≤ 255 , the compressed data to be transmitted or stored.

The disadvantage of using such a 2-D FFT method to compute the CSI scheme for τ= 3 is that it requires more computational time due to the increase in the size of the data from

171x171 pixels to 256x256 pixels.

Since the size of the compressed data to be transmitted for τ.= 3 is not always equal to a power of two, in order to compress such an odd-sized image, a 9-point Winograd DFT

(WDFT) described below is used to implement the FCSI scheme. First since 171 = 9x19 the Y_{j j2} and B_S , where 0 ≤ j_x,j₂,s_x,s₂ ≤ 110 , can be divided into 19x19 sub-images, each of size 9 9 pixels, shown in FIG. 11. One might think that the 2-D Fourier transform of each sub-image could be implemented directly by the use of the 9x9 WDFT algorithm, and that the 2-D Fourier transforms Y_{m n} and B_{m n} of the Y_{J ι} and B_{S >S} of size

17X171 could be obtained, respectively. In addition, the X_{m n} could be calculated from the 7_OTj„ divided by the B_{m n} .

Finally, one might expect that the inverse 2-D Fourier transform of X_{m n} could be obtained by the use of the 9x9 WDFT algorithm. In other words, the reconstructed image of size 17X171 might be obtained. However, this inverse transform is found experimentally to have serious artifacts. This is seen in the reconstructed image, shown in FIG. 12. FIG. 12 illustrates a reconstructed image with serious artifacts that used FCSI method implemented by the direct use of the 9 x 9 Winograd DFT for compression. FIG. 11 shows a 19^χ19 sub-images of size 9x9 in an exemplary image of size 171x171 To remove the artifacts found in FIG. 12 by the direct use of the 9x9 WDFT compression described above, a novel type of overlap-save sub-image technique is applied to the FCSI method. An illustrative example of the FCSI algorithm is implemented first by the simplified 5x5 WDFT described below. This new type of overlap-save sub-image method is illustrated in this example, the block diagram of which is depicted in FIG. 13. FIG. 13 is a simple illustrative example of the FCSI algorithm implemented by the

5^χ5 Winograd DFT and the overlap-save sub-image method. The block diagram shown in FIG. 13 is separated into two parts as indicated by broken lines. The first part, labeled by "I", is a FCSI encoder that uses the 5x5 WDFT algorithm with the overlap-save sub- image technique. The second part, labeled by "II", is a FCSI decoder. In the first part, consider the source image data to be of size 24x24 pixels as shown in FIG. 13(a). The use of (25) with zonal filter 4 for τ= 3 yields the reduced image needed to obtain the coefficients Y_h of size 8x8 pixels as depicted in FIG. 13(b). Next these coefficients

Y _h are divided into four overlapping 5 5 sub-images as shown in FIG. 13(c). Note that each sub-image of size 5x5 pixels overlaps each adjacent sub-image with a border of width 2. It is shown experimentally that this border can be used to solve the boundary condition between the two neighboring sub-images when using the Winograd DFT algorithm.

As illustrated in FIG. 13, one uses the next three steps to transform the four I overlapping 5^χ5 sub-images of 7, ₂ to the corresponding four overlapping 5x5 sub- images of _{k] ι} . The first step is to take the 5x5 WDFT of the four overlapping 5x5 sub- images of Y_{j J2} to obtain the four transformed 5x5 corresponding sub-images of Ϋ_{m n} .

5 The second step is to divide these four 5 5 sub-images of Y_{m n} by the 5x5 WDFT of B_{s s} in (33) or B_m „ for τ= 3 to obtain the four transforms of the 5x5 sub-images of x_m>n . The third step is to take the inverse 5x5 WDFT for these four 5x5 sub-images of X_{m n} to obtain finally the corresponding four overlapping 5x5 sub-images of X_kιιkι as shown in

FIG. 13(d). 10 Because some of the pixels in the overlapping border of 5x5 sub-images of X_kιΛι appear in the other adjacent 5 5 sub-images of X_{k tkl} , the duplicated pixels in the four overlapping 5x5 sub-images of χ_{k >kι} are deleted or removed. By this means the four overlapping 5x5 sub-images of X_{k Jc} become the four non-overlapping 4x4 sub-images of

X_k .₂ . To illustrate this, in FIG. 13(d) each 5 5 sub-image has an overlapping border of

15 width 2; FIG. 13(e) illustrates the remaining samples of each sub-image obtained by the use of this overlap-save method. In FIG. 13(d), first consider all four sub-images with overlapping border in the column direction. Because of the end-effects, the last columns of sub-images 1 and 3 are the duplicated columns in the overlapping border to be removed. 0 However, in sub-images 2 and 4 in FIG. 13(d), the first columns of these two sub- images are also duplicated columns in the overlapping border that need to be deleted. Finally, the above overlap-save method in the row direction is accomplished similarly to that of the column direction. A combination of these four non-overlapping 4x4 sub- images of X_{k Ji2} yields the entire 8x8 image of X_{k >kϊ} shown in FIG. 13(f). These X_kukι 5 image data are the compressed data to be transmitted or stored. In the second part, using the entire 8x8 X_{k J} compressed data, the 24x24 reconstructed data, shown in FIG. 13(g), are obtained by means of the cubic-spline reconstruction function, given in (32).

Thus the steps included in compressing an image are defining a cubic-spline filter; correlating the filter with the signal to obtain a correlated signal; autocorrelating the filter 0 to obtain autocorrelated filter coefficients; computing a transform of the correlated signal and the autocorrelated filter coefficients; dividing the transform of the correlated signal by the transform of the autocorrelated filter coefficients to obtain a transform of a compressed signal; and computing an inverse transform of the transform of the compressed signal to obtain the compressed signal. ^ In order to compress an actual image of size 512x512 by such an overlap-save sub- image method, the 9x9 WDFT instead of the 5x5 WDFT is used for the FCSI scheme. Then, the FCSI encoding algorithm using the 9x9 Winograd DFT and the overlap-save sub-image method is summarized in the following steps: • Choose τ= 3. The compression ratio is roughly τ² = 9.

• Apply (25) with zonal filter 4 to the original image of size 512x512 to find all of the 171x171 coefficients, 7, , . • Also apply (33) to find B_{s t} . Then take the 9x9 WDFT of B_{S S} to obtain the 9x9 coefficients, B_{m n} .

• Divide all of the coefficients 7 ■ into the appropriate overlapping 9^χ9 sub-images of Y: : with a border of width 2.

• Take the 9x9 WDFT of all overlapping 9x9 sub-images of Y_jχJι to obtain the transforms of the corresponding 9x9 sub-images of _m „ .

• Compute X_{m n} = Ϋ_{m n} IB_{m n} to obtain each 9x9 sub-image of X„_hn .

• Take the inverse 9x9 WDFT of all of the 9x9 sub-images of X_mjt to obtain the 9 9 sub-images of X_k^ .

• Remove the duplicate pixels in the overlapping border of two adjacent 9x9 sub-images of X_{k kϊ} . By this means the overlapping sub-images of X_{k kι} become the non- overlapping sub-images of X_k ^ .

• Combine each non-overlapping sub-image of X_{k k} to obtain the entire image of X_k fo , where 0 ≤ k_x,k₂ < 511. These image data X^ are the compressed data to be transmitted or stored. In the FCSI decoding algorithm, the compressed data X_{k k} constitutes the encoding algorithm; the reconstructed data are obtained by means of the cubic-spline reconstruction function in (32). Hence, by means of the above-described overlap-save sub-image technique, both the compressed and the reconstructed image are obtained and found experimentally to have no apparent artifacts, as shown in FIG. 14. FIG. 14 is a reconstructed image with no apparent artifacts that used the FCSI implemented by the 9x9 Winograd DFT and the overlap-save method for compression. VII. A MODIFIED JPEG ENCODER-DECODER

In this section, a modified JPEG encoder-decoder is presented for image data compression. This algorithm applies the CSI or FCSI scheme with a compression ratio of r² to 1 as the pre-processing step of the JPEG encoder for τ= 2, 3 as shown in FIG. 15. As a consequence the cubic-spline reconstruction with a ratio of 1 to τ² is used for the post-processing step of this JPEG decoder for τ= 2, 3 as shown in FIG. 16. For this compression algorithm, an original image in the RGB (Red, Green and Blue) color space is converted into another preliminary image in YUN color space prior to the CSI or FCSI pre-processing.

This YUV image is followed by the CCIR 601 color space with format 4:1:1. The size of the original RGB image is assumed to be 512x512x3 =786,432 bytes, i.e., each set of 512x512 = 262,144 bytes is for one of the red, green, and blue colors. After color- space conversion one set of 512x512 bytes is used for Y, and two sets of 256x256=65,536 bytes are used for the U and V color components. It follows from Phillip E. Mattison, Practical Digital Video with Programming Examples in C, John Wiley & Sons, Inc., 1994 [19], the contents of which are hereby incoφorated be reference, that the equations for the conversion from RGB to YUV are Y = 0.299R + 0.587C7 + 0.1145, U = 0.493(5 - Y) = 0.463E - 0.147R - 0.289G and V= 0.877(R - Y) = 0.615R - 0.515 - 0.1005.

There are two processing steps in the encoder phase. The first step is the preprocessing that uses the CSI or FCSI scheme with a compression ratio of τ² to 1 for each of Y, U, and V images. In this procedure, the input image is a Y image of size 512 512 bytes, and the output image is an encoded image of size [^~512/r^~|χ [512/r^"j bytes, where denotes the least integer greater than or equal to x. For the U and V images the input image has 256x 256 bytes so that the output image to be encoded is 256/r χ [^"256/τ^"| bytes. In other words, there are two cases used in this step. For the first case τ= 2, the CSI scheme implemented directly with the FFT algorithm is used for the original Y, U and V images and the output images are 256x256 bytes for Y image and 128x128 bytes for U and V images.

For the second case τ= 3, the FCSI scheme implemented by the 9x9 WDFT with overlap-save method is used for the original Y, U and V images and the output images are 17X171 bytes for Y image and 85x85 bytes for U and V images. At the end of the CSI or FCSI algorithm the three separate Y, U, and V images are combined into one YUV image. The second step is to use the JPEG DCT-based encoding algorithm [9]. The image after this step is called the compressed image. This compressed image has now a very small number of pixels when compared to the original image. The resulting image still has the standard JPEG format. As a consequence, this compressed image can use the standard JPEG decoder, also save on storage and decrease the transfer time for a communication.

Also in the modified JPEG decoder, there are two processes used which are reversed in some of the encoding steps. The first step is the JPEG DCT-based decoding algorithm [9]. After this step, the image file is separated into three separate Y, U, and V images. The second step is the post-processing step that uses the cubic-spline reconstruction with a ratio of 1 to τ² for Y, U, and V image. This step uses only the cubic-spline function to reconstruct the image data. After this inteφolation the size of the Y image is therefore converted from f512/τ]^χ f512/r], i.e., 256x256 for τ= 2 or 171x171 for τ= 3, to 512x 512 bytes, and the U and V images are increased from f^"256/τ]χ[256/r , i.e., 128x128 for τ= 2 or 85x85 for τ= 3, to 256x 256 bytes. Then the three Y, U, and V images are combined again into one YUV format. Finally, this YUV image is converted into the reconstructed RGB image. Again, it follows from [19] that the equations for the conversion from YUN to RGB ar R = Y+ IΛAOV, G = Y- 0.395U- 0.581 F, and 5 = 7+

2.032L .

VIII. EXPERIMENTAL RESULTS

Let X(i,j) and S(i,j) be the original and reconstructed images, respectively, where 0 < i < M - 1 and 0 ≤ j ≤ N-l are the index numbers separated in the vertical and horizontal directions of the images. The mean-square error (MSE) of the 2-D signal of an image is given by

Thus, the PSNR of the 2-D signal are defined by

where MSE- γ > MSE_V and MSE_V are the MSE of color component Y, U and V, respectively. The total PSNR of the color component Y, U and JXs defined by

255'

PSNR_T(dB) = 10log_l0 - (40)

(MSE_Y + MSE_υ + MSE_V ) 13 Experimental results for the 2-D signal are presented using linear inteφolation, linear spline inteφolation, cubic-convolution inteφolation, cubic B-spline inteφolation and the CSI scheme. Using the encoding method described in Section II and the decoding method described in Section III for the 2-D signal with x= 2, from (38) and (39), the PSNR values of the 2-D signal with a compression ratio of 4:1 (τ = 2 ) are calculated for some gray images of size 512 by 512. The experimental results of the 2-D signal using the above five inteφolation schemes are shown in Table III. The results in this table show that the CSI scheme obtains the best PSNR over all other compared inteφolation methods.

TABLE III THE PSNR (dB) OF 2-D FUNCTION WITH COMPRESSION RATIO 4 : 1 ( τ = 2 ).

In Table IV, the experimental results with a compression ratio of 9:1 (τ - 3) are presented for the FCSI algorithm that uses zonal filter 1, zonal filter 2, zonal filter 3 and zonal filter 4 for the same gray images of size 512 by 512. It is observed from this table that the FCSI scheme with zonal filter 4 obtains a PSNR that is similar to the other three zonal filters. In addition, the experimental results with a compression ratio of 9:1 (r = 3) that use linear inteφolation, linear-spline inteφolation, cubic-convolution inteφolation and the FCSI scheme with zonal filter 4 are shown in Table V for the same gray images of size 512 by 512. It follows from this table that the FCSI scheme with zonal filter 4 obtains the best PSNR of the four inteφolation methods.

TABLE IV THE PSNR (dB) WITH COMPRESSION RATIO OF 9: 1 ( τ ~ 3 ) FOR THE FCSI SCHEME WITH

DIFFERENT ZONAL FILTERS.

TABLE V

THE PSNR (dB) WITH COMPRESSION RATIO OF 9: 1 ( τ = 3 ) FOR FOUR INTERPOLATION

SCHEMES.

Table VI lists the PSNR values of the reconstructed color test (Lena) image of size 512 by 512 at different compression ratios for the JPEG method [13], the AOL algorithm [8], and the modified JPEG encoder-decoder for τ= 2 described in FIG. 15 and FIG. 16, also see [10]. For the same compression ratios the PSNR of the Lena image, obtained by the modified JPEG encoder-decoder for τ= 2 are higher than those of both the JPEG method and the AOL algorithm. Furthermore, using (40), Table VII lists the total PSNR values of the colored reconstructed Lena image of size 512 by 512 at higher compression ratios for the JPEG method, the AOL algorithm, and the modified JPEG encoder-decoder for τ= 2 and τ= 3. For the same compression ratio of 250:1 the total PSNR value of the Lena image, obtained by the modified JPEG encoder-decoder for τ= 3, is higher than the JPEG method, the AOL algorithm and the modified JPEG encoder-decoder for τ.= 2.

TABLE VI PSNR (dB) VALUES OF THE RECONSTRUCTED COLOR LENA IMAGE OF SIZE 512 BY 512 AT DIFFERENT COMPRESSION RATIOS FOR JPEG, AOL, AND MODIFIED JPEG ENCODER- DECODER FOR τ= 2.

TABLE Vπ TOTAL PSNR (dB) VALUES OF THE RECONSTRUCTED COLOR LENA IMAGE OF SIZE 512 BY 512 AT DIFFERENT COMPRESSION RATIOS FOR JPEG, AOL, MODIFIED JPEG ENCODER- DECODER FOR τ= 2 AND τ= 3.

For the gray Lena image of size 512 by 512, the computational times of the CSI and FCSI schemes for τ= 3 were implemented on a 400-MHz Intel Pentium II personal computer using C-code. In the encoder, the FCSI scheme using the 9-point WDFT with overlap-save sub-images requires around 0.15 sec when compared with around 0.57 sec for the CSI using the FFT. Therefore, the FCSI scheme is faster than that of the CSI scheme.

Finally, the JPEG method, the AOL algorithm, and the modified JPEG encoder- decoder for τ= 2 and τ= 3 were also implemented on the same 400-MHz Intel Pentium II personal computer using C-code. The computational time of the color Lena image of size 512 by 512 at the compression ratio of 200:1 for these four algorithms are given in Table VIII. In the encoding and decoding, the modified JPEG encoder-decoder for τ= 3 requires only 0.42 sec and 0.27 sec when compared with 1.13 sec and 0.34 sec for the modified JPEG encoder-decoder for τ= 2, 1.09 sec and 0.30 sec for the AOL algorithm, and 0.80 sec and 0.65 sec for the JPEG method, respectively. It is by this means the computational time of the modified JPEG encoder for τ= 3 requires, respectively, 0.71 sec, 0.67 sec, and 0.38 sec less time than the modified JPEG encoder for τ= 2, the AOL encoder, and the JPEG encoder. Also, the computational time of the modified JPEG decoder for τ= 3 requires, respectively, 0.07 sec, 0.03 sec, and 0.38 sec less time than the modified JPEG decoder for τ= 2, the AOL decoder, and the JPEG decoder.

TABLE NIII THE COMPUTATIONAL TIME (sec) OF THE COLOR LENA IMAGE OF SIZE 512 BY 512 AT THE COMPRESSION RATIO OF 200:1 IMPLEMENTED ON A 400-MHz INTEL PENTIUM II PERSONAL

COMPUTER.

FIG. 17 shows the reconstructed image of Lena at the same compression ratio of 100:1, using the JPEG method, the AOL algorithm, and the modified JPEG encoder- decoder for τ- 2. Clearly the Lena image using the modified JPEG encoder-decoder for τ= 2 indicates a better subjective quality of reconstructed image than both the JPEG method and the AOL algorithm. Also FIG. 18 shows the reconstructed image of Lena at a higher compression ratio of 200:1, using the JPEG method, the AOL algorithm, and the modified JPEG encoder-decoder for τ= 2 and τ= 3. In this figure, the Lena image using the modified JPEG encoder-decoder for τ= 3 indicates a better subjective quality of reconstructed images than the JPEG algorithm, the AOL algorithm and the modified JPEG encoder-decoder for τ= 2.

FIG. 17(a) shows an original Lena image, while FIG. 17(b) depicts a reconstructed image by the JPEG method with PSNRγ= 30.76 dB, PSNRrj = 33.90 dB, and PSNRγ = 33.90 dB. FIG. 17(c) Illustrates a reconstructed image by the AOL algorithm with PSNRy = 30.91 dB, PSNRj = 35.60 dB, and PSNRr = 35.75 dB; and FIG. 17(d) represents a reconstructed image by the modified JPEG encoder-decoder for τ- 2 , PSNRy = 31.20 dB, PSNRJ = 35.66 dB, PSNRy^ 35.85 dB.

FIG. 18(a) shows a reconstructed image by the JPEG algorithm with encoding time = 0.80 sec, decoding time = 0.65 sec, and PSNRf = 27.90 dB, while FIG. 18(b) depicts a reconstructed image by the AOL algorithm with encoding time = 1.09 sec, decoding time = 0.30 sec, and PSNR = 31.06 dB. FIG. 18(c) illustrates a reconstructed image by the modified JPEG encoder-decoder for x~2, encoding time =1.13 sec, decoding time = 0.34 sec, and PSNR ⁷ = 31.18 dB; and FIG. 18(d) represents a reconstructed image by the modified JPEG encoder-decoder for τ= 3, Encoding time = 0.42 sec, Decoding time = 0.27 sec, PSNR = 31.19 dB.

9-POINT WINOGRAD DISCRETE FOURIER TRANSFORM

The algorithm developed in Dean P. Kolba and Thomas W. Parks, "A prime factor FFT algorithm using high-speed convolution," IEEE Trans, on Acoustic, Speech, and Signal Processing, vol. ASSP-25, No.4, pp.281-294, Aug. 1977 [16], the contents of which are hereby incoφorated by reference, for computing the 9-point Winograd discrete Fourier transform (DFT) is shown by

X(k)=∑x(n)W^nk, A = 0,-,8, (41) n=0 where W = _e ^~j(-^2π/9^ is the 9-th root of unity in the complex number field, and j = ΛXT . Algorithm for the 9-point Winograd DFT : a_λ = x(l) + x(8) , a₂ = x(l) - x($) , a₃ = x(2) + x(l) , α₄ = x(2) - x(l) , a₅ = x(4) + x(5) , ^a6 ^{= x}(4) ^_ *(5) , a₇ = x(3) + x(6) , α₈ =x(3)-x(6),

a_x j = — α₃ + αj , a_X2 ^{= a}2 ~ ^a6 ■> ^α13 ^{= a}2 ^{+ α}43 ^Ω14 ^{= ~a}4 ^{~~ a}6 ' ^a\S ^{~ a}\ ^{+ α}3 ^"*^{" a}5 > ^a6 ^{= a}2 ^{~ a}4 ^{+ a}6 J ^a\l ^{~ X}(Q) + «15 + aη , '«ι =0.19740α₉, m₂ = 0.56858₁₀ , m₃ =0.37111 _n, m₄ = 0.54253α₁₂ , m₅ =0.10026α₁₃, ma = 0.44228α₁₄ , mη =-^ ₇, w₈ =0.86603α , trig =^a_X5, m_XQ = 0.86603α₁₆ , c = x(0) - niη , c₂ = m - ₃ , c₃ = m_x + ?w₃ , c =m_x+ rø , c₅ =c-_[+C2-c₃, eg = + c₃ + c ,

c_u=c_&+ c_x0 - ₈ , c₁₃ = ~c₉ + c₁₀ + ₈ , c₁₄ = x(0) + ₇ - τn₉ , X(0) = a_xl, (l) = c₅-y^'c_π, (2) = c₆ -y^'c₁₂ , X(3) = c_u - jm_l0 , X(4) = c₇-jc_u,

X(5) = c_η + yc₁₃ , X(6) = c₁₄ + y^'m₁₀ , X(7) = c₆+ jc_l2 , X($) = c₅+ jc_{x x} .

Thus, the 9-point Winograd DFT requires only 8 multiplications, 49 additions, and 2 shifts, a substantially smaller number of computations than other known algorithms. THE 5-POLNT WINOGRAD DISCRETE FOURIER TRANSFORM The algorithm developed in [16] for computing the 5 -point Winograd discrete Fourier transform (DFT) is given by:

X(k) = ∑x(n)W"^k, /c = 0,-,4, (42)

«=o where W = _e ^'Λ2π/5) is the 5-th root of unity in the complex number field, and j = -ϊ ■ Algorithm for the 5-point Winograd DFT is given by: α, = x(l) + x(4) , a₂ = x(l) - x(4) , a₃ = x(2) + x(3) , a₄ = x(2) - x(3) , a₅=a₂+a₄, a₆=a_x-a₃, a_η = a_x + a₃ , a₈ = x(0) + α₇ , m_x =0.95106α₅, m₂ =1.53884α₂, m₃ =0.36327₄, m₄ = 0.55902α₆ , m₅ =}.,, c_x = x(0) — m₅ , c₂ = C_j + m₄ , c₃ = c - m₄ , c₄ = m_x - m₃ , c₅ = m₂ - m_x ,

X(0) = a_s , X(\) = c₂ - jc₄ , X(2) = c₃ - jc₅ , X(3) = c₃ + jc₅ , X(4) = c₂ + jc₄ .

Thus, the 5 -point Winograd DFT requires only 4 multiplications, 17 additions, and 1 shift, a substantially smaller number of computations than other known algorithms.

In this invention, a new CSI scheme based on the least-squares method with the cubic-spline function has been proposed to compress the image data. It is shown that the CSI scheme implemented by the FFT algorithm yields a better PSNR performance than all other inteφolation methods for the reconstructed image. In addition a fast CSI, called FCSI, is developed for image compression. Such a FCSI scheme requires fewer additions and multiplications in the decimation process than the original CSI scheme. In the FCSI, a fast 9-point WDFT algorithm is used to compute the CSI scheme with zonal filter 4 and an overlap-save sub-image technique is utilized to solve the serious boundary artifacts between any two neighboring sub-images of the actual image. It is shown by computer runs that the FCSI scheme requires 0.42 sec less time than the original CSI scheme for τ= 3.

Finally, this FCSI scheme for τ= 3 is used with the JPEG encoder-decoder to speed up the modified JPEG encoder-decoder for τ= 2 in color image coding. A computer simulation shows that the modified JPEG encoder-decoder for τ= 3 obtains a better subjective quality and PSNR of the reconstructed image than the JPEG method for high compression ratios. Also, it requires less computational time than the JPEG method, the AOL algorithm and the modified JPEG encoder-decoder for τ= 2 in both the encoding and decoding processes.

It is understood that the exemplary schemes and the respective implementations described herein and shown in the drawings represents only exemplary embodiments of the present invention. Indeed, various modifications and additions may be made to such embodiment without departing from the spirit and scope of the invention. For example, the present invention can be implemented utilizing a computer program, a dedicated electronic circuit, or firmware for an electronic image processing chip. Furthermore, those skilled in the art will appreciate that the method of the present invention is applicable to images having various different sizes with different values for τ and different sizes for Winograd DFT. Thus, the description herein of a 9 x 9, or 5 x 5 of Winograd^'DFT is by way of illustration only, and not by way of limitation. Further, those skilled in the art will appreciate that various different VLSI implementations of the present invention are contemplated. Thus, these and other modifications and additions may be obvious to those skilled in the art and may be implemented to adapt the present invention for use in a variety of different applications including streaming of video signals, combining with MPEG IV standard for still pictures, applying to JPEG 2000 standard, reducing the size of digital pictures (such as in digital cameras), and the like.

Claims

WHAT IS CLAIMED IS:

1. A method performed by a computer for coding a 1-D signal comprising the steps of: defining a 1-D cubic-spline filter by

(3/2)|t|³ -(5/2)|t|² ₊ l ,0<|t|<l R(t)H -(l/2)|t|³+(5/2)|t|²-4|tj + 2 ,l<|t|<2 (1) ,2<|t| applying the filter to an input signal Xm with Y_j= ∑X(m + jτ)R(m) (9)

»ι=-2r+l to compute Y_j : applying

2r-l

A _k= ∑ R(m + (j-k)τ)R(m) m=-2r+l

(10)

where A_jk can be expressed in circulant form by letting

^Ai,k - ^B(k-j)„ (11) and where (k - j)_n denotes the residue (k - j) modulo n , and B₀= ,B_x=;β,B₂=_r,B₃=δ,B₄=0, B_n-₄=0,B_n.₃ =δ,B_n_₂=_r, (12)

B_n__x = β , to compute _9₇ ; computing, the FFT of Y. and B_j to obtain Ϋ_m and B_m , respectively; computing X_m = YjB_m ; and computing the inverse FFT of X_m to obtain a compressed data X_k .

2. The method of claim 1 further comprising the step of:

«-l applying X_k in ∑ A_J>kX_k = 7₇ for 0≤j≤n-l , (5) to perform a convolution of

(3/2)|r -(5/2)|t|²+l ,0≤|t|<l

X_k and R(t) given by R(t) = (l/2)|t|³+(5/2)|t|²-4|t| + 2 ,l≤|t|<2 (1) 0 ,2<)t|

«-l to obtain S(t) given by S(t) = ∑X_kΨ_k (t) = ∑X^_kR(t-kτ) (3). i=0 4 Σ=0

3. The method of claim 1 , wherein x= 2.

4. The method of claim 1 , wherein τ= 3.

5. An encoder for data compression comprising: a 1-D cubic-spline filter defined by

(3/2) |t|³ -(5/2)|t|² + ι ,o≤|t|<ι

Λ( = -(l/2)|t]³+(5/2)|t|²-4|t| + 2 ,l<|t|<2 (1)

,2<t means for applying the filter to an input signal Xm with

2r-l

Y_j = ∑ X(m + jτ)R(m) (9) to compute 7_y : m=-2r+l means for applying

2τ-l

= ∑ R(m + (j-k)τ)R(m) m=-2τ+l τ² -^x _2r ^R{m) = if j-k≡ Omodw

∑XX₊i ^R(^{m + τ})^R(m) = β if j - k ≡ ±1 mod n (10)

∑²X-_2r+ι ^R(^m + 2r)Λ(m) = γ if j - k ≡ ±2 mod n

∑_m=-_2r+\ ^R(^m + 3τ)R(m) = δ if y - k ≡ ±3 mod n

0 otherwise where A_Jζ can be expressed in circulant form by letting

^Aj,k ~ ^B(.k-j)„ J (ID and where (k - j)_n denotes the residue (k - j) modulo n , and

-_t B_n_₄ = 0 , B„_₃ = δ , B_n_₂ = γ , (12) B„__x = , to compute B_j ; means for computing the FFT of 7- and B_j to obtain Y_m and B_m respectively; means for computing X_m - YJB_m ; and means for computing the inverse FFT of X_m to obtain a compressed data X_k

6. The encoder of claim 5 further comprising a decoder having means for

«-l applying X_k in ∑A_jkX_k =Y. for O≤j≤n-l , (5) to perform a convolution of l

X_k and R(t) given by R(t) 2 (1)

to obtain S(t) given by S(t) = ∑ X_kΨ_k (t) = ∑ X_kR(t - kτ) (3) k=0 k=0

7. The encoder of claim 5, wherein τ= 2.

8. The encoder of claim 5, wherein τ= 3.

A method performed by a computer for coding a 2-D signal comprising the steps of: defining a2-D cubic-spline filter by R(t,t₂) = R(t_x)-R(t₂), (17) where R(t_x) and R(t₂) are 1-D cubic-spline functions, respectively; applying the filter to an input signal (t_j,t₂) with periods nτ and n₂τ , with

2r-l 2r-l

7 Λ,.Λ = Σ Σ X(m_x+j_xτ,m₂+j₂τ)R(m_x,m₂), (25) m,=-2τ+l m,=-2r+l to compute 7_{7 h}

2τ-\ 2τ-\ applying A_JJιAkι = Σ ∑ R(m + (j_x - k )τ, m₂ + (j₂ - k₂ )τ)R(m_x , m₂ )

»ι,=-2r+l ₂=-2r+l if y^' _j - k_x ≡ 0 mod n_x and y^' ₂ - k₂ ≡ 0 mod «₂ if y^'ι —k ≡±lmodz._j andy^' ₂ —k₂ ≡0mod;;₂ or y'_j - k ≡ 0 mod n_x and j₂-k₂ ≡±l mod n₂ if y^'i -k_x ≡±2modκ_j andy^' ₂ -k₂ ≡0 od/J₂ or y^' _j - k_x ≡ 0 mod n_x and y^' ₂ —k₂ ≡ ±2 mod n₂ if j_x —k_x ≡±3mod«[ and ^' ₂ -A₂ ≡0mod«₂ or ^' _j -k_x ≡Omodn_j andy₂ — k₂ ≡±3modn₂ if y^' _j -k ≡±lmodH_j andy^' ₂ —k₂ ≡±l mod n₂ if y^' _j —k ≡±2mod«j andy^' ₂ —k₂ ≡±lmod7z₂ ory^' _j -k_x ≡±lmod/Z_j andy^' ₂ -k₂ ≡±2modzz₂ if J, -k_x ≡±3mod72_j andy^' ₂ -k₂ ≡±lmodzz₂ ory^' _j —k_x ≡±lmodzz_j andy^' ₂ —k₂ ≡±3modzι₂ if y^' _j - k_x ≡ ±2 mod zz_j and y^' ₂ - k₂ ≡ ±2 mod zz₂ if y^' _j - k_x ≡ ±3 mod zz_j and y^' ₂ - k₂ ≡ ±2 mod zz₂ ory^' _j —k_x ≡±2mod/Z_j andy^' ₂ —k₂ ≡ ±3 mod n₂

otherwise where, the array A_{hh^k k} can be expressed in 2-D circulant form as

^ΛjJιMϊ ^=B{ -J ,_h, ⁽*2-Λ)„₂ ' (²⁷) and where (k_t-j_t) denotes the residue (k_t -j,) modulo n_t for t = l,2 and J

where 0 < s_; ≤ n, - 1 for i - 1,2 , to compute 5 computing the 2-D FFT of Y _h and B_{S Λ} to obtain Y_mn and B_mn , respectively; computing „_Vi = Y_mn /B_m>n ; and computing the inverse FFT of X_m „ to obtain a compressed image X_{k kl} ■

10. The method of claim 9 further comprising the step of:

Applying X_kχ in ∑ " ,_Λ^ _2,^₂ =Y_λιh , ≤j ≤n, -1 and ι = l,2. (21) i,=0 A₂=0 to perform a 2-D convolution of X_kι and R(t_x,t₂) given by R(t_x,t₂) = R(t_x ) ■ R(t₂ ) to obtain a 2-D reconstructed image S(t_j , t₂ ) .

11. The method of claim 9, wherein σ = 0.

12. The method of claim 9, wherein σ = μ = 0.

13. The method of claim 9, wherein η = λ = p = μ = σ = 0.

14. The method of claim 9, wherein τ= 2.

15. The method of claim 9, wherein τ= 3.

16. An encoder for data compression comprising: a 2-D cubic-spline filter defined by R(t_x,t₂) = R(t_x)-R(t₂) , (17) where R(t_x) and R(t₂) are 1-D cubic-spline functions, respectively; means for applying the filter to an input signal X(t_x,t₂) with periods nτ and

2r-l 2r-I with ^YJ,J, = Σ Σ ^χ-"h + 7X ^m2 + J2^T)^R(^tnι . ^m2 ) , (25)

-2r+l m₂=-2r+l to obtain 7 J,lJl ' means for applying

2r-l 2τ-l

A JJ2M2 = Σ Σ R(m_x+(j -k_x)τ,m₂+(j₂-k₂)τ)R(m_x,m₂) m,=-2r+l m,=-2r+l

where, the array A, . _{k k} can be expressed in 2-D circulant form as

A ^jJi ⁽*₂-Λ⁾»₂ ' (27)

and where (k_t -j.) denotes the residue (k_t - y^',) modulo n , for z^' = l,2 and

where 0 < s, < n_; - 1 for i = 1,2 , to obtain 5, , ; means for computing the 2-D FFT of 7₇ -₂ and 5 to obtain Y_{m n} and 7i_{m n}, respectively; means for computing X_{m n} = 7_ra „ /5„_M, ; and means for computing the inverse FFT of X_{m n} to obtain a compressed image

X

17. The encoder of claim 16 further comprising a decoder having means for applying X in l, 2. (21)

to perform a 2-D convolution of _{k k} and R(t_x,t₂) given by R(t, , t₂ ) = R(t ) ^■ R(t₂ ) to obtain a 2-D reconstructed image S(t_x ,t₂) .

18. The encoder of claim 16, wherein τ= 2.

19. The encoder of claim 16, wherein τ= 3.

20. The encoder of claim 16, wherein σ = 0.

21. The encoder of claim 16, wherein σ = μ = 0.

22. The encoder of claim 16, wherein η = λ = p = μ = σ = Q .

23. A method performed by a computer for coding an image X(t_x,t₂) comprising the steps of: applying a zonal filter R(m,m₂) to the image N(tj,t₂) where

2τ-\ 2r-l

^YJ,_J, ^~ ∑ ∑ X(m_x+j_xτ,m₂+j₂τ)R(m_x,m₂), (25) wi_[=-2r+l ₂=-2r+l to compute 7 Mi '

2τ-\ 2r-l applying A,J,2 2 ∑ ∑ R(mχ + (j, - K )τ, m₂ + (y₂ - k₂ )τ)R(m ,m₂) m,=-2r+l τιι₂=-2r+l

where, the array A j\jl ' \ 2 can be expressed in 2-D circulant form as

A = E

J lA^k2 (*ι-Λ)n„ (^k -J2)„ (27) and where (k_t- -y^',)„_. denotes the residue (k_t - j,) modulo 77 ,. for z^' = l,2 and

-#0,«₂-2 ^■°0,71_j-l -#l,«₂-2 B 1,71,-1 -#2,»₂-2 2,«₂-l ^■"3,71,-2 3,n,-l

3,n₂-2 -#71,-3,1₂-1 -"7l,-2,71₂-2 - ■"S7„I,-2,71₂-1

5 n, -l,«,-2 B «|— I./I₇-I

where 0 ≤s_i ≤n₍ -1 for z^' = l,2, to compute B ; where constants α, β, γ, δ, η, λ, p, μ, and σ are autocorrelation coefficients between zonal filter R(m_x,m₂) ; dividing all of the coefficients B into respective overlapping sub- images of B with a border of width 2; computing a Winograd discrete Fourier transform (WDFT) of all overlapping sub-images of B to obtain coefficients, B_{m n} ; dividing all of the coefficients Y_{j h} into respective overlapping sub-images of Y; i with a border of width 2; computing a WDFT of all overlapping sub-images of Y_{J ι} to obtain transforms of respective sub-images of Y_{m n} ; computing X_{m n} = Y_{m tl} I B_{m n} to obtain each sub-image of X_m„ ; computing an inverse WDFT of all of the sub-images of X_m„ to obtain sub-images of X_k ; removing duplicate pixels in an overlapping border of two adjacent sub- images of X_{k kι} ; and combining each non-overlapping sub-image of X_{k k} to obtain the entire image of X .

24. The method of claim 23, wherein τ= 2.

25. The method of claim 23, wherein τ= 3.

26. The method of claim 23, wherein the autocorrelation coefficient between zonal filter σ = 0.

27. The method of claim 23, wherein autocorrelation coefficients between zonal filter σ = μ = 0.

28. The method of claim 23, wherein autocorrelation coefficients between zonal filter η = Z = p = μ = σ = 0.

29. The method of claim 23, wherein the WDFT is a 9 x 9 WDFT given by

X(k)=∑x(n)W"^k, λ = 0,-,8, (41) and the sub-images of 7_{Λ 2} are 9 x 9.

«=o

30. The method of claim 23, wherein the WDFT is a 5 x 5 WDFT given by,

4

X(k) = ∑χ(ή) W"^k , k = 0, • ^■ • , 4 , (42) and the sub-images of 7_ΛΛ are 5 x 5.

31. The method of claim 23 further comprising the step of: applying X_k^ in ∑ ∑ X_kl,_k2 ^A _λj2,_kk2= _J\,_j2 > 0<j) ≤n_t -1 and i = l,2. (21) τt,=0 k₂=0 to perform a 2-D convolution of X_k>k and R(t,t₂) to obtain a 2-D reconstructed image S(t,,t₂).

32. A modified JPEG encoder for image compression comprising: means for applying a zonal filter R(m_x,m₂) to an image Z(tj,t₂) where

^Y _A,j2 = X(m_x+j_xτ,m₂±j₂τ)R(m,m₂), (25)

to compute 7, _h ; means for applying

^Λfj₂,k_xk₂ = + Ul - )^T> ^ftl2 + U2 - ^k2 )R(fl J ^m2 )

where, the array A, ■ _{k k} can be expressed in 2-D circulant form as

= B, (⁴ι^_Λ)ιι„ (^k2-J2 n₂ (27) and where (k_l - y^',),,_. denotes the residue (k_t - _ ,) modulo 77 ₍- for z^' = l,2 and

where 0 < s . < n, - 1 for i = 1,2 , to compute B_{s S2} ; where constants α, β, γ, δ, η, λ, p, μ, and σ are autocorrelation coefficients between zonal filter R(m_x,m₂) ; ι means for dividing all of the coefficients B into respective overlapping sub-images of B with a border of width 2; means for computing a Winograd discrete Fourier transform (WDFT) of all 5 overlapping sub-images of B to obtain coefficients, B_m „ ; means for dividing all of the coefficients 7,- ₇ into respective overlapping sub-images of Y_{j ι} with a border of width 2; means for computing a WDFT of all overlapping sub-images of 7,- -₂ to obtain transforms of respective sub-images of Ϋ_m>n ; 10 means for computing X_m _„ = 7_m „ /i?„ _! to obtain each sub-image of X_{m n} ; means for computing an inverse WDFT of all of the sub-images of X_m „ to obtain sub-images of X_{k >kι} ; means for removing duplicate pixels in an overlapping border of two adjacent sub-images of X_{k ι} ; and

15 means for combining each non-overlapping sub-image of X_{k ι} to obtain the entire image of X_{k >k} .

33. The encoder of claim 32, wherein τ= 2.

0 34. The encoder of claim 32, wherein τ= 3.

35. The encoder of claim 32 wherein the autocorrelation coefficient between zonal filter σ = 0.

5 36. The encoder of claim 32, wherein autocorrelation coefficients between zonal filter σ = μ = 0.

37. The encoder of claim 32, wherein autocorrelation coefficients between zonal filter η = λ = p = μ = σ = 0 .

30

38. The encoder of claim 32, wherein the WDFT is a 9 x 9 WDFT given by

8

X(k) = ∑ x(n) W^nk, k = 0, • • • , 8 , (41) and the sub-images of Y_{J J2} are 9 x 9.

77=0

35 39. The encoder of claim 32, wherein the WDFT is a 5 x 5 WDFT given by,

4 ,

X(k) = ∑χ(n) W^nk, * = o, • • ■ , 4 , (42) and the sub-images of Y_{Λ h} are 5 x 5.

1 40. The encoder of claim 32 further comprising a decoder having means for applying X_{k kι} in ∑ ∑ X_{k kι}A_{hh kι} = Y_{J h} , O ≤ j, ≤ n_t -1 and * = 1, 2. to

A-,=0 k₂=0 c perform a 2-D convolution of X_{k kl} and R(t ,t₂) to obtain a 2-D reconstructed image.

41. A method performed by a computer for coding a signal comprising the steps of: o defining a cubic-spline filter; correlating the filter with the signal to obtain a correlated signal; autocorrelating the filter to obtain autocorrelated filter coefficients; computing a transform of the correlated signal and the autocorrelated filter coefficients; 15 dividing the transform of the correlated signal by the transform of the autocorrelated filter coefficients to obtain a transform of a compressed signal; and computing an inverse transform of the transform of the compressed signal to obtain the compressed signal.

20

42. The method of claim 41, wherein the signal, the filter, and the transforms are one dimensional.

43. The method of claim 41, wherein the signal, the filter, and the transforms are 25 two dimensional.

44. The method of claim 41, further comprising the step of computing a convolution of the compressed signal and the filter to obtain a reconstructed signal.

30

45. The method of claim 41, wherein the step of computing a transform comprises computing a fast Fourier transform (FFT) and the step of computing an inverse transform comprises computing an inverse FFT.

35 46. The method of claim 41, wherein the step of computing a transform comprises computing a Winograd discrete Fourier transform (WDFT) with an overlap- save scheme and the step of computing an inverse transform comprises computing an inverse WDFT.

47. The method of claim 41, further comprising the step of defining a zonal filter to simplify the steps of correlating and autocorrelating.

48. A modified JPEG encoder for data compression comprising: a cubic-spline filter; means for correlating the filter with an input signal to obtain a correlated signal; means for autocorrelating the filter to obtain autocorrelated filter coefficients; means for computing a transform of the correlated signal and the autocorrelated filter coefficients; means for dividing the transform of the correlated signal by the transform of the autocorrelated filter coefficients to obtain a transform of a compressed signal; and means for computing an inverse transform of the transform of the compressed signal to obtain the compressed signal.

49. The encoder of claim 48, wherein the signal, the filter, and the transforms are one dimensional.

50. The encoder of claim 48, wherein the signal, the filter, and the transforms are two dimensional.

51. The encoder of claim 48, further comprising a decoder having means for computing a convolution of the compressed signal and the filter to obtain a reconstructed signal.

52. The encoder of claim 48, wherein the means for computing a transform comprises means for computing a fast Fourier transform (FFT) and the means for computing an inverse transform comprises means for computing an inverse FFT.

53. The encoder of claim 48, wherein the means for computing a transform comprises means for computing a Winograd discrete Fourier transform

(WDFT) with an overlap-save scheme and the means for computing an inverse transform comprises means for computing an inverse WDFT.

4. The encoder of claim 48, further comprising means for defining a zonal filter to simplify the means for correlating and the means for autocorrelating.