(* Content-type: application/mathematica *) (*** Wolfram Notebook File ***) (* http://www.wolfram.com/nb *) (* CreatedBy='Mathematica 6.0' *) (*CacheID: 234*) (* Internal cache information: NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 145, 7] NotebookDataLength[ 37432, 1172] NotebookOptionsPosition[ 31893, 990] NotebookOutlinePosition[ 33856, 1056] CellTagsIndexPosition[ 33785, 1051] WindowFrame->Normal ContainsDynamic->False*) (* Beginning of Notebook Content *) Notebook[{ Cell[CellGroupData[{ Cell["Denoising with Wavelet Transforms", "Title", CellChangeTimes->{{3.4084937095*^9, 3.40849371428125*^9}}], Cell["\<\ Wavelet Workshop June 4-7, 2008 University of St. Thomas\ \>", "Subtitle", CellChangeTimes->{{3.421426484890625*^9, 3.42142648821875*^9}}], Cell["\<\ Catherine Beneteau Caroline Haddad David Ruch Patrick Van Fleet\ \>", "Subsubtitle", CellChangeTimes->{{3.42144927371875*^9, 3.42144928703125*^9}}], Cell[CellGroupData[{ Cell["Objectives", "Section"], Cell["\<\ The purpose of this notebook is to introduce you to some basic wavelet-based \ denoising methods.\ \>", "Text", CellChangeTimes->{3.421755817067687*^9}] }, Open ]], Cell[CellGroupData[{ Cell["Conventions", "Section"], Cell[TextData[{ "This notebook uses the package ", StyleBox["DiscreteWavelets", FontColor->RGBColor[1, 0, 0]], " (written by Patrick Van Fleet). All commands from the ", StyleBox["DiscreteWavelets", FontColor->RGBColor[1, 0, 0]], " library will be denoted in ", StyleBox["red", FontColor->RGBColor[1, 0, 0]], ". Help is available for every command in ", StyleBox["the package", FontColor->GrayLevel[0]], ". Click on Help and then Documentation Center. At the bottom-right of the \ page is a link for Installed AddOns. Click this link and one of the options \ is DiscreteWavelets. Click this link to go to the Help Browser. Like all ", StyleBox["Mathematica", FontSlant->"Italic"], " help screens, the help is \"live\" - you can either execute the commands \ in the help to see the effects of the command or cut and paste them into your \ own notebook.\n\nComments are useful within cells of code. Any code enclosed \ by (* *) is a comment and ignored by the ", StyleBox["Mathematica", FontSlant->"Italic"], " kernel." }], "Text", CellChangeTimes->{{3.4085631561875*^9, 3.408563201375*^9}, { 3.4085632433125*^9, 3.408563315796875*^9}, {3.40856339940625*^9, 3.4085634034375*^9}, {3.408565167890625*^9, 3.40856516825*^9}}] }, Open ]], Cell[CellGroupData[{ Cell["Load DiscreteWavelets and Audio", "Section", CellChangeTimes->{ 3.408563584328125*^9, {3.4217558514279413`*^9, 3.4217558531154847`*^9}}], Cell[BoxData[{ RowBox[{"<<", "DiscreteWavelets`DiscreteWavelets`"}], "\[IndentingNewLine]", RowBox[{"<<", "Audio`"}]}], "Input", CellChangeTimes->{{3.408563455078125*^9, 3.408563467234375*^9}, { 3.408565207546875*^9, 3.40856520834375*^9}, {3.40856604671875*^9, 3.408566047*^9}, {3.4217558554436693`*^9, 3.4217558589593844`*^9}, { 3.4217575361585693`*^9, 3.421757536299198*^9}}] }, Open ]], Cell[CellGroupData[{ Cell["The Shrinkage Function", "Section"], Cell[TextData[{ "The shrinkage function is defined in ", StyleBox["DiscreteWavelets ", FontColor->RGBColor[1, 0, 0]], "and called ", StyleBox["ShrinkageFunction.", FontFamily->"Arial", FontColor->RGBColor[1, 0, 0]], " Here is the shrinkage function plotted for \[Lambda]=2." }], "Text"], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{ StyleBox["ShrinkageFunction", FontColor->RGBColor[1, 0, 0]], "[", RowBox[{"v", ",", "2"}], "]"}], ",", RowBox[{"{", RowBox[{"v", ",", RowBox[{"-", "5"}], ",", "5"}], "}"}], ",", RowBox[{"Ticks", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"-", "4"}], ",", RowBox[{"-", "3"}], ",", RowBox[{"-", "2"}], ",", RowBox[{"-", "1"}], ",", "1", ",", "2", ",", "3", ",", "4"}], "}"}], ",", RowBox[{"{", RowBox[{ RowBox[{"-", "3"}], ",", RowBox[{"-", "2"}], ",", RowBox[{"-", "1"}], ",", "1", ",", "2", ",", "3"}], "}"}]}], "}"}]}], ",", RowBox[{"PlotStyle", "\[Rule]", "Blue"}]}], "]"}]], "Input", CellChangeTimes->{{3.4086583225*^9, 3.408658322703125*^9}}], Cell["Let's try it on an example", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"v", "=", RowBox[{"{", RowBox[{ RowBox[{"-", "3.2"}], ",", RowBox[{"-", "2.5"}], ",", "1", ",", "1", ",", "1.1", ",", "5", ",", "3.2", ",", "1.9", ",", RowBox[{"-", "4"}], ",", RowBox[{"-", "6"}]}], "}"}]}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"lambda", "=", "2"}], ";"}], "\[IndentingNewLine]", RowBox[{"Map", "[", RowBox[{ RowBox[{ RowBox[{ StyleBox["ShrinkageFunction", FontColor->RGBColor[1, 0, 0]], "[", RowBox[{"#", ",", "lambda"}], "]"}], "&"}], ",", "v"}], "]"}]}], "Input",\ CellChangeTimes->{{3.408658505125*^9, 3.408658505203125*^9}}] }, Open ]], Cell[CellGroupData[{ Cell["A Sample Signal ", "Section"], Cell["\<\ For our example, we will use a sample signal suggested by Donoho. He calls \ this function the Heavisine function. \ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{"Heavisine", "[", "t_", "]"}], ":=", RowBox[{ RowBox[{"4", "*", RowBox[{"Sin", "[", RowBox[{"4", "*", "Pi", "*", "t"}], "]"}]}], "-", RowBox[{"Sign", "[", RowBox[{"t", "-", ".3"}], "]"}], "-", RowBox[{"Sign", "[", RowBox[{".72", "-", "t"}], "]"}]}]}], ";"}]], "Input", CellChangeTimes->{{3.40865799753125*^9, 3.408657998125*^9}, { 3.408658503453125*^9, 3.408658503546875*^9}}], Cell[TextData[{ "Let's sample this function to create our true signal ", StyleBox["v", FontWeight->"Bold"], " and then ", StyleBox["ListPlot", FontWeight->"Bold"], " it." }], "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"n", "=", "2048"}], ";"}], "\n", RowBox[{ RowBox[{"v", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"Heavisine", "[", FractionBox["k", "n"], "]"}], ",", RowBox[{"{", RowBox[{"k", ",", "0", ",", RowBox[{"n", "-", "1"}]}], "}"}]}], "]"}]}], ";"}], "\n", RowBox[{"TruSig", "=", RowBox[{"ListPlot", "[", RowBox[{"v", ",", RowBox[{"PlotStyle", "\[Rule]", "Red"}]}], "]"}]}]}], "Input", CellChangeTimes->{{3.40847228309375*^9, 3.4084722835*^9}, { 3.408473703734375*^9, 3.408473710734375*^9}, {3.408658004390625*^9, 3.408658005859375*^9}, {3.40865850084375*^9, 3.4086585009375*^9}}] }, Open ]], Cell[CellGroupData[{ Cell["Creating Noise", "Section"], Cell["Now we create some Gaussian white noise.", "Text", CellChangeTimes->{{3.4084732154375*^9, 3.408473216375*^9}}], Cell[BoxData[{ RowBox[{ RowBox[{"nd", "=", RowBox[{"NormalDistribution", "[", RowBox[{"0.", ",", "1."}], "]"}]}], ";"}], "\n", RowBox[{ RowBox[{"SeedRandom", "[", "]"}], ";"}], "\n", RowBox[{ RowBox[{"noise", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"Random", "[", "nd", "]"}], ",", RowBox[{"{", RowBox[{"k", ",", "1", ",", "n"}], "}"}]}], "]"}]}], ";"}], "\n", RowBox[{"ListPlot", "[", "noise", "]"}]}], "Input", CellChangeTimes->{{3.408658026953125*^9, 3.4086580284375*^9}, { 3.40865849909375*^9, 3.408658499203125*^9}}], Cell[TextData[{ "If you want, you can actually \"play\" the noise. ", StyleBox["Please be considerate of your neighbors ", FontSlant->"Italic"], ":-)", StyleBox[" ", FontSlant->"Italic"] }], "Text", CellChangeTimes->{{3.408473125390625*^9, 3.40847316109375*^9}}], Cell[BoxData[ RowBox[{"ListPlay", "[", RowBox[{"noise", ",", RowBox[{"SampleRate", "\[Rule]", "11025"}]}], "]"}]], "Input", CellChangeTimes->{ 3.408472640296875*^9, {3.40865849715625*^9, 3.4086584973125*^9}}], Cell["\<\ We create the noisy vector by picking a noise level and simply adding v to a \ scaled version of the noise.\ \>", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"sigma", "=", "0.5"}], ";"}], "\n", RowBox[{ RowBox[{"y", "=", RowBox[{"v", "+", RowBox[{"sigma", " ", "noise"}]}]}], ";"}], "\n", RowBox[{"plt", "=", RowBox[{"ListPlot", "[", RowBox[{"y", ",", RowBox[{"Ticks", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{ "0", ",", "0.2", ",", "0.4", ",", "0.6", ",", "0.8", ",", "1"}], "}"}], ",", RowBox[{"{", RowBox[{ RowBox[{"-", "6"}], ",", RowBox[{"-", "4"}], ",", RowBox[{"-", "2"}], ",", "0", ",", "2", ",", "4", ",", "6"}], "}"}]}], "}"}]}], ",", RowBox[{"PlotStyle", "\[Rule]", "Brown"}]}], "]"}]}]}], "Input", CellChangeTimes->{{3.40865804496875*^9, 3.4086580564375*^9}, { 3.408658495171875*^9, 3.408658495265625*^9}}] }, Open ]], Cell[CellGroupData[{ Cell["Estimating the Noise Level", "Section"], Cell[TextData[{ "To estimate the noise level, we use the ", StyleBox["MAD", FontColor->RGBColor[1, 0, 0]], " function in the ", StyleBox["DiscreteWavelets", FontColor->RGBColor[1, 0, 0]], " package. We have to first compute the wavelet transform. We'll use the \ Daubechies 4 filter for this example and compute 2 iterations." }], "Text", CellChangeTimes->{{3.4086580816875*^9, 3.408658094296875*^9}}], Cell[BoxData[{ RowBox[{ RowBox[{"h", "=", RowBox[{"N", "[", RowBox[{ StyleBox["Daub", FontColor->RGBColor[1, 0, 0]], "[", "4", "]"}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"Print", "[", RowBox[{ "\"\\"", ",", "h", ",", "\"\<.\>\""}], "]"}], ";"}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"i", "=", "5"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"wt", "=", RowBox[{ StyleBox["WT1D", FontColor->RGBColor[1, 0, 0]], "[", RowBox[{"y", ",", "h", ",", RowBox[{"NumIterations", "\[Rule]", "i"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ StyleBox["WaveletVectorPlot", FontColor->RGBColor[1, 0, 0]], "[", RowBox[{"wt", ",", RowBox[{"NumIterations", "\[Rule]", "i"}]}], "]"}]}], "Input", CellChangeTimes->{{3.408658126984375*^9, 3.408658129203125*^9}, { 3.40865849325*^9, 3.40865849340625*^9}}], Cell[TextData[{ "We need to extract the first highpass portion. There is a command in ", StyleBox["DiscreteWavelets", FontFamily->"Arial", FontColor->RGBColor[1, 0, 0]], " to do this: ", StyleBox["WaveletVectorToList", FontColor->RGBColor[1, 0, 0]], " takes a wavelet vector built using i iterations and creates a list of \ length i+1. The first element in the output is the lowpass portion, the \ second is the ith highpass portion,..., the last is the first highpass \ portion." }], "Text", CellChangeTimes->{{3.4217559984473305`*^9, 3.4217560609489307`*^9}}], Cell[BoxData[{ RowBox[{ RowBox[{"wtlist", "=", RowBox[{"WaveletVectorToList", "[", RowBox[{"wt", ",", RowBox[{"NumIterations", "\[Rule]", "i"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"lowpass", "=", RowBox[{"First", "[", "wtlist", "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"hp1", "=", RowBox[{"Last", "[", "wtlist", "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{"ListPlot", "[", "hp1", "]"}]}], "Input", CellChangeTimes->{{3.40847269415625*^9, 3.40847269503125*^9}, { 3.408658492140625*^9, 3.40865849221875*^9}, {3.4217559299455767`*^9, 3.421755988197068*^9}}], Cell["Now we estimate the noise:", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"sigmahat", "=", RowBox[{ RowBox[{ StyleBox["MAD", FontColor->GrayLevel[0.500008]], "[", "hp1", "]"}], "/", ".6745"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Print", "[", RowBox[{ "\"\\"", ",", "sigmahat", ",", "\"\<.\>\""}], "]"}], ";"}]}], "Input", CellChangeTimes->{{3.40865849109375*^9, 3.408658491203125*^9}, 3.4217560822776012`*^9}] }, Open ]], Cell[CellGroupData[{ Cell["VISUShrink", "Section"], Cell["\<\ To perform VISUShrink, we compute the universal tolerance and use it in the \ shrink functions on the highpass portions. \ \>", "Text", CellChangeTimes->{{3.4217560884652596`*^9, 3.421756091090327*^9}}], Cell[BoxData[{ RowBox[{ RowBox[{"wtlist", "=", RowBox[{ StyleBox["WaveletVectorToList", FontColor->RGBColor[1, 0, 0]], "[", RowBox[{"wt", ",", RowBox[{"NumIterations", "\[Rule]", "i"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"lowpass", "=", RowBox[{"First", "[", "wtlist", "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"highpass", "=", RowBox[{"Flatten", "[", RowBox[{"Drop", "[", RowBox[{"wtlist", ",", "1"}], "]"}], "]"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ "Here", " ", "is", " ", "the", " ", "VisuShrink", " ", "tolerance"}], " ", "*)"}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{"lambda", "=", RowBox[{"sigmahat", "*", RowBox[{"Sqrt", "[", RowBox[{"2", "*", RowBox[{"Log", "[", RowBox[{"Length", "[", "highpass", "]"}], "]"}]}], "]"}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"Print", "[", RowBox[{ "\"\\"", ",", "lambda", ",", "\"\<.\>\""}], "]"}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ "Build", " ", "the", " ", "new", " ", "highpass", " ", "portions"}], " ", "*)"}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"newhighpass", "=", RowBox[{"Map", "[", RowBox[{ RowBox[{ RowBox[{ StyleBox["ShrinkageFunction", FontColor->RGBColor[1, 0, 0]], "[", RowBox[{"#", ",", "lambda"}], "]"}], "&"}], ",", "highpass"}], "]"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ "Join", " ", "the", " ", "new", " ", "highpass", " ", "with", " ", "the", " ", "original", " ", "lowpass"}], " ", "*)"}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{"newwt", "=", RowBox[{"Join", "[", RowBox[{"lowpass", ",", "newhighpass"}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ StyleBox["WaveletVectorPlot", FontColor->RGBColor[1, 0, 0]], "[", RowBox[{"newwt", ",", RowBox[{"NumIterations", "\[Rule]", "i"}], ",", RowBox[{"PlotRange", "\[Rule]", "All"}]}], "]"}]}], "Input", CellChangeTimes->{ 3.408472720359375*^9, {3.4086522890625*^9, 3.408652290453125*^9}, { 3.408658489296875*^9, 3.40865848940625*^9}, {3.421756105934457*^9, 3.421756134263307*^9}}], Cell["\<\ We apply the inverse wavelet transform to obtain our estimate of v:\ \>", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"vhat", "=", RowBox[{ StyleBox["IWT1D", FontColor->RGBColor[1, 0, 0]], "[", RowBox[{"newwt", ",", "h", ",", RowBox[{"NumIterations", "\[Rule]", "i"}]}], "]"}]}], ";"}], "\n", RowBox[{"SignalEst", "=", RowBox[{"ListPlot", "[", RowBox[{"vhat", ",", RowBox[{"PlotStyle", "\[Rule]", "Brown"}]}], "]"}]}]}], "Input", CellChangeTimes->{{3.4084736729375*^9, 3.4084736773125*^9}, { 3.40865842678125*^9, 3.40865842690625*^9}}], Cell["Here is a visual comparison:", "Text", CellChangeTimes->{{3.40847384796875*^9, 3.408473853546875*^9}}], Cell[BoxData[ RowBox[{"GraphicsGrid", "[", RowBox[{ RowBox[{"{", RowBox[{"{", RowBox[{"TruSig", ",", "plt", ",", "SignalEst"}], "}"}], "}"}], ",", RowBox[{"ImageSize", "\[Rule]", "1000"}]}], "]"}]], "Input", CellChangeTimes->{{3.408473777875*^9, 3.408473821078125*^9}, { 3.40865231675*^9, 3.40865233215625*^9}, {3.408658425625*^9, 3.40865842571875*^9}}] }, Open ]], Cell[CellGroupData[{ Cell["Extra Credit Problem", "Section"], Cell[TextData[{ "Try denoising another of 4 basic functions. Simply replace the Heavisine \ function above with the function of your choice. All functions should be \ sampled on the interval [0,1]. Try changing the number of iteration of the \ transformation. You can also use ", StyleBox["Daub", FontColor->RGBColor[1, 0, 0]], "[6], ", StyleBox["Daub", FontColor->RGBColor[1, 0, 0]], "[8], ", StyleBox["Daub", FontColor->RGBColor[1, 0, 0]], "[10], ", StyleBox["Coif", FontColor->RGBColor[1, 0, 0]], "[1], ", StyleBox["Coif", FontColor->RGBColor[1, 0, 0]], "[2] as possible filters." }], "Text", CellChangeTimes->{{3.40847300734375*^9, 3.40847303146875*^9}, { 3.408658260984375*^9, 3.408658302203125*^9}, {3.408658429390625*^9, 3.40865845778125*^9}}], Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{"Doppler", "[", "t_", "]"}], ":=", RowBox[{ RowBox[{"Sqrt", "[", RowBox[{"t", "*", RowBox[{"(", RowBox[{"1", "-", "t"}], ")"}]}], "]"}], "*", RowBox[{"Sin", "[", RowBox[{"2", "*", "Pi", "*", RowBox[{"1.05", "/", RowBox[{"(", RowBox[{"t", "+", ".05"}], ")"}]}]}], "]"}]}]}], ";"}], "\n", RowBox[{ RowBox[{ RowBox[{"SquareWave", "[", "t_", "]"}], ":=", RowBox[{"Sign", "[", RowBox[{"Cos", "[", RowBox[{"8", "*", "Pi", "*", "t"}], "]"}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"Id", "[", "t_", "]"}], ":=", RowBox[{"Which", "[", RowBox[{ RowBox[{"0", "\[LessEqual]", "t", "\[LessEqual]", RowBox[{"1", "/", "4"}]}], ",", RowBox[{"t", "-", RowBox[{"1", "/", "8"}]}], ",", "True", ",", "0"}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"SawTooth", "[", "t_", "]"}], ":=", RowBox[{"Sum", "[", RowBox[{ RowBox[{"Id", "[", RowBox[{"t", "-", RowBox[{"k", "/", "4"}]}], "]"}], ",", RowBox[{"{", RowBox[{"k", ",", "0", ",", "3"}], "}"}]}], "]"}]}], ";"}]}], "Input", CellChangeTimes->{{3.4086523955625*^9, 3.40865240671875*^9}, { 3.4086524385*^9, 3.40865243884375*^9}, {3.40865842375*^9, 3.408658423921875*^9}}], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Doppler", "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "1"}], "}"}], ",", RowBox[{"PlotRange", "\[Rule]", "All"}]}], "]"}]], "Input", CellChangeTimes->{{3.4086524195*^9, 3.4086524989375*^9}, { 3.408658422015625*^9, 3.408658422234375*^9}}] }, Open ]], Cell[CellGroupData[{ Cell["Audio Denoising", "Section", CellChangeTimes->{{3.4217561678735423`*^9, 3.421756170154851*^9}}], Cell["\<\ Let's list the audio names and information. If the following cell fails, execute the cell below:\ \>", "Text", CellChangeTimes->{{3.4217582669272766`*^9, 3.4217582768806567`*^9}}, CellTags->"b:0.2.1"], Cell[BoxData[{ RowBox[{"AudioList", "[", "]"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"afiles", "=", RowBox[{"AudioNames", "[", "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"filename", "=", RowBox[{"afiles", "[", RowBox[{"[", "1", "]"}], "]"}]}], ";"}]}], "Input", CellChangeTimes->{{3.4217582873027983`*^9, 3.421758306084529*^9}}], Cell["\<\ If the cell above fails, use the cell below to create the audio file name.\ \>", "Text", CellChangeTimes->{{3.421758308959603*^9, 3.42175832466313*^9}}], Cell[BoxData[{ RowBox[{ RowBox[{ "dirname", "=", "\"\\""}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"audiofile", "=", "\"\\""}], ";"}], "\[IndentingNewLine]", RowBox[{"(*", RowBox[{ RowBox[{"audiofile", "=", "\"\\""}], ";"}], "*)"}], "\[IndentingNewLine]", RowBox[{"(*", RowBox[{ RowBox[{"audiofile", "=", "\"\\""}], ";"}], "*)"}], "\[IndentingNewLine]", RowBox[{"(*", RowBox[{ RowBox[{"audiofile", "=", "\"\\""}], ";"}], "*)"}], "\[IndentingNewLine]", RowBox[{"(*", RowBox[{ RowBox[{"audiofile", "=", "\"\\""}], ";"}], "*)"}], "\[IndentingNewLine]", RowBox[{"(*", RowBox[{ RowBox[{"audiofile", "=", "\"\\""}], ";"}], "*)"}], "\[IndentingNewLine]", RowBox[{"(*", RowBox[{ RowBox[{"audiofile", "=", "\"\\""}], ";"}], "*)"}], "\[IndentingNewLine]", RowBox[{"(*", RowBox[{ RowBox[{"audiofile", "=", "\"\\""}], ";"}], "*)"}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{"filename", "=", RowBox[{"dirname", "<>", "audiofile"}]}], ";"}]}], "Input", CellChangeTimes->{{3.4023988865*^9, 3.402398886703125*^9}, { 3.40293321892874*^9, 3.4029332190695133`*^9}, {3.42175645195894*^9, 3.421756461990447*^9}, {3.421756496960092*^9, 3.4217565569147515`*^9}, { 3.421756587478034*^9, 3.4217566628549633`*^9}, {3.4217575285958757`*^9, 3.4217575287208796`*^9}, {3.4217583321164455`*^9, 3.4217583331789727`*^9}}, CellTags->"b:0.2.1"], Cell[BoxData[ RowBox[{ RowBox[{"data", " ", "=", " ", RowBox[{"ReadSoundFile", "[", RowBox[{"filename", ",", RowBox[{"PrintHeader", "\[Rule]", "True"}]}], "]"}]}], ";"}]], "Input", CellChangeTimes->{{3.40239888846875*^9, 3.402398888609375*^9}, { 3.4029332175835705`*^9, 3.402933218037174*^9}, {3.421756427598941*^9, 3.4217564359741554`*^9}, {3.4217566699957714`*^9, 3.4217566730896006`*^9}, { 3.4217575265645742`*^9, 3.421757526642701*^9}}, CellTags->"b:0.2.1"], Cell[BoxData[{ RowBox[{ RowBox[{"data", "=", RowBox[{"ChopVector", "[", RowBox[{"data", ",", RowBox[{"PowersOfTwo", "\[Rule]", "4"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Print", "[", RowBox[{"\"\\"", ",", RowBox[{"Length", "[", "data", "]"}], ",", "\"\<.\>\""}], "]"}], ";"}]}], "Input", CellChangeTimes->{{3.40239888846875*^9, 3.402398888609375*^9}, { 3.4029332175835705`*^9, 3.402933218037174*^9}, 3.421756427598941*^9, { 3.4217575257051773`*^9, 3.4217575258458056`*^9}}, CellTags->"b:0.2.1"], Cell["The cell below plays the clip.", "Text", CellTags->"b:0.2.1"], Cell[BoxData[{ RowBox[{ RowBox[{"srate", "=", "22050"}], ";"}], "\[IndentingNewLine]", RowBox[{"ListPlay", "[", RowBox[{"data", ",", RowBox[{"SampleRate", "\[Rule]", "srate"}]}], "]"}]}], "Input", CellChangeTimes->{ 3.402327795953125*^9, {3.40239889065625*^9, 3.402398890921875*^9}, { 3.40293319922045*^9, 3.402933200768958*^9}, {3.421757523595748*^9, 3.4217575236895003`*^9}}, CellTags->"b:0.2.1"], Cell["\<\ We now generate some Gaussian white noise and add it to the sound clip.\ \>", "Text", CellTags->"b:0.2.1"], Cell[BoxData[{ RowBox[{ RowBox[{"n", "=", RowBox[{"Length", "[", "data", "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"nd", "=", RowBox[{"NormalDistribution", "[", RowBox[{"0.", ",", "1."}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"SeedRandom", "[", "]"}], ";"}], "\n", RowBox[{ RowBox[{"noise", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"Random", "[", "nd", "]"}], ",", RowBox[{"{", RowBox[{"k", ",", "1", ",", "n"}], "}"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"sigma", "=", "1000."}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"v", "=", RowBox[{"data", "+", RowBox[{"sigma", "*", "noise"}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{"ListPlay", "[", RowBox[{"v", ",", RowBox[{"SampleRate", "\[Rule]", "srate"}]}], "]"}]}], "Input", CellChangeTimes->{ 3.40232780490625*^9, {3.402398892640625*^9, 3.402398892765625*^9}, { 3.4029331981411867`*^9, 3.4029331982193937`*^9}, {3.4217575223300905`*^9, 3.4217575223925924`*^9}}, CellTags->"b:0.2.1"], Cell["\<\ We next compute four iterations of the wavelet transform. We'll use the D6 \ filter.\ \>", "Text", CellChangeTimes->{3.4217562091871*^9}, CellTags->"b:0.2.1"], Cell[BoxData[{ RowBox[{ RowBox[{"its", "=", "4"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"wt", "=", RowBox[{"WT1D", "[", RowBox[{ RowBox[{"N", "[", "v", "]"}], ",", RowBox[{"N", "[", RowBox[{"Daub", "[", "6", "]"}], "]"}], ",", RowBox[{"NumIterations", "\[Rule]", "its"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{"WaveletVectorPlot", "[", RowBox[{"wt", ",", RowBox[{"NumIterations", "\[Rule]", "its"}], ",", RowBox[{"PlotRange", "\[Rule]", "All"}], ",", RowBox[{"PlotLabel", "->", "\"\\""}], ",", RowBox[{"UseColors", "->", "True"}]}], "]"}]}], "Input", CellChangeTimes->{{3.40232781109375*^9, 3.4023278406875*^9}, { 3.402398894375*^9, 3.402398894515625*^9}, {3.4029326863000326`*^9, 3.40293268919051*^9}, {3.4029331965770364`*^9, 3.4029331966239605`*^9}, { 3.4217562124371834`*^9, 3.4217562180935783`*^9}, {3.4217575207675505`*^9, 3.4217575208613033`*^9}}, CellTags->"b:0.2.1"], Cell["\<\ Our first task is to estimate the noise level. Of course, we know the noise \ level but it is good to see how the estimator works.\ \>", "Text", CellTags->"b:0.2.1"], Cell[BoxData[ RowBox[{ RowBox[{"wtlist", "=", RowBox[{"WaveletVectorToList", "[", RowBox[{"wt", ",", RowBox[{"NumIterations", "\[Rule]", "its"}]}], "]"}]}], ";"}]], "Input", CellChangeTimes->{{3.402398896609375*^9, 3.4023988968125*^9}, { 3.4029326954089413`*^9, 3.4029326955026865`*^9}, {3.402933194371585*^9, 3.4029331944341507`*^9}, {3.421757519251887*^9, 3.421757519345639*^9}}, CellTags->"b:0.2.1"], Cell[BoxData[{ RowBox[{ RowBox[{"hp1", "=", RowBox[{"Last", "[", "wtlist", "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"sigmahat", "=", RowBox[{ RowBox[{"MAD", "[", "hp1", "]"}], "/", ".6745"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Print", "[", RowBox[{ "\"\\"", ",", "sigmahat", ",", "\"\<.\>\""}], "]"}], ";"}]}], "Input", CellChangeTimes->{{3.402398897828125*^9, 3.40239889796875*^9}, { 3.4029331931046233`*^9, 3.402933193214114*^9}, {3.421757516673696*^9, 3.4217575167830734`*^9}}, CellTags->"b:0.2.1"], Cell[TextData[{ "We will use Donoho's SUREShrink method for picking a different tolerance \ \[Lambda] for each highpass portion. \n\nThe variable ", StyleBox["lambda", FontWeight->"Bold"], " below is a vector if length its. The first element is the tolerance that \ is used for the first highpass portion, and so on. \n\nFeel free to pick \ your own ", StyleBox["lambda", FontWeight->"Bold"], " vector." }], "Text", CellChangeTimes->{{3.421756872094695*^9, 3.4217570075356627`*^9}}, CellTags->"b:0.2.1"], Cell[BoxData[{ RowBox[{ RowBox[{"lambda", "=", RowBox[{"sigmahat", "*", RowBox[{"Reverse", "[", RowBox[{"Map", "[", RowBox[{"DonohoSure", ",", "hp"}], "]"}], "]"}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"Print", "[", RowBox[{ "\"\\"", ",", "k", ",", "\"\< is \>\"", ",", RowBox[{"lambda", "[", RowBox[{"[", "k", "]"}], "]"}], ",", "\"\<.\>\""}], "]"}], ",", RowBox[{"{", RowBox[{"k", ",", "1", ",", "its"}], "}"}]}], "]"}], ";"}]}], "Input", CellChangeTimes->{{3.402327934375*^9, 3.402327934734375*^9}, { 3.402398902296875*^9, 3.4023989025*^9}, {3.4029327327664037`*^9, 3.4029327328913975`*^9}, {3.4029331900388894`*^9, 3.402933190117097*^9}, { 3.4217575154861655`*^9, 3.4217575155642924`*^9}}, CellTags->"b:0.2.1"], Cell["\<\ We now apply the shrinkage function (with the appropriate tolerance) to each \ highpass portion of the transform.\ \>", "Text", CellTags->"b:0.2.1"], Cell[BoxData[{ RowBox[{ RowBox[{"newhp", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"Map", "[", RowBox[{ RowBox[{ RowBox[{"ShrinkageFunction", "[", RowBox[{"#", ",", RowBox[{"lambda", "[", RowBox[{"[", "k", "]"}], "]"}]}], "]"}], "&"}], ",", RowBox[{"hp", "[", RowBox[{"[", "k", "]"}], "]"}]}], "]"}], ",", RowBox[{"{", RowBox[{"k", ",", "1", ",", "its"}], "}"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"lp", "=", RowBox[{"First", "[", "wtlist", "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"newwt", "=", RowBox[{"Join", "[", RowBox[{"lp", ",", RowBox[{"Flatten", "[", "newhp", "]"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{"WaveletVectorPlot", "[", RowBox[{"wt", ",", RowBox[{"NumIterations", "\[Rule]", "its"}], ",", RowBox[{"PlotRange", "\[Rule]", "All"}], ",", RowBox[{"PlotLabel", "->", "\"\\""}], ",", RowBox[{"UseColors", "->", "True"}]}], "]"}], "\[IndentingNewLine]", RowBox[{"WaveletVectorPlot", "[", RowBox[{"newwt", ",", RowBox[{"NumIterations", "\[Rule]", "its"}], ",", RowBox[{"PlotRange", "\[Rule]", "All"}], ",", RowBox[{"PlotLabel", "->", "\"\\""}], ",", RowBox[{"UseColors", "->", "True"}]}], "]"}]}], "Input", CellChangeTimes->{{3.40232789015625*^9, 3.402327898765625*^9}, { 3.40232793159375*^9, 3.402327938890625*^9}, {3.402398943328125*^9, 3.402398943609375*^9}, {3.40293273595374*^9, 3.4029327654834785`*^9}, { 3.4029331887250032`*^9, 3.4029331888501353`*^9}, {3.4217575141267557`*^9, 3.4217575142517586`*^9}}, CellTags->"b:0.2.1"], Cell["\<\ Finally, we apply the inverse wavelet transform and play the denoised audio \ clip.\ \>", "Text", CellTags->"b:0.2.1"], Cell[BoxData[ RowBox[{ RowBox[{"denoised", "=", RowBox[{"IWT1D", "[", RowBox[{"newwt", ",", RowBox[{"Daub", "[", "6", "]"}], ",", RowBox[{"NumIterations", "\[Rule]", "its"}]}], "]"}]}], ";"}]], "Input", CellChangeTimes->{{3.40232792953125*^9, 3.402327930015625*^9}, { 3.402398945375*^9, 3.40239894553125*^9}, {3.402932775545463*^9, 3.4029327776703544`*^9}, {3.402933186769816*^9, 3.4029331868793063`*^9}, { 3.421757512970476*^9, 3.421757513079854*^9}, {3.421758210441456*^9, 3.421758214410308*^9}}, CellTags->"b:0.2.1"], Cell["Here is the original noisy clip:", "Text", CellTags->"b:0.2.1"], Cell[BoxData[ RowBox[{"ListPlay", "[", RowBox[{"v", ",", RowBox[{"SampleRate", "\[Rule]", "srate"}]}], "]"}]], "Input", CellChangeTimes->{{3.40232790690625*^9, 3.40232792846875*^9}, { 3.4023989475625*^9, 3.402398948046875*^9}, {3.402933185643628*^9, 3.402933185737477*^9}, {3.421757512064203*^9, 3.421757512189206*^9}}, CellTags->"b:0.2.1"], Cell["Here is the denoised clip.", "Text", CellTags->"b:0.2.1"], Cell[BoxData[ RowBox[{"ListPlay", "[", RowBox[{"denoised", ",", RowBox[{"SampleRate", "\[Rule]", "srate"}]}], "]"}]], "Input", CellChangeTimes->{{3.40232791184375*^9, 3.402327925859375*^9}, { 3.402398950234375*^9, 3.4023989505*^9}, {3.4029331845487227`*^9, 3.4029331846894965`*^9}, {3.4217575108766727`*^9, 3.4217575111110535`*^9}}, CellTags->"b:0.2.1"] }, Open ]], Cell[CellGroupData[{ Cell["Exercise - For Fun", "Section", CellChangeTimes->{{3.4217570299581113`*^9, 3.421757033661331*^9}}], Cell[TextData[{ "Open the Laurel and Hardy audio file above. It is very large so you might \ you might want to use the ", StyleBox["Take", FontWeight->"Bold"], " command to only take say one quarter of the file. Play the audio file - \ it is noisy since it is an old recording.\n\nYour task: By whatever means \ possible (but hopefully using wavelets!), denoise the audio clip!" }], "Text", CellChangeTimes->{{3.4217572501981244`*^9, 3.421757272729951*^9}, { 3.421757303293234*^9, 3.4217574432186904`*^9}}] }, Open ]] }, Open ]] }, WindowSize->{1272, 902}, WindowMargins->{{-3, Automatic}, {Automatic, 2}}, FrontEndVersion->"6.0 for Microsoft Windows (32-bit) (April 20, 2007)", StyleDefinitions->FrontEnd`FileName[{"Creative"}, "NaturalColor.nb", CharacterEncoding -> "WindowsANSI"] ] (* End of Notebook Content *) (* Internal cache information *) (*CellTagsOutline CellTagsIndex->{ "b:0.2.1"->{ Cell[18707, 612, 213, 6, 65, "Text", CellTags->"b:0.2.1"], Cell[19459, 636, 1684, 47, 222, "Input", CellTags->"b:0.2.1"], Cell[21146, 685, 488, 10, 41, "Input", CellTags->"b:0.2.1"], Cell[21637, 697, 597, 15, 62, "Input", CellTags->"b:0.2.1"], Cell[22237, 714, 68, 1, 29, "Text", CellTags->"b:0.2.1"], Cell[22308, 717, 421, 10, 62, "Input", CellTags->"b:0.2.1"], Cell[22732, 729, 117, 3, 29, "Text", CellTags->"b:0.2.1"], Cell[22852, 734, 1093, 31, 162, "Input", CellTags->"b:0.2.1"], Cell[23948, 767, 171, 5, 29, "Text", CellTags->"b:0.2.1"], Cell[24122, 774, 992, 23, 82, "Input", CellTags->"b:0.2.1"], Cell[25117, 799, 177, 4, 29, "Text", CellTags->"b:0.2.1"], Cell[25297, 805, 427, 9, 41, "Input", CellTags->"b:0.2.1"], Cell[25727, 816, 614, 17, 82, "Input", CellTags->"b:0.2.1"], Cell[26344, 835, 517, 13, 101, "Text", CellTags->"b:0.2.1"], Cell[26864, 850, 896, 23, 62, "Input", CellTags->"b:0.2.1"], Cell[27763, 875, 159, 4, 29, "Text", CellTags->"b:0.2.1"], Cell[27925, 881, 1718, 43, 122, "Input", CellTags->"b:0.2.1"], Cell[29646, 926, 129, 4, 29, "Text", CellTags->"b:0.2.1"], Cell[29778, 932, 555, 12, 41, "Input", CellTags->"b:0.2.1"], Cell[30336, 946, 70, 1, 29, "Text", CellTags->"b:0.2.1"], Cell[30409, 949, 355, 7, 41, "Input", CellTags->"b:0.2.1"], Cell[30767, 958, 64, 1, 29, "Text", CellTags->"b:0.2.1"], Cell[30834, 961, 368, 7, 41, "Input", CellTags->"b:0.2.1"]} } *) (*CellTagsIndex CellTagsIndex->{ {"b:0.2.1", 32261, 1001} } *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[590, 23, 110, 1, 73, "Title"], Cell[703, 26, 149, 5, 96, "Subtitle"], Cell[855, 33, 158, 6, 89, "Subsubtitle"], Cell[CellGroupData[{ Cell[1038, 43, 29, 0, 75, "Section"], Cell[1070, 45, 163, 4, 29, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[1270, 54, 30, 0, 75, "Section"], Cell[1303, 56, 1259, 28, 101, "Text"] }, Open ]], Cell[CellGroupData[{ Cell[2599, 89, 145, 2, 75, "Section"], Cell[2747, 93, 391, 7, 62, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[3175, 105, 41, 0, 75, "Section"], Cell[3219, 107, 298, 9, 29, "Text"], Cell[3520, 118, 876, 27, 41, "Input"], Cell[4399, 147, 42, 0, 29, "Text"], Cell[4444, 149, 685, 21, 102, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[5166, 175, 35, 0, 75, "Section"], Cell[5204, 177, 144, 4, 47, "Text"], Cell[5351, 183, 468, 13, 41, "Input"], Cell[5822, 198, 191, 8, 29, "Text"], Cell[6016, 208, 675, 18, 96, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[6728, 231, 33, 0, 75, "Section"], Cell[6764, 233, 117, 1, 29, "Text"], Cell[6884, 236, 577, 16, 102, "Input"], Cell[7464, 254, 272, 8, 29, "Text"], Cell[7739, 264, 219, 5, 41, "Input"], Cell[7961, 271, 131, 3, 29, "Text"], Cell[8095, 276, 844, 25, 82, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[8976, 306, 45, 0, 75, "Section"], Cell[9024, 308, 413, 10, 29, "Text"], Cell[9440, 320, 1000, 31, 142, "Input"], Cell[10443, 353, 575, 13, 47, "Text"], Cell[11021, 368, 640, 16, 102, "Input"], Cell[11664, 386, 42, 0, 29, "Text"], Cell[11709, 388, 459, 14, 62, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[12205, 407, 29, 0, 75, "Section"], Cell[12237, 409, 214, 4, 29, "Text"], Cell[12454, 415, 2445, 71, 302, "Input"], Cell[14902, 488, 91, 2, 29, "Text"], Cell[14996, 492, 493, 13, 62, "Input"], Cell[15492, 507, 109, 1, 29, "Text"], Cell[15604, 510, 381, 9, 41, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[16022, 524, 39, 0, 75, "Section"], Cell[16064, 526, 781, 23, 47, "Text"], Cell[16848, 551, 1382, 43, 102, "Input"], Cell[18233, 596, 332, 8, 41, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[18602, 609, 102, 1, 75, "Section"], Cell[18707, 612, 213, 6, 65, "Text", CellTags->"b:0.2.1"], Cell[18923, 620, 367, 9, 82, "Input"], Cell[19293, 631, 163, 3, 29, "Text"], Cell[19459, 636, 1684, 47, 222, "Input", CellTags->"b:0.2.1"], Cell[21146, 685, 488, 10, 41, "Input", CellTags->"b:0.2.1"], Cell[21637, 697, 597, 15, 62, "Input", CellTags->"b:0.2.1"], Cell[22237, 714, 68, 1, 29, "Text", CellTags->"b:0.2.1"], Cell[22308, 717, 421, 10, 62, "Input", CellTags->"b:0.2.1"], Cell[22732, 729, 117, 3, 29, "Text", CellTags->"b:0.2.1"], Cell[22852, 734, 1093, 31, 162, "Input", CellTags->"b:0.2.1"], Cell[23948, 767, 171, 5, 29, "Text", CellTags->"b:0.2.1"], Cell[24122, 774, 992, 23, 82, "Input", CellTags->"b:0.2.1"], Cell[25117, 799, 177, 4, 29, "Text", CellTags->"b:0.2.1"], Cell[25297, 805, 427, 9, 41, "Input", CellTags->"b:0.2.1"], Cell[25727, 816, 614, 17, 82, "Input", CellTags->"b:0.2.1"], Cell[26344, 835, 517, 13, 101, "Text", CellTags->"b:0.2.1"], Cell[26864, 850, 896, 23, 62, "Input", CellTags->"b:0.2.1"], Cell[27763, 875, 159, 4, 29, "Text", CellTags->"b:0.2.1"], Cell[27925, 881, 1718, 43, 122, "Input", CellTags->"b:0.2.1"], Cell[29646, 926, 129, 4, 29, "Text", CellTags->"b:0.2.1"], Cell[29778, 932, 555, 12, 41, "Input", CellTags->"b:0.2.1"], Cell[30336, 946, 70, 1, 29, "Text", CellTags->"b:0.2.1"], Cell[30409, 949, 355, 7, 41, "Input", CellTags->"b:0.2.1"], Cell[30767, 958, 64, 1, 29, "Text", CellTags->"b:0.2.1"], Cell[30834, 961, 368, 7, 41, "Input", CellTags->"b:0.2.1"] }, Open ]], Cell[CellGroupData[{ Cell[31239, 973, 105, 1, 75, "Section"], Cell[31347, 976, 518, 10, 65, "Text"] }, Open ]] }, Open ]] } ] *) (* End of internal cache information *)