Solutions (9)

Assigned: March 20
Due: April 3

Fourier Spectrum and Power Spectrum

Compute and display the Fourier transform (using ListPlotComplex[…]) and the power spectrum (using ListPlot[Log[Abs[…]]] or ListPlot[Abs[…]]) for the following data sequences.

nX = 2^10 (* 1k Samples *)

x1 = Table[Mod[θ 0.2123,1], {θ, 0, nX -1}];

x2 = Table[If[Mod[θ 0.2123,2π] < π, 1, -1], {θ, 0, nX -1}];

x3 = Table[Sqrt[Abs[Sin[θ 0.2123]]], {θ, 0, nX -1}];

x1a = Table[Mod[θ 0.02123,1], {θ, 0, nX -1}];

x2a = Table[If[Mod[θ 0.02123,2π] < π, 1, -1], {θ, 0, nX -1}];

x3a = Table[Sqrt[Abs[Sin[θ 0.02123]]], {θ, 0, nX -1}];

Solution

ListPlotComplex[x_List,options___] := ListPlot[{Re[x], Im[x]},PlotStyle-> {Blue,Red}, options]

DoFourierPlots[x_List,xName_] := Module[{xf},
xf = Fourier[x];
GraphicsColumn[{ListPlot[x, PlotRange -> All, PlotLabel -> xName],
ListPlotComplex[xf, PlotRange -> All, PlotLabel -> xName],
ListPlot[Abs[xf], PlotRange -> All, PlotLabel -> xName, PlotStyle -> PointSize[0.01]]}, Frame -> All]
]

DoFourierPlots[#1, #2]& @@@ {{x1, "x1"}, {x2, "x2"}, {x3, "x3"}}

DoFourierPlots[#1, #2]& @@@ {{x1a, "x1a"}, {x2a, "x2a"}, {x3a, "x3a"}}

If I Had a Million Dollars (In-Class)

Using the same procedure used in class, import the sound track "Million_Dollar_64kb.aif" and create a frequency spectrogram.

Solution

Directory[]

SetDirectory["~/Desktop/"]

s = Import["Million_Dollar_64kb.aif"]

Show[s]

rate = Part[s,1,2](* Sample Rate: Samples Per Second *)

sData = Drop[Part[s,1,1,1],1]; (* Make exactly 2^16 = 64 kS *)

nData = Length[sData]

2^16

tfd[x_List,windowSize_:200, padSize_:1] := Module[{i,j,n,pad, window, xw, xf, result},
n = Length[x];
pad = Table[0, {i, 1, windowSize padSize}];
window = Table[If[i ≤ windowSize/2, (i-1) ,(windowSize - i)], {i, 1, windowSize}]/(0.5 windowSize);
result = {};
Do[xw = Take[x, {j,j+windowSize-1}]; 
xw = xw - (Plus @@ xw)/n; (* remove average *)
xw = Join[pad , (xw window), pad];
xf = Fourier[xw];
result = Append[result, Abs[Drop[xf,-Floor[Length[xf]/2]]]],
{j, 1, n-windowSize -1, windowSize/2}];
Transpose[result]]

sTFD = tfd[sData];

Dimensions[sTFD]

Print["[Max, Min] = ",Max[sTFD],", ",Min[sTFD]];ListDensityPlot[sTFD[[All,Range[1,654/2]]], PerformanceGoal -> "Quality",Mesh->False, DataRange -> {{0, nData/rate},{0,rate/2}}]

Band Pass Filter

Apply a band-pass filter to the signal above, and re-display the spectrogram.

Solution

sf = Fourier[sData];

bandpass[j_,lp_:5000,hp_:1000] := If[(hp < j < lp) || (nData - lp < j < nData - hp), 1.0, 0.0]

bpFilter = Table[bandpass[i],{i, 1, nData}];

sfBP = sf bpFilter;

(sfBP . Conjugate[sfBP] // Chop) / (sData .sData) (* energy fraction retained *)

xfBP = InverseFourier[sfBP];

ListPlay[Re[xfBP], SampleRate -> rate]

sBPTFD = tfd[Re[xfBP]];

ArrayPlot[Reverse[sTFD,1],Mesh->False,PlotRange->{0,0.02}]
ArrayPlot[Reverse[sBPTFD,1],Mesh->False,PlotRange->{0,0.02}]

Information Entropy

Introductory Comments

An important measure in communitcation and data theory 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, 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]))]

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 100 bits of data) is

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

Tally[i1]

Sort[Tally[i1]]

ListPlot[Sort[Tally[i1]], Filling-> Axis, AxesLabel -> {"Value", "Count"}]

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

First[Transpose[Reverse[Sort[Tally[i1]],2]]]/Length[i1]

Plus @@ %

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

Entropy[ ] is also a built-in function in Mathematica…

Entropy[2,i1] // N

For another example, a simple list with low entropy  (and requiring at least 30 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[…].

Assignment

What is the Entropy of the six signals in the first part of this assignment?

What is the Entropy of the sound track "Million_Dollar_64kb.aif" before and after you applied the bandpass filter?

For each of these problem, make sure that you quantize your data before applying the Entropy[…] function.

Solution

Part 1

Since "x1" and "x3" are real data lists, we first quantize these data into 8-bit integers.

x1 // Short

x2 // Short

x3 // Short

x1Quantized = Round[255 (x1 - Min[x1])/(Max[x1]-Min[x1])];
x1aQuantized = Round[255 (x1a - Min[x1a])/(Max[x1a]-Min[x1a])];

x1Quantized //Short

x3Quantized = Round[255 (x3 - Min[x3])/(Max[x3]-Min[x3])];
x3aQuantized = Round[255 (x3a - Min[x3a])/(Max[x3a]-Min[x3a])];

{#2, N[Entropy[2,#1]]}&@@@ {{x1Quantized, "x1"}, {x2, "x2"}, {x3Quantized, "x3"}}

{#2, N[Entropy[2,#1]]}&@@@ {{x1aQuantized, "x1a"}, {x2a, "x2a"}, {x3aQuantized, "x3a"}}

The entropy of "x2" is 7 times less than either "x1" or "x3". Is this expected?

Part 2

sData // Short

{Min[#], Max[#]}& [sData]

sDataQuantized = Round[255 (sData - Min[sData])/(Max[sData]-Min[sData])];

sDataQuantized // Short

N[Entropy[2,sDataQuantized]]

{Min[#], Max[#]}& [Re[xfBP]]

xfBPQuantized = Round[255 (Re[xfBP] - Min[Re[xfBP]])/(Max[Re[xfBP]]-Min[Re[xfBP]])];

N[Entropy[2,xfBPQuantized]]

The bandpass filter increased the information entropy! Part of the reason for this can be explained by the peak reduction created by filtering the high frequencies. If we quantize the signal to the original data scale, then we see that the information entropy is reduced.

xfBPQuantizedB = Round[255 (Re[xfBP] - Min[sData])/(Max[sData]-Min[sData])];

N[Entropy[2,xfBPQuantizedB]]

Information in FFT

Finally, if we compute the information entropy for the Fourier transformed data, then filtering of certain quasiharmonic data can significantly reduce the file size.  This is primarily because filtering reduces the bandwidth, or file size.

Compare data size…

{Min[#], Max[#]}& [Re[sf]]
{Min[#], Max[#]}& [Im[sf]]
{Min[#], Max[#]}& [Abs[sf]]
Print["Percent of data equal to zero = ",100.0 Count[sf, z_/; Abs[z] == 0]/Length[sf]]

{Min[#], Max[#]}& [Re[sfBP]]
{Min[#], Max[#]}& [Im[sfBP]]
{Min[#], Max[#]}& [Abs[sfBP]]
Print["Percent of data equal to zero = ",100.0 Count[sfBP, z_/; Abs[z] == 0]/Length[sfBP]]

Compare information entropy…

{Entropy[2,Re[sf]] ,Entropy[2,Im[sf]] }// N

{Entropy[2,Re[sfBP]] ,Entropy[2,Im[sfBP]] }// N