Discrete Cosine Transform and wavelets

Definitions

SetDirectory["/Users/knutm/Work/Calculus book/2011/sound"]

"/Users/knutm/Work/Calculus book/2011/sound"

cm=72/2.54

Disrete Cosine Transform

dct[u_]:=Module[{n=Length[u],c},
c=Sqrt[2/n]Table[Sum[u[[s+1]]Cos[Pi(2s+1)r/(2n)],{s,0,n-1}],{r,0,n-1}];
c[[1]]=c[[1]]/Sqrt[2];
c]

invdct[v_]:=Module[{n=Length[v],c},
c[r_]:=If[r==0,Sqrt[1/n],Sqrt[2/n]];Table[Sum[c[r]v[[r+1]]Cos[Pi(2s+1)r/(2n)],{r,0,n-1}],{s,0,n-1}]]

fam[n_,r_]:=Table[Cos[Pi(2s+1)r/(2n)],{s,0,n-1}];

inner[u_,v_]:=u.v

Wavelet Transforms

waveletdecompose[snd_,eps_]:=
Module[{c=snd,d,w},
If [EvenQ[length[c]],c=Drop[c,-1]];
d=Take[c,{1,-1,2}];
w=Take[ListConvolve[{-0.5,1,-0.5},c],{1,-1,2}];
w=Map[If[Abs[#]<eps,0,#]&,w];
{d,w}
];

waveletreconstruct[{d_,w_}]:=Riffle[d,ListConvolve[{0.5,0.5},d]+w];

Computations

invdct[dct[{0,1,2,3,4}]]

N[%]

Example with sound (making figures)

Import part of Here Comes the Sun

snd=Import["sun.wav"];

Dimensions[snd[[1,1,2]]]

snd0=Take[snd[[1,1,2]],{1000000,1022050}];

g0=ListPlot[snd0,PlotJoined->True,DataRange->{0,0.5},AxesOrigin->{0,-0.3}]

Export["sun0.pdf",g0,ImageSize->12cm]

snd00=Take[snd[[1,1,2]],{1000000,1001254}];

g00=ListPlot[snd00,PlotJoined->True,DataRange->{0,0.5/32},AxesOrigin->{0,-0.3}]

Export["sun1.pdf",g00,ImageSize->8cm]

snd000=Take[snd[[1,1,2]],{1000000,1000156}];

g000=ListPlot[snd000,PlotJoined->True,DataRange->{0,0.5/256},AxesOrigin->{0,-0.3}]

Export["sun2.pdf",g000,ImageSize->8cm]

snd1=Take[snd[[1,1,2]],{1000000,1000400}];

Using the Discrete Cosine Transform

g=ListPlot[snd1,PlotJoined->True,ImageSize->6cm]

Export["data.pdf",g]

fouriersnd=FourierDCT[snd1];

g=ListPlot[fouriersnd,PlotRange->All,PlotJoined->True,ImageSize->6cm]

Export["dct.pdf",g]

Compression

fouriersndeps=Map[If[Abs[#]<0.02,0,#]&,fouriersnd];

g=ListPlot[fouriersndeps,PlotRange->All,PlotJoined->True,ImageSize->6cm]

Export["dcteps.pdf",g]

Count[fouriersndeps,0]

Inverse transform

snd2=FourierDCT[fouriersndeps,3];

g=ListPlot[snd2,PlotJoined->True,ImageSize->6cm]

Export["datadcteps.pdf",g]

g=ListPlot[snd1-snd2,PlotJoined->True,ImageSize->6cm]

Export["dcterr.pdf",g]

Wavelet transform

Decompose

ListPlot[snd1,PlotJoined->True]

ListPlay[snd1,SampleRate->44100]

{d,w}=waveletdecompose[snd1,0];

g=ListPlot[d,PlotJoined->True,ImageSize->6cm]

Export["wavd.pdf",g]

g=ListPlot[w,PlotJoined->True,ImageSize->6cm]

Export["wavw.pdf",g]

{d,weps}=waveletdecompose[snd1,0.02];

Count[weps,0]

Length[weps]

g=ListPlot[weps,PlotJoined->True,ImageSize->6cm]

Export["wavweps.pdf",g]

cc=waveletreconstruct[{d,weps}];

g=ListPlot[cc,PlotJoined->True,ImageSize->6cm]

Export["wavcceps.pdf",g]

g=ListPlot[cc-snd1,PlotJoined->True,PlotRange->All,ImageSize->6cm]

Export["waverr.pdf",g]

Example with sound

Import part of Here Comes the Sun

snd=Import["/Users/knutm/Work/Calculus book/2011/sound/sun.wav"];

snd1=Take[snd[[1,1,2]],{1000000,1300000}];

ListPlay[snd1,SampleRate->44100]

Discrete Cosine Transform

fouriersnd=FourierDCT[snd1];

ListPlot[fouriersnd,PlotRange->All,PlotJoined->True]

ListPlay[fouriersnd,SampleRate->44100]

file=$TemporaryPrefix<>"dct";
BinaryWrite[file,fouriersnd,"Real32"];
Close[file]

Set small coefficients to zero

fouriersndeps=Map[If[Abs[#]<0.01,0,#]&,fouriersnd];

ListPlot[fouriersndeps,PlotJoined->True,PlotRange->All,PlotRange->All]

Count[fouriersndeps,0]

file=$TemporaryPrefix<>"dcteps";
BinaryWrite[file,fouriersndeps,"Real32"];
Close[file]

Inverse transform

snd2=FourierDCT[fouriersndeps,3];

ListPlay[snd2,SampleRate->44100]

Max[Abs[snd2-snd1]]

Wavelet transform

Write to file

file=$TemporaryPrefix<>"lyd";
BinaryWrite[file,snd1,"Real32"];
Close[file]

Decompose

{d,w}=waveletdecompose[snd1,0];

Write to file

file=$TemporaryPrefix<>"waveletlyd";
BinaryWrite[file,d,"Real32"];
BinaryWrite[file,w,"Real32"];
Close[file]

Play the sound

{ListPlay[d,SampleRate->22050],ListPlay[w,SampleRate->22050]}

Allow error in decomposition

{d,weps}=waveletdecompose[snd1,0.02];

Number of points that are zero

Count[weps,0]

Length[weps]

ListPlay[weps,SampleRate->22050]

Write compressed sound

file=$TemporaryPrefix<>"waveletepslyd";
BinaryWrite[file,d,"Real32"];
BinaryWrite[file,weps,"Real32"];
Close[file]

Reconstruct

cc=waveletreconstruct[{d,weps}];

Compute error

Max[Abs[cc-snd1]]

Play compressed sound

ListPlay[cc,SampleRate->44100]

Some figures

gs=Plot[101325.+If[t<0.1,0,Min[1.,3.78(t-0.1)^2]Sin[32.Pi (t-0.1)]],{t,0,1}]

Export["sounda.pdf",gs,ImageSize->8cm];

gs1=Plot[0.25Sin[40Pi t]+Cos[4Pi t],{t,0,2}]

Export["soundb.pdf",gs1,ImageSize->8cm]

g0=Plot[Sin[t],{t,0,2Pi}]

Export["sin0.pdf",g0,ImageSize->8cm]

g1=Plot[Sin[2Pi t],{t,0,1}]

Export["sin1.pdf",g1,ImageSize->8cm]

g2=Plot[Sin[2Pi 2 t],{t,0,1}]

Export["sin2.pdf",g2,ImageSize->8cm]

g3=Plot[Sin[2Pi 5 t],{t,0,1}]

Export["sin3.pdf",g3,ImageSize->8cm]