(*^ ::[ Information = "This is a Mathematica Notebook file. It contains ASCII text, and can be transferred by email, ftp, or other text-file transfer utility. It should be read or edited using a copy of Mathematica or MathReader. If you received this as email, use your mail application or copy/paste to save everything from the line containing (*^ down to the line containing ^*) into a plain text file. On some systems you may have to give the file a name ending with ".ma" to allow Mathematica to recognize it as a Notebook. The line below identifies what version of Mathematica created this file, but it can be opened using any other version as well."; FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2"; MacintoshStandardFontEncoding; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8, 24, "Times"; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "Times"; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6, 14, "Times"; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 18, "Times"; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 14, "Times"; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 12, "Times"; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-2, 12, "Courier"; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-2, 12, "Courier"; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R32768, L-2, 12, "Courier"; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-2, 12, "Courier"; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B32768, L-2, 12, "Courier"; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier"; fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, 10, "Geneva"; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = leftheader, inactive, L2, 12, "Times"; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; fontset = leftfooter, inactive, L2, 12, "Times"; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; paletteColors = 128; currentKernel; ] :[font = title; inactive; preserveAspect; startGroup] Linear Methods of Applied Mathematics Filtering and signal processing ;[s] 3:0,0;38,1;39,0;71,-1; 2:2,25,18,Times,1,24,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = text; inactive; preserveAspect; plain; bold] (c) Copyright 1994-1997 by Evans M. Harrell II and James V. Herod. All rights reserved. :[font = text; inactive; Cclosed; preserveAspect; plain; bold; startGroup] Notes for the instructor. ;[s] 2:0,1;25,0;26,-1; 2:1,13,9,Times,1,12,0,0,0;1,11,8,Times,1,9,0,0,0; :[font = text; inactive; preserveAspect; endGroup] This contains calculations and examples which correlate with an appendix of the WWW text by Harrell and Herod. Students can be encouraged to cut and paste from this notebook to do homework. :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Instructions :[font = text; inactive; preserveAspect; endGroup] This notebook uses Mathematica to perform calculations for Harrell and Herod's hypertext book, Linear Methods of Applied Mathematics. The student needs only a basic knowledge of Mathematica to use the notebook, which is designed both to show how to work problems in the text and to provide a template for using Mathematica to work other problems of the student's own design. Calculations will be performed when the reader presses enter in a given calculation cell (bold print). It is best to activate the cells in order, so that Mathematica will be able to call on operators and functions defined in earlier cells. Red color coding is used to warn the reader when a given calculation relies on an earlier one. ;[s] 11:0,0;19,1;30,0;96,1;133,0;181,1;192,0;315,1;326,0;538,1;549,0;721,-1; 2:6,13,9,Times,0,12,0,0,0;5,13,9,Times,2,12,0,0,0; :[font = text; inactive; preserveAspect] An intuitive way to understand projections of functions, which has many applications in signal processing, is in terms of filtering. The word "filtering" refers to an attempt to extract the important part of some data while eliminating random contributions called "noise" or other unwanted features which obscure the ones that matter. To illustrate the procedure, we choose a simple selection of features that are regarded as good and bad, to create "low-pass" and "high-pass" filters. Suppose that we have a signal f(t) which contains some pure frequencies we are interested in, say corresponding to pure vibratory functions تتت sin(2¹t), sin(4¹t), ..., sin(20¹t), cos(2¹t), cos(4¹t), ..., cos(20¹t), تتتتتتتتتت (f1) with frequencies 1,2,..., 10, plus some noise fnoise(t). (The frequency of a sinusoidal function sin(t) differs from the angular frequency by the factor 2 .) We write the full function as: تتت f(t) = a1 cos(2 t) + a2 cos(4 t) + ... a10 cos(20 t) تتتتتتتت + b1 sin(2 t) + b2 sin(4 t) + ... b10 sin(20 t) تتتتتتتت + fnoise(t) تتتتتت =: ffiltered(t) + fnoise(t) Often the constant term (a0 in the sections on Fourier series) is left out of the expansion, since a signal usually oscillates around average zero. Example F.1 Let us use Mathematica to see the effect of filtering as applied to some interesting functions on the interval 0 < x < 1. We begin with the function sin(1/(x-.05)), which oscillates strongly near x=0: ;[s] 27:0,0;773,1;778,0;930,1;931,0;944,1;945,0;962,1;964,0;988,1;989,0;1002,1;1003,0;1020,1;1022,0;1046,1;1051,0;1067,1;1075,0;1082,1;1087,0;1118,1;1119,0;1241,2;1252,0;1264,3;1275,0;1455,-1; 4:14,13,9,Times,0,12,0,0,0;11,21,13,Times,64,12,0,0,0;1,13,9,Times,1,12,0,0,0;1,13,9,Times,2,12,0,0,0; :[font = input; preserveAspect] f[x_] := Sin[1/(x - .05)] :[font = input; preserveAspect] Plot[Sin[1/(x - .05)], {x, 0, 1}] :[font = text; inactive; preserveAspect] The integrals to calculate the Fourier series for this function cannot be done in closed form, so we calculate them numerically, by slightly modifying the code we used in the notebook for Chapter II: :[font = input; dontPreserveAspect] Ave[g_, {a_,b_}] := (1/(b-a)) NIntegrate[(g /. x -> intvar1), \ {intvar1, a, b}] A[g_,m_, {a_,b_}] := (2/(b-a)) NIntegrate[ \ (Cos[2 m Pi intvar2/(b-a)] g /. x -> intvar2), \ {intvar2, a, b}] B[g_,n_, {a_,b_}] := (2/(b-a)) NIntegrate[ \ (Sin[2 n Pi intvar3/(b-a)] g /. x -> intvar3), \ {intvar3, a, b}] :[font = input; preserveAspect] LowPassFilter[f_,N_,{a_,b_}] := Ave[f,{a,b}] + \ Sum[Cos[2 n Pi x/(b-a)] A[f,n,{a,b}], {n,1,N}] + \ Sum[Sin[2 n Pi x/(b-a)] B[f,n,{a,b}], {n,1,N}] ;[s] 5:0,0;84,1;85,0;139,1;140,0;162,-1; 2:3,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = input; Cclosed; preserveAspect; startGroup] LowPassFilter[f[x],3, {0,1}] ;[s] 4:0,1;13,0;14,1;15,0;29,-1; 2:2,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = message; inactive; preserveAspect] NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 7 recursive bisections in intvar1 near intvar1 = 0.0664062. :[font = message; inactive; preserveAspect] NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 7 recursive bisections in intvar2 near intvar2 = 0.0664062. :[font = message; inactive; preserveAspect] NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 7 recursive bisections in intvar2 near intvar2 = 0.0664062. :[font = message; inactive; preserveAspect] General::stop: Further output of NIntegrate::ncvb will be suppressed during this calculation. :[font = output; output; inactive; preserveAspect; endGroup] 0.4587191403209957 - 0.0804973066364554*Cos[2*Pi*x] + 0.1923333905677756*Cos[4*Pi*x] - 0.1738189362009336*Cos[6*Pi*x] - 0.7232503640527357*Sin[2*Pi*x] + 0.03135682128579393*Sin[4*Pi*x] + 0.005857476439925394*Sin[6*Pi*x] ;[o] 0.458719 - 0.0804973 Cos[2 Pi x] + 0.192333 Cos[4 Pi x] - 0.173819 Cos[6 Pi x] - 0.72325 Sin[2 Pi x] + 0.0313568 Sin[4 Pi x] + 0.00585748 Sin[6 Pi x] :[font = text; inactive; preserveAspect] We notice that there were convergence problems associated with the wildness of this function. We shall ignore them and assume that the calculations are sufficiently accurate. :[font = input; preserveAspect] pic1 = Plot[%, {x,0,1}, PlotStyle -> RGBColor[1,0,0]] pic2 = Plot[f[x], {x,0,1}] Show[{pic1,pic2}] ;[s] 7:0,0;12,1;13,0;89,1;93,0;94,1;98,0;101,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; preserveAspect] For comparison, let's look at the high-frequency part of this function. A high-pass filter would throw away the result of the low-pass filter and keep the rest: :[font = input; preserveAspect] HighPassFilter[g_,N_, {a_,b_}] := g - LowPassFilter[g,N, {a,b}] ;[s] 3:0,0;38,1;51,0;64,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = input; Cclosed; preserveAspect; startGroup] HighPassFilter[f[x], 3, {0,1}] :[font = message; inactive; preserveAspect] NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 7 recursive bisections in intvar1 near intvar1 = 0.0664062. :[font = message; inactive; preserveAspect] NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 7 recursive bisections in intvar2 near intvar2 = 0.0664062. :[font = message; inactive; preserveAspect] NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 7 recursive bisections in intvar2 near intvar2 = 0.0664062. :[font = message; inactive; preserveAspect] General::stop: Further output of NIntegrate::ncvb will be suppressed during this calculation. :[font = output; output; inactive; preserveAspect; endGroup] -0.4587191403209957 + 0.0804973066364554*Cos[2*Pi*x] - 0.1923333905677756*Cos[4*Pi*x] + 0.1738189362009336*Cos[6*Pi*x] + Sin[(-0.05 + x)^(-1)] + 0.7232503640527357*Sin[2*Pi*x] - 0.03135682128579393*Sin[4*Pi*x] - 0.005857476439925394*Sin[6*Pi*x] ;[o] -0.458719 + 0.0804973 Cos[2 Pi x] - 0.192333 Cos[4 Pi x] + 1 0.173819 Cos[6 Pi x] + Sin[---------] + -0.05 + x 0.72325 Sin[2 Pi x] - 0.0313568 Sin[4 Pi x] - 0.00585748 Sin[6 Pi x] :[font = input; preserveAspect] pic3 = Plot[%, {x,0,1}, PlotStyle -> RGBColor[1,0,0]] Show[{pic2,pic3}] ;[s] 7:0,0;12,1;13,0;61,1;65,0;66,1;70,0;73,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; preserveAspect] Example F.2 This time we perform signal processing on the function sin(1/(x-.05)), in order to enhance the low frequencies, with the integral kernel: ;[s] 2:0,1;11,0;152,-1; 2:1,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] IKernel[x_,t_] := 6 Sum[Cos[2 Pi n x] Cos[2 Pi n t] /n, {n,1,10}] + \ 6 Sum[Sin[2 Pi m x] Sin[2 Pi m t] /m, {m,1,10}] :[font = input; preserveAspect] EnhanceLow[g_] := 6 Sum[Cos[2 Pi n x] NIntegrate[(g /. x -> t) Cos[2 Pi n t], {t,0,1}] /n, {n,1,10}] + \ 6 Sum[Sin[2 Pi m x] NIntegrate[(g /. x -> t) Sin[2 Pi m t], {t,0,1}] /m, {m,1,10}] :[font = input; Cclosed; preserveAspect; startGroup] EnhanceLow[f[x]] ;[s] 4:0,1;10,0;11,1;12,0;17,-1; 2:2,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = message; inactive; preserveAspect] NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 7 recursive bisections in t near t = 0.0664062. :[font = message; inactive; preserveAspect] NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 7 recursive bisections in t near t = 0.0664062. :[font = message; inactive; preserveAspect] NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 7 recursive bisections in t near t = 0.0664062. :[font = message; inactive; preserveAspect] General::stop: Further output of NIntegrate::ncvb will be suppressed during this calculation. :[font = message; inactive; preserveAspect] NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, oscillatory integrand, or insufficient WorkingPrecision. :[font = output; output; inactive; preserveAspect; endGroup] 6*(-0.04024865331822768*Cos[2*Pi*x] + 0.04808334764194389*Cos[4*Pi*x] - 0.02896982270015559*Cos[6*Pi*x] - 0.0092612346534609*Cos[8*Pi*x] + 0.01148150974601692*Cos[10*Pi*x] + 0.006593230035384617*Cos[12*Pi*x] - 0.004515971681640377*Cos[14*Pi*x] - 0.006398706836345607*Cos[16*Pi*x] - 0.001280388080518789*Cos[18*Pi*x] + 0.002865763840281879*Cos[20*Pi*x]) + 6*(-0.3616251820263678*Sin[2*Pi*x] + 0.007839205321448484*Sin[4*Pi*x] + 0.000976246073320899*Sin[6*Pi*x] - 0.02674924034044754*Sin[8*Pi*x] - 0.01261375492307642*Sin[10*Pi*x] + 0.003354205014612682*Sin[12*Pi*x] + 0.003522670160376309*Sin[14*Pi*x] - 0.003249359079018527*Sin[16*Pi*x] - 0.006085396178419817*Sin[18*Pi*x] - 0.00391887721688457*Sin[20*Pi*x]) ;[o] 6 (-0.0402487 Cos[2 Pi x] + 0.0480833 Cos[4 Pi x] - 0.0289698 Cos[6 Pi x] - 0.00926123 Cos[8 Pi x] + 0.0114815 Cos[10 Pi x] + 0.00659323 Cos[12 Pi x] - 0.00451597 Cos[14 Pi x] - 0.00639871 Cos[16 Pi x] - 0.00128039 Cos[18 Pi x] + 0.00286576 Cos[20 Pi x]) + 6 (-0.361625 Sin[2 Pi x] + 0.00783921 Sin[4 Pi x] + 0.000976246 Sin[6 Pi x] - 0.0267492 Sin[8 Pi x] - 0.0126138 Sin[10 Pi x] + 0.00335421 Sin[12 Pi x] + 0.00352267 Sin[14 Pi x] - 0.00324936 Sin[16 Pi x] - 0.0060854 Sin[18 Pi x] - 0.00391888 Sin[20 Pi x]) :[font = subsubsection; inactive; preserveAspect; startGroup] Tip in case you have trouble with this command :[font = text; inactive; preserveAspect; endGroup] This involved some heavy-duty integration. If you do not have a powerful computer, you may prefer to have fewer than 20 terms in this calculation. You can do this by modifying the definition of EnhanceLow. ;[s] 3:0,0;196,1;206,0;208,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,10,Courier,0,12,0,0,0; :[font = input; preserveAspect; endGroup] pic4 = Plot[%, {x,0,1}, PlotStyle -> RGBColor[1,0,0]] Show[{pic2,pic4}] ;[s] 7:0,0;12,1;13,0;61,1;65,0;66,1;70,0;73,-1; 2:4,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,0; ^*)