Compression and Reconstruction Accuracy of Digital Data

Student Projects in AP1603 investigating Fourier transforms.

All projects must have…
Name:____________
Email:_____@columbia.edu

Introduction

The second general theme of our Mathematica investigations is the Fourier transform. In our first lecture, we examined how a complete expansion of a function in terms of simple Cosine and Sine functions allowed the solution of a partial differential equation. The solution was represented as an infinite sum of functions. However, we learned that if we take only the most significant (or largest amplitude) functions, the difference between the actual signal and the approximation can be relatively small. 

In our next two lectures, we examined the application of Fourier transforms to digital data. We explored the active field of digital signal processing. Since the invention of the computer and the FFT, large amounts of data can be transformed rapidly. Hardware implementations of these algorithms have lead to the many commercial products found in commercial electronics (like the DVD, streaming media, and accurate compression of large photographs.) 

For your second research project, I suggest that you determine the relationship between the degree of data compression and the accuracy of the reconstructed data. You may experiment with either audio or graphical data or both. You may look at the files posted on our web site, or import or create your own data sets. Try to limit your analysis to one type of data file so that your project explores a topic in depth.

In order to carry out a quantitative analysis, you must define measurements of data compression and of data reconstruction accuracy. Below, I will define three measures that you should feel free to use. The first is the information entropy. It is a measure of the disorder of the data. Multiplying the entropy by the number of data points is the minimum possible file size required to contain the information within the data. The second is the root mean square relative error (RMSE) of the reconstruction relative to the "real" data. It can be related to the signal to noise ratio (SNR). When the RMSE is large (or the SNR is small), the reconstruction becomes inaccurate. These two measurements can be used to quantify the trade-off between file size and accuracy. The third is the absolute error. The absolute error is, perhaps, the easiest to understand; it represents the average difference between the original and reconstructed image or sound over the entire data set.

Off[General::spell1];

Information Entropy

The first measure that you can use in your project is called "information entropy". Entropy is a concept first developed for statistical physics. It is a measure of the disorder of a system. When dealing with digital data, the entropy measures the disorder (or complexity) of the data. When the entropy is small, the data is not disordered, and only a few digital values exist within the entire data set. For example, if you had 1000 data points with values all equal to 1, then the information entropy would be exactly zero. The data is not at all disordered, and (theoretically) you could represent the information contained within this data set by one bit. When the entropy is large, the data is very disordered, and many digital values exist in the data set.

The Mathematica function for information entropy is written as

iEntropy[list_List]:=-Plus@@(N[#1 Log[2,#1]]&)/@((First[Transpose[Reverse[Sort[Tally[list]],2]]])/(Length[list]))

Characters["mississippi"]

(((First[Transpose[Reverse[Sort[Tally[#1]],2]]])/(Length[#1]))&)@@{Characters["mississippi"]}

iEntropy[Characters["mississippi"]]

This function computes the sum of the logs of the "probability" of each value in the data.

For example, a simple list with high entropy (and requiring at least 1000 bits of data) is

i1=RandomInteger[{1,10},30]

Reverse[Sort[Tally[i1]],2]

Print["Entropy = ",iEntropy[i1], "; Bits Required = ", iEntropy[i1] Length[i1]];

For example, a simple list with low entropy  (and requiring at least 60 bits of data) is

i2 = Table[If[i<15,1, 2], {i, 1, 30}]

Reverse[Sort[Tally[i2]],2]

Print["Entropy = ",iEntropy[i2], "; Bits Required = ", iEntropy[i2] Length[i2]];

As the data becomes more orderly, the information entropy decreases.

For two-dimensional data, your must first Flatten[…] the data list before using iEntropy[…]. This converts the multiple-dimensional data to a one-dimensional list having the form acceptable to Frequencies[…].

You should be very aware of the difference between Real and Integer numbers. When you import an audio file (.aif), the integer (8-bit) data is converted to Real numbers. The data is still "quantized" into 256 different real numbers. The iEntropy[…] function works with quantized real numbers as well as integers. However, once you transform this data (using either a FFT or a DCT), the transform will not be quantized. The information entropy will grow considerably. Therefore, in examining your data whether you use JPEG-weighted quantization or a uniform quantization, you must first quantize your data before you measure the information content with iEntropy[…].

Root Mean Square Relative Error

The root mean square relative error (RMSE) is a good measure of the difference between two data sequences. The function to compute the RMSE is defined by

rmsError[l1_List,l2_List] := Sqrt[Plus @@ (#^2& /@ (l1 - l2))] / Sqrt[Plus @@ (l1 l1)]

rmsError[{a,b,c},{a2,b2,c2}]

(Remember, the "energy" within a signal is defined by "Plus @@ (data data)".)

noise=0.1`;
i1=N[Table[Sin[((i 3 π)/100)] ((Sin[((i π)/100)])^2),{i,1,100}]];
i2=i1+RandomReal[{-noise,noise},100];
Show[ListPlot[i1,DisplayFunction->Identity],ListPlot[i2,DisplayFunction->Identity,PlotStyle->RGBColor[1,0,0]],DisplayFunction->$DisplayFunction]

rmsError[i1,i2]

noise=0.2`;
i1=N[Table[Sin[((i 3 π)/100)] ((Sin[((i π)/100)])^2),{i,1,100}]];
i2=i1+RandomReal[{-noise,noise},100];
Show[ListPlot[i1,DisplayFunction->Identity],ListPlot[i2,DisplayFunction->Identity,PlotStyle->RGBColor[1,0,0]],DisplayFunction->$DisplayFunction]

rmsError[i1,i2]

The RMSE is related to the signal to noise ratio (SNR). This formula is another measure of relative error expressed in "decibels".

snr[l1_List, l2_List] := - 20 Log[rmsError[l1,l2]]

snr[i1,i2]

As the signal to noise ratio increases, then there is less distortion or degradation of the digital signal.

Average Absolute Error

The average absolute error is the average difference between the original and reconstructed data.

absError[l1_List,l2_List] := Plus @@ Abs /@ (l1 - l2) / Length[l1]

absError[i1,i2]

Question for Study

One of the most important questions in digital processing is how to compact data so that essential information is not loss. In practice, the trade off between data size and information accuracy is subjective. Choices are made by a "trial and test" approach. Usually, users can select the degree of compaction corresponding to application requirements.

There are many methods for data compression, but we have discussed only two methods in class. These are (1) filtering and (2) quantization. (Another rather basic method to reduce file-size is "downsampling". Downsampling corresponds to a reduction in data rate or image resolution.)

I suggest that you investigate filtering and quantization on the information entropy of data files. As you change filter parameters or the degree of quantization, the entropy (and therefore the data size) can decrease significantly. What is the resulting increase in the RMSE or the decrease in the signal to noise as the file becomes more compact.

In the following, filter and quantization methods are summarized below. We use only the Discrete Cosine Transform (DCT) since both filtering and quantization can be performed more simply than with the FFT.

Definitions

In order to help with your projects, the following summarizes useful definitions for the discrete Cosine transform (DCT).

DCT of 1D (audio) Data

dct1D[x_List] := Module[{n,e0, tFactors, xeven, xodd},
n = Length[x];
{xeven, xodd} = Transpose[Partition[x, 2]];
e0 = Exp[-I π/(2 n)];
tFactors = {1} ~ Join~ Table[Sqrt[2] e0^i, {i, 1, n-1}];
Re[tFactors InverseFourier[Join[xeven, Reverse[xodd]]]]]

idct1D[x_List] :=  Module[{n,e0,tFactors, xt, xfirst, xlast},
n = Length[x];
e0 = Exp[-I π/(2 n)];
tFactors = {1} ~ Join~ Table[Sqrt[2] e0^i, {i, 1, n-1}];
tFactors = Conjugate[N[tFactors]];
xt = Re[Fourier[tFactors x]];
xfirst = Take[xt, n/2];
xlast = Take[xt,-n/2];
Flatten[Transpose[{xfirst, Reverse[xlast]}]]]

Blocked DCT of 2D (image) Data

bDCTfactors = N[{1} ~ Join~ Table[Sqrt[2] Exp[-I π i/(2 8)], {i, 1, 8-1}]];
bIDCTfactors = Conjugate[bDCTfactors];
bDCT[list_?(Length[#] == 8 &)] := Re[ bDCTfactors InverseFourier[ N[list[[{1,3,5,7,8,6,4,2}]]]]];
bIDCT[list_?(Length[#] == 8 &)] := Re[Fourier[bIDCTfactors list]][[{1,8,2,7,3,6,4,5}]];
bDCT[array_?MatrixQ] := Transpose[bDCT /@ Transpose[bDCT /@ array]];
bIDCT[array_?MatrixQ] := Transpose[bIDCT /@ Transpose[bIDCT /@ array]];
bDCT[list_?(Length[#] > 8&)] := Join @@ (bDCT /@ Partition[list,8]);
bIDCT[list_?(Length[#] > 8&)] := Join @@ (bIDCT /@ Partition[list,8]);

Quantization

You can quantize your data with the Round[…] command. Uniform quantization simply rounds all data values to equally-spaced values ranging from the maximum to the minimum of your data. JPEG quantization is an example of a non-uniform quantization. Some of the 2D DCT transform coefficients are weighted more strongly during the quantization process.

A useful command to perform uniform quantization of Real data is…

quantizeUniform[x_List,bits_:8] := Module[{q},
q = 2^bits/(Max[x] - Min[x]);
Round[q x]/q];

The quantization matrix for JPEG quantization is given by

qLum =
{{16, 11, 10, 16, 24, 40, 51, 61},
{12, 12, 14, 19, 26, 58, 60, 55},
{14, 13, 16, 24, 40, 57, 69, 56},
{14, 17, 22, 29, 51, 87, 80, 62},
{18, 22, 37, 56, 68,109,103, 77},
{24, 35, 55, 64, 81,104,113, 92},
{49, 64, 78, 87,103,121,120,101},
{72, 92, 95, 98,112,100,103, 99}};
qLum // MatrixForm

This is applied to the DCT of 8 x 8 blocks of an image.

Functions used to perform block quantization of graphical images…

BlockImage[image_, blocksize_:{8, 8}] :=Partition[image, blocksize] /;And @@ IntegerQ /@ (Dimensions[image]/blocksize)

UnBlockImage[blocks_] := Partition[Flatten[Transpose[blocks, {1, 3, 2}]],{Times @@ Dimensions[blocks][[{2, 4}]]}]

bQuantize[image_, qMatrix_] := Map[(#/qMatrix)&, BlockImage[image, Dimensions[qMatrix]],{2}] //UnBlockImage // Round;

bUnQuantize[image_, qMatrix_] := Map[(# qMatrix)&, BlockImage[image, Dimensions[qMatrix]],{2}] //UnBlockImage

Filtering (using the DCT)

The approximate "frequencies" associated with the DCT are (m × rate)/(2 × N); where rate  is the sampling rate (samples per second), N is the total number os samples, and m ranges from 0 to N – 1.

A generalized bandpass filter can be written as

bpFilter[xDCT_List,fLow_,fHigh_,rate_:41000.0] := Module[{n, filter},
n = Length[xDCT];
filter = Table[If[fLow < i rate/(2 n)  < fHigh, 1, 0], {i, 1, n}];
xDCT filter];

To use this filter, first take the DCT of a one-dimensional signal, and then apply the filter. The iDCT reconstructs the filtered signal. The function bpFilter[…] takes an optional argument pre-set to the sample rate found in audio CDs.

Note that the information entropy of the filtered DCT may be smaller than the reconstructed signal. For example, Dylan Thomas's voice has the following properties.

Example filter of Dylan Thomas's voice

SetDirectory["~/Desktop/"]; (* Set path to data file *)
s = Import["Do_Not_Go_Gentle_64kb.aif"];
rate = Part[s, 1, 2];
sData = Drop[Part[s,1,1,1],1]; (* Make exactly 2^16 = 64 kS *)

sDCT = dct1D[sData];

For unfiltered data, the DCT compacts the file by about a factor of three with 8 bit quantization…

Print["Data Entropy = ",iEntropy[sData],"; DCT (REAL) Entropy = ",iEntropy[sDCT], "; DCT (quantized) Entropy = ", iEntropy[quantizeUniform[sDCT]]];

sfDCT = quantizeUniform[bpFilter[sDCT,200.0, 4000.0]]; (* filter from 200 Hz to 5 kHz and 8-bit quantization *)

Print["Filtered DCT Entropy = ", iEntropy[sfDCT]];

sReconstructed = idct1D[sfDCT];

rmsError[sData,sReconstructed]

snr[sData,sReconstructed]

absError[sData,sReconstructed]

Surprisingly, these measurements show a great variability. The RMSE (and the SNR) indicate very, very poor reconstruction accuracy; whereas, the average absolute error indicates a relatively good reconstruction. Judge for yourself. Try ListPlay[…] and listen to the differences.

ListPlay[sData,SampleRate -> rate]

ListPlay[sReconstructed,SampleRate -> rate]

Image Compression using Blocked DCT

In class, we performed  8 x 8 block DCT of two-dimensional images. An example from class is presented below. Our class algorithm works only for images that can be evenly divided into 8 x 8 blocks.

Dylan Thomas's Portrait

SetDirectory["~/Desktop/"];
g = Import["D_Thomas.tif"];
gData = ImageData[g,"Byte"];

ArrayPlot[gData,Mesh->False,AspectRatio->Automatic, ColorFunction -> GrayLevel]

Perform the block DCT…

bgData = bDCT[gData];

Quantize…

bqgData = bQuantize[bgData,qLum];

Print["Data Entropy = ",iEntropy[Flatten[gData]],"; DCT (REAL) Entropy = ",iEntropy[Flatten[bgData]], "; DCT (quantized) Entropy = ", iEntropy[Flatten[bqgData]]];

Finally, we perform an inverse transform in order to check reconstruction accuracy…

gReconstructed=bIDCT[bUnQuantize[bqgData,qLum]];
ArrayPlot[gReconstructed,Mesh->False,AspectRatio->Automatic, ColorFunction -> GrayLevel]

rmsError[Flatten[gData],Flatten[gReconstructed]]

absError[Flatten[gData],Flatten[gReconstructed]]

For image data, the RMSE (and SNR) appear to be better measurements for the study of blocked data.