%23%20%2F%2F%2F%20script%0A%23%20requires-python%20%3D%20%22%3E%3D3.11%22%0A%23%20dependencies%20%3D%20%5B%22marimo%3E%3D0.23.9%22%2C%20%22numpy%3E%3D2.0%22%2C%20%22matplotlib%3E%3D3.9%22%2C%20%22pandas%3E%3D2.2%22%2C%20%22scipy%3E%3D1.14%22%2C%20%22urllib3%3E%3D2.2%22%5D%0A%23%20%2F%2F%2F%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.25.2%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Data%20filtering%20in%20signal%20processing%0A%0A%20%20%20%20%3E%20Marcos%20Duarte%2C%20Renato%20Naville%20Watanabe%3Cbr%3E%0A%20%20%20%20%3E%20%5BLaboratory%20of%20Biomechanics%20and%20Motor%20Control%5D(https%3A%2F%2Fbmclab.pesquisa.ufabc.edu.br%2F)%3Cbr%3E%0A%20%20%20%20%3E%20Federal%20University%20of%20ABC%2C%20Brazil%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20How%20to%20use%20this%20guide%0A%0A%20%20%20%20This%20notebook%20is%20an%20introduction%20to%20data%20filtering%20and%20to%20the%20most%20basic%20filters%20used%20to%20process%20biomechanical%20data%3A%20the%20moving%20average%2C%20the%20Butterworth%20filter%2C%20the%20critically%20damped%20filter%20and%20a%20few%20others.%20It%20starts%20from%20a%20concrete%20problem%2C%20getting%20an%20acceleration%20out%20of%20measured%20positions%2C%20and%20keeps%20coming%20back%20to%20it%20until%20the%20problem%20is%20solved%20reasonably%20well.%0A%0A%20%20%20%20You%20should%20be%20familiar%20with%20the%20%5Bbasic%20properties%20of%20signals%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fsignal_basic_properties.py)%20before%20proceeding%3A%20frequency%2C%20harmonics%2C%20sampling%2C%20noise%20and%20the%20signal-to-noise%20ratio.%20Finite%20differences%2C%20as%20used%20in%20%5BKinematics%20of%20a%20particle%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fkinematics_particle.py)%2C%20will%20also%20come%20up.%0A%0A%20%20%20%20Read%20it%20in%20order%20and%20run%20each%20cell%20as%20you%20reach%20it.%20Work%20through%20the%20distributed%20Exercises%20before%20opening%20their%20solutions.%20The%20final%20Problems%20have%20collapsed%20final%20answers%20for%20independent%20practice.%20Where%20you%20find%20a%20**Challenge**%20or%20a%20set%20of%20**Guiding%20questions**%2C%20stop%20and%20answer%20on%20a%20scratchpad%20before%20moving%20on.%20Several%20of%20them%20ask%20you%20to%20predict%20a%20number%20*before*%20the%20code%20prints%20it%3B%20the%20prediction%20is%20the%20point%2C%20and%20being%20wrong%20is%20the%20most%20useful%20thing%20that%20can%20happen%20to%20you%20here.%0A%0A%20%20%20%20**Challenge%200.**%20Think%20of%20a%20quantity%20you%20would%20like%20to%20know%20but%20can%20only%20get%20by%20differentiating%20a%20measured%20position%3A%20the%20acceleration%20of%20a%20body%20segment%2C%20the%20velocity%20of%20a%20marker%2C%20the%20jerk%20of%20a%20reaching%20movement.%20If%20the%20position%20has%20an%20error%20of%201%20mm%2C%20guess%20how%20large%20the%20error%20in%20the%20acceleration%20will%20be.%20Write%20the%20guess%20down%3B%20you%20will%20be%20able%20to%20check%20its%20order%20of%20magnitude%20on%20real%20data%20in%20a%20moment.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Python%20setup%0A%0A%20%20%20%20NumPy%20and%20SciPy%20do%20the%20computation%20(SciPy's%20%60signal%60%20module%20has%20the%20filters)%2C%20pandas%20reads%20the%20data%20files%2C%20and%20Matplotlib%20draws%20the%20plots.%20These%20packages%20are%20declared%20in%20the%20notebook%20metadata.%20Open%20the%20file%20in%20marimo%3B%20the%20Pezzack%2C%20EMG%20and%20ball-toss%20examples%20need%20access%20to%20the%20repository%20data.%20No%20motion-capture%20hardware%20or%20video%20software%20is%20needed%20to%20run%20the%20supplied%20examples.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_()%3A%0A%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20from%20notebook_resources%20import%20repo_file%0A%20%20%20%20except%20ModuleNotFoundError%20as%20_error%3A%0A%20%20%20%20%20%20%20%20if%20_error.name%20!%3D%20%22notebook_resources%22%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20raise%0A%20%20%20%20%20%20%20%20raise%20ModuleNotFoundError(%0A%20%20%20%20%20%20%20%20%20%20%20%20%22This%20notebook%20needs%20notebook_resources.py.%20Download%20it%20from%3A%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22https%3A%2F%2Fraw.githubusercontent.com%2FBMClab%2FBMC%2Fmaster%2Ffunctions%2Fnotebook_resources.py%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22Save%20it%20in%20the%20same%20folder%20as%20this%20notebook%2C%20then%20rerun%20this%20cell.%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20name%3D%22notebook_resources%22%2C%0A%20%20%20%20%20%20%20%20)%20from%20None%0A%20%20%20%20return%20(repo_file%2C)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20timeit%0A%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20pandas%20as%20pd%0A%20%20%20%20import%20matplotlib%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20from%20scipy%20import%20interpolate%2C%20signal%0A%0A%20%20%20%20matplotlib.rc(%22axes%22%2C%20labelsize%3D13%2C%20titlesize%3D14)%0A%20%20%20%20matplotlib.rc(%22xtick%22%2C%20labelsize%3D11)%0A%20%20%20%20matplotlib.rc(%22ytick%22%2C%20labelsize%3D11)%0A%20%20%20%20matplotlib.rc(%22legend%22%2C%20fontsize%3D11)%0A%20%20%20%20return%20interpolate%2C%20np%2C%20pd%2C%20plt%2C%20signal%2C%20timeit%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20The%20problem%3A%20an%20acceleration%20from%20measured%20positions%0A%0A%20%20%20%20In%201977%2C%20Pezzack%2C%20Norman%20and%20Winter%20published%20a%20paper%20investigating%20the%20effects%20of%20differentiation%20and%20filtering%20on%20experimental%20data.%20They%20filmed%20a%20bar%20being%20rotated%20by%20hand%20in%20a%20horizontal%20plane%2C%20digitized%20its%20angle%20from%20the%20film%2C%20and%20at%20the%20same%20time%20measured%20its%20angular%20acceleration%20directly%2C%20with%20an%20accelerometer.%20Since%20then%20these%20data%20have%20become%20a%20benchmark%20for%20testing%20new%20algorithms%20(they%20are%20also%20available%20from%20the%20%5BISB%20website%5D(https%3A%2F%2Fisbweb.org%2Fdata%2Fpezzack%2Findex.html)).%20We%20will%20use%20the%20accelerometer's%20acceleration%20as%20an%20independent%20reference%2C%20recognizing%20that%20it%20too%20has%20measurement%20uncertainty%2C%20and%20ask%3A%20can%20we%20recover%20it%20from%20the%20angle%20alone%3F%0A%0A%20%20%20%20The%20most%20direct%20answer%20is%20to%20differentiate%20the%20angle%20twice%20with%20finite%20differences.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20how%20the%20double-differentiated%20angle%20will%20compare%20with%20the%20measured%20acceleration.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20pd%2C%20plt%2C%20repo_file)%3A%0A%20%20%20%20_data%20%3D%20pd.read_csv(%0A%20%20%20%20%20%20%20%20repo_file('data%2FPezzack.txt')%2C%0A%20%20%20%20%20%20%20%20sep%3D%22%5Ct%22%2C%0A%20%20%20%20%20%20%20%20header%3DNone%2C%0A%20%20%20%20%20%20%20%20skiprows%3D6%2C%0A%20%20%20%20).to_numpy()%0A%20%20%20%20pz_time%2C%20pz_disp%2C%20pz_disp_noisy%2C%20pz_aacc%20%3D%20_data.T%0A%20%20%20%20pz_dt%20%3D%20np.mean(np.diff(pz_time))%0A%0A%20%20%20%20%23%20acceleration%20by%20the%20second-order%20finite%20difference%20(2%20samples%20shorter)%0A%20%20%20%20pz_aacc_diff%20%3D%20np.diff(pz_disp%2C%202)%20%2F%20pz_dt**2%0A%20%20%20%20_rmse%20%3D%20np.sqrt(np.mean((pz_aacc_diff%20-%20pz_aacc%5B1%3A-1%5D)%20**%202))%0A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(11%2C%204))%0A%20%20%20%20_axs%5B0%5D.plot(pz_time%2C%20pz_disp%2C%20%22b.-%22)%0A%20%20%20%20_axs%5B0%5D.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_axs%5B0%5D.set_ylabel(%22Angular%20displacement%20%5Brad%5D%22)%0A%20%20%20%20_axs%5B1%5D.plot(pz_time%2C%20pz_aacc%2C%20%22g%22%2C%20linewidth%3D2%2C%20label%3D%22Accelerometer%20(reference)%22)%0A%20%20%20%20_axs%5B1%5D.plot(pz_time%5B1%3A-1%5D%2C%20pz_aacc_diff%2C%20%22r%22%2C%20label%3D%22Angle%20differentiated%20twice%22)%0A%20%20%20%20_axs%5B1%5D.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_axs%5B1%5D.set_ylabel(%22Angular%20acceleration%20%5Brad%2Fs%24%5E2%24%5D%22)%0A%20%20%20%20_axs%5B1%5D.legend(frameon%3DFalse%2C%20loc%3D%22upper%20left%22)%0A%20%20%20%20plt.suptitle(%22Pezzack's%20benchmark%20data%22%2C%20fontsize%3D16)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%0A%20%20%20%20print(f%22Sampling%20interval%3A%20%7Bpz_dt%3A.4f%7D%20s%20(%7B1%20%2F%20pz_dt%3A.1f%7D%20Hz)%2C%20%7Bpz_time.size%7D%20samples%22)%0A%20%20%20%20print(f%22RMSE%20of%20the%20differentiated%20acceleration%3A%20%7B_rmse%3A.1f%7D%20rad%2Fs2%22)%0A%20%20%20%20return%20pz_aacc%2C%20pz_disp%2C%20pz_dt%2C%20pz_time%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20angle%20looks%20clean%3A%20a%20smooth%20curve%20with%20no%20visible%20noise.%20Its%20second%20derivative%20is%20not.%20The%20red%20curve%20follows%20the%20reference%20acceleration%20only%20roughly%2C%20with%20a%20jagged%20error%20of%20tens%20of%20rad%2Fs%C2%B2%20in%20places%20and%20an%20RMSE%20(root-mean-square%20error)%20of%20about%2011%20rad%2Fs%C2%B2.%0A%0A%20%20%20%20The%20noise%20comes%20from%20small%20random%20errors%20in%20the%20digitization%20of%20each%20film%20frame%2C%20a%20fraction%20of%20a%20degree.%20Invisible%20in%20the%20angle%2C%20they%20dominate%20the%20acceleration.%20This%20notebook%20is%20about%20why%20that%20happens%20and%20what%20to%20do%20about%20it.%0A%0A%20%20%20%20**Guiding%20questions%200.**%0A%0A%20%20%20%201.%20The%20errors%20in%20the%20angle%20are%20far%20too%20small%20to%20see%20in%20the%20left%20panel.%20Why%20do%20they%20become%20so%20large%20in%20the%20right%20one%3F%0A%20%20%20%202.%20Is%20the%20error%20in%20the%20red%20curve%20larger%20where%20the%20acceleration%20is%20large%2C%20or%20is%20it%20spread%20over%20the%20whole%20recording%3F%0A%20%20%20%203.%20Would%20sampling%20the%20film%20faster%2C%20with%20the%20same%20digitization%20error%20per%20frame%2C%20make%20the%20acceleration%20better%20or%20worse%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Why%20differentiation%20amplifies%20noise%0A%0A%20%20%20%20Removing%20noise%20from%20a%20signal%20is%20rarely%20trivial%2C%20and%20the%20problem%20gets%20worse%20with%20numerical%20differentiation%3A%20noise%20at%20frequencies%20higher%20than%20those%20of%20the%20signal%20is%20amplified%20by%20differentiation%2C%20so%20the%20signal-to-noise%20ratio%20(SNR)%20drops%20with%20each%20derivative.%0A%0A%20%20%20%20To%20see%20why%2C%20consider%20the%20following%20function%20representing%20some%20experimental%20data%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20f%20%3D%20%5Csin(%5Comega%20t)%20%2B%200.1%5Csin(10%5Comega%20t)%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20first%20component%2C%20with%20a%20large%20amplitude%20(1)%20and%20a%20low%20frequency%20(1%20Hz)%2C%20represents%20the%20signal%3B%20the%20second%2C%20with%20a%20small%20amplitude%20(0.1)%20and%20a%20high%20frequency%20(10%20Hz)%2C%20represents%20the%20noise.%20The%20SNR%20of%20these%20data%20is%20%24(1%2F0.1)%5E2%20%3D%20100%24.%0A%0A%20%20%20%20**Before%20you%20read%20on**%2C%20differentiate%20%24f%24%20twice%20by%20hand%20and%20predict%20the%20SNR%20of%20%24f'%24%20and%20of%20%24f''%24.%0A%0A%20%20%20%20%24%24%0A%20%20%20%20f%5C%2C'%20%3D%20%5Comega%20%5Ccos(%5Comega%20t)%20%2B%20%5Comega%20%5Ccos(10%5Comega%20t)%0A%20%20%20%20%24%24%0A%0A%20%20%20%20%24%24%0A%20%20%20%20f%5C%2C''%20%3D%20-%5Comega%5E2%20%5Csin(%5Comega%20t)%20-%2010%5Comega%5E2%20%5Csin(10%5Comega%20t)%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Each%20differentiation%20multiplies%20a%20component%20of%20frequency%20%24f%24%20by%20%242%5Cpi%20f%24%2C%20so%20the%2010%20Hz%20noise%20gains%20a%20factor%20of%2010%20on%20the%201%20Hz%20signal%20every%20time.%20For%20the%20first%20derivative%2C%20SNR%20%3D%201%2C%20and%20for%20the%20second%2C%20SNR%20%3D%200.01%3A%20the%20noise%20now%20has%20a%20hundred%20times%20the%20power%20of%20the%20signal.%20The%20following%20plots%20illustrate%20the%20problem%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt)%3A%0A%20%20%20%20_t%20%3D%20np.arange(0%2C%201%2C%200.01)%0A%20%20%20%20_w%20%3D%202%20*%20np.pi%20*%201%20%20%23%201%20Hz%0A%20%20%20%20%23%20signal%20and%20noise%2C%20and%20their%20derivatives%0A%20%20%20%20_s%2C%20_n%20%3D%20np.sin(_w%20*%20_t)%2C%200.1%20*%20np.sin(10%20*%20_w%20*%20_t)%0A%20%20%20%20_sd%2C%20_nd%20%3D%20_w%20*%20np.cos(_w%20*%20_t)%2C%20_w%20*%20np.cos(10%20*%20_w%20*%20_t)%0A%20%20%20%20_sdd%2C%20_ndd%20%3D%20-_w%20*%20_w%20*%20np.sin(_w%20*%20_t)%2C%20-_w%20*%20_w%20*%2010%20*%20np.sin(10%20*%20_w%20*%20_t)%0A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(3%2C%201%2C%20sharex%3DTrue%2C%20figsize%3D(9%2C%206.5))%0A%20%20%20%20_axs%5B0%5D.set_title(%22Differentiation%20of%20signal%20and%20noise%22)%0A%20%20%20%20for%20_ax%2C%20_sig%2C%20_noi%2C%20_label%20in%20zip(%0A%20%20%20%20%20%20%20%20_axs%2C%20(_s%2C%20_sd%2C%20_sdd)%2C%20(_n%2C%20_nd%2C%20_ndd)%2C%20(%22f%22%2C%20%22f%20'%22%2C%20%22f%20''%22)%0A%20%20%20%20)%3A%0A%20%20%20%20%20%20%20%20_ax.plot(_t%2C%20_sig%2C%20%22b.-%22%2C%20linewidth%3D1%2C%20label%3D%22signal%22)%0A%20%20%20%20%20%20%20%20_ax.plot(_t%2C%20_noi%2C%20%22g.-%22%2C%20linewidth%3D1%2C%20label%3D%22noise%22)%0A%20%20%20%20%20%20%20%20_ax.plot(_t%2C%20_sig%20%2B%20_noi%2C%20%22r.-%22%2C%20linewidth%3D2%2C%20label%3D%22signal%2Bnoise%22)%0A%20%20%20%20%20%20%20%20_ax.set_ylabel(_label)%0A%20%20%20%20_axs%5B0%5D.legend(frameon%3DFalse%2C%20fontsize%3D10%2C%20loc%3D%22upper%20right%22)%0A%20%20%20%20_axs%5B2%5D.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20In%20the%20top%20panel%20the%20noise%20is%20a%20small%20ripple%20on%20the%20signal.%20In%20the%20bottom%20panel%20the%20signal%20is%20the%20small%20one.%20That%20is%20the%20Pezzack%20acceleration%20in%20miniature%2C%20and%20it%20points%20to%20the%20remedy%3A%20if%20the%20noise%20lives%20at%20higher%20frequencies%20than%20the%20signal%2C%20remove%20those%20frequencies%20*before*%20differentiating.%20That%20is%20what%20a%20low-pass%20filter%20does.%0A%0A%20%20%20%20%23%23%20Filter%20and%20smoothing%0A%0A%20%20%20%20In%20data%20acquisition%20with%20an%20instrument%2C%20the%20noise%20often%20has%20higher%20frequencies%20and%20lower%20amplitudes%20than%20the%20signal%20of%20interest.%20Removing%20it%20is%20called%20filtering%20or%20smoothing.%0A%0A%20%20%20%20%3Ca%20href%3D%22https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FFilter_(signal_processing)%22%3EFiltering%3C%2Fa%3E%20is%20a%20process%20that%20attenuates%20some%20unwanted%20component%20or%20feature%20of%20a%20signal.%20A%20filter%20usually%20removes%20certain%20frequency%20components%20of%20the%20data%2C%20according%20to%20its%20frequency%20response.%20The%20%5Bfrequency%20response%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FFrequency_response)%20is%20the%20quantitative%20measure%20of%20the%20output%20spectrum%20of%20a%20system%20in%20response%20to%20a%20stimulus%2C%20and%20is%20used%20to%20characterize%20the%20dynamics%20of%20the%20system.%0A%0A%20%20%20%20%5BSmoothing%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FSmoothing)%20is%20the%20removal%20of%20local%20(short-scale)%20fluctuations%20from%20the%20data%20while%20preserving%20a%20more%20global%20pattern%20(those%20local%20variations%20could%20be%20noise%2C%20or%20just%20a%20short-scale%20phenomenon%20that%20is%20not%20of%20interest).%20A%20filter%20with%20a%20low-pass%20frequency%20response%20performs%20smoothing.%0A%0A%20%20%20%20With%20respect%20to%20its%20implementation%2C%20a%20filter%20can%20be%20an%20%5Banalog%20filter%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FPassive_analogue_filter_development)%20or%20a%20%5Bdigital%20filter%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FDigital_filter).%20An%20analog%20filter%20is%20an%20electronic%20circuit%20that%20filters%20an%20input%20electrical%20signal%20and%20outputs%20a%20filtered%20electrical%20signal%3B%20a%20simple%20one%20can%20be%20built%20with%20a%20resistor%20and%20a%20capacitor.%20A%20digital%20filter%20is%20a%20system%20that%20filters%20digital%20(discrete-time)%20data.%20The%20anti-aliasing%20filter%20of%20an%20acquisition%20system%20must%20be%20analog%2C%20because%20it%20acts%20before%20sampling%3B%20everything%20in%20this%20notebook%20is%20digital%2C%20and%20acts%20on%20data%20already%20recorded.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Example%3A%20the%20moving-average%20filter%0A%0A%20%20%20%20A%20simple%20low-pass%20(smoothing)%20filter%20is%20the%20moving%20average%3A%20the%20arithmetic%20mean%20of%20successive%20subsequences%20of%20%24m%24%20terms%20of%20the%20data.%20For%20instance%2C%20the%20moving%20averages%20with%20window%20sizes%20%24m%24%20equal%20to%202%20and%203%20are%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bl%7D%0A%20%20%20%20y_%7BMA(2)%7D%20%3D%20%5Cfrac%7B1%7D%7B2%7D%5Bx_1%2Bx_2%2C%5C%3B%20x_2%2Bx_3%2C%5C%3B%20%5Ccdots%2C%5C%3B%20x_%7Bn-1%7D%2Bx_n%5D%20%5C%5C%0A%20%20%20%20y_%7BMA(3)%7D%20%3D%20%5Cfrac%7B1%7D%7B3%7D%5Bx_1%2Bx_2%2Bx_3%2C%5C%3B%20x_2%2Bx_3%2Bx_4%2C%5C%3B%20%5Ccdots%2C%5C%3B%20x_%7Bn-2%7D%2Bx_%7Bn-1%7D%2Bx_n%5D%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20with%20the%20general%20formula%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20y%5Bi%5D%20%3D%20%5Cfrac%7B1%7D%7Bm%7D%5Csum_%7Bj%3D0%7D%5E%7Bm-1%7D%20x%5Bi%2Bj%5D%20%5Cquad%20%5Ctext%7Bfor%7D%20%5Cquad%20i%3D1%2C%20%5C%3B%20%5Cdots%2C%20%5C%3B%20n-m%2B1%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24n%24%20is%20the%20number%20of%20data.%20Here%20is%20a%20naive%20implementation%2C%20with%20a%20loop%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.function%0Adef%20moving_average(x%2C%20window)%3A%0A%20%20%20%20%22%22%22Moving%20average%20of%20%60x%60%20with%20window%20size%20%60window%60%20(a%20naive%20loop).%22%22%22%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20y%20%3D%20np.empty(len(x)%20-%20window%20%2B%201)%0A%20%20%20%20for%20i%20in%20range(len(y))%3A%0A%20%20%20%20%20%20%20%20y%5Bi%5D%20%3D%20np.sum(x%5Bi%20%3A%20i%20%2B%20window%5D)%20%2F%20window%0A%20%20%20%20return%20y%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Let's%20test%20it%20on%20a%20signal%20that%20jumps%20from%200%20to%201%20and%20back%2C%20plus%20random%20noise%20with%20a%20standard%20deviation%20of%200.1.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20what%20a%20window%20of%2011%20samples%20will%20do%20to%20the%20noise%2C%20and%20what%20it%20will%20do%20to%20the%20two%20jumps.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt)%3A%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(seed%3D42)%0A%20%20%20%20_signal%20%3D%20np.zeros(300)%0A%20%20%20%20_signal%5B100%3A200%5D%20%2B%3D%201%0A%20%20%20%20_x%20%3D%20_signal%20%2B%20_rng.standard_normal(300)%20%2F%2010%0A%0A%20%20%20%20_y%20%3D%20moving_average(_x%2C%2011)%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(9%2C%204))%0A%20%20%20%20_ax.plot(_x%2C%20%22b.-%22%2C%20linewidth%3D1%2C%20label%3D%22raw%20data%22)%0A%20%20%20%20_ax.plot(_y%2C%20%22r.-%22%2C%20linewidth%3D2%2C%20label%3D%22moving%20average%22)%0A%20%20%20%20_ax.legend(frameon%3DFalse%2C%20loc%3D%22upper%20right%22)%0A%20%20%20%20_ax.set_xlabel(%22Sample%22)%0A%20%20%20%20_ax.set_ylabel(%22Amplitude%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22Standard%20deviation%20of%20the%20noise%2C%20before%3A%20%7Bnp.std(_x%5B%3A90%5D)%3A.3f%7D%2C%20%22%0A%20%20%20%20%20%20%20%20f%22after%3A%20%7Bnp.std(_y%5B%3A85%5D)%3A.3f%7D%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20noise%20is%20reduced%20several%20times.%20Averaging%20%24m%24%20independent%20values%20reduces%20their%20standard%20deviation%20by%20%24%5Csqrt%7Bm%7D%24%20on%20average%2C%20about%203.3%20for%20%24m%20%3D%2011%24%3B%20an%20estimate%20from%20a%20stretch%20of%20data%20this%20short%20scatters%20a%20good%20deal%20around%20that%20value.%20The%20price%20is%20paid%20at%20the%20jumps%3A%20each%20one%20is%20now%20a%20ramp%2011%20samples%20long.%20Every%20low-pass%20filter%20makes%20this%20trade%20between%20noise%20and%20sharpness%3B%20the%20rest%20of%20the%20notebook%20is%20about%20making%20it%20well.%0A%0A%20%20%20%20Look%20also%20at%20*where*%20the%20jumps%20are.%20The%20raw%20data%20jump%20at%20samples%20100%20and%20200%2C%20the%20filtered%20ones%20about%205%20samples%20earlier.%20The%20output%20%60y%5Bi%5D%60%20is%20the%20mean%20of%20%60x%5Bi%5D%60%20to%20%60x%5Bi%2B10%5D%60%2C%20centred%20on%20sample%20%60i%2B5%60%2C%20but%20we%20plotted%20it%20at%20%60i%60.%20We%20will%20come%20back%20to%20this%2C%20and%20to%20faster%20implementations%20of%20the%20moving%20average%2C%20later.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Exercise%202.**%20Apply%20a%20three-point%20moving%20average%20to%20%24%5B0%2C0%2C3%2C0%2C0%5D%24.%20What%20happens%20to%20a%20constant%20sequence%20of%203%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%202%0A%20%20%20%20The%20valid%20windows%20have%20sums%203%2C%203%20and%203%2C%20giving%20%24%5B1%2C1%2C1%5D%24.%20A%20constant%20sequence%20has%20window%20sums%209%2C%20giving%203%20at%20every%20complete%20window%3A%20the%20coefficients%20sum%20to%20one.%20The%20isolated%20peak%20is%20spread%20in%20time%2C%20explaining%20both%20the%20smoothing%20and%20the%20loss%20of%20sharpness.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Digital%20filters%0A%0A%20%20%20%20In%20signal%20processing%2C%20a%20digital%20filter%20is%20a%20system%20that%20performs%20mathematical%20operations%20on%20a%20signal%20to%20modify%20certain%20aspects%20of%20it.%20A%20digital%20filter%20(more%20precisely%2C%20a%20causal%2C%20linear%2C%20time-invariant%20(LTI)%20digital%20filter)%20can%20be%20seen%20as%20the%20implementation%20of%20the%20following%20difference%20equation%20in%20the%20time%20domain%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bl%20l%7D%0A%20%20%20%20y_n%20%26%3D%20b_0x_n%20%2B%20b_1x_%7Bn-1%7D%20%2B%20%5Ccdots%20%2B%20b_Mx_%7Bn-M%7D%20-%20a_1y_%7Bn-1%7D%20-%20%5Ccdots%20-%20a_Ny_%7Bn-N%7D%20%5C%5C%0A%20%20%20%20%26%3D%20%5Csum_%7Bk%3D0%7D%5EM%20b_kx_%7Bn-k%7D%20-%20%5Csum_%7Bk%3D1%7D%5EN%20a_ky_%7Bn-k%7D%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20the%20output%20%24y%24%20is%20the%20filtered%20version%20of%20the%20input%20%24x%24%2C%20%24a_k%24%20and%20%24b_k%24%20are%20the%20filter%20coefficients%20(real%20values)%2C%20and%20the%20order%20of%20the%20filter%20is%20the%20larger%20of%20%24N%24%20and%20%24M%24.%0A%0A%20%20%20%20This%20general%20equation%20describes%20a%20recursive%20filter%3A%20the%20filtered%20value%20%24y_n%24%20is%20computed%20from%20current%20and%20previous%20values%20of%20%24x%24%20*and*%20from%20previous%20values%20of%20%24y%24%2C%20its%20own%20output%3B%20it%20is%20a%20system%20with%20feedback.%20A%20filter%20that%20does%20not%20reuse%20its%20outputs%20as%20inputs%20(a%20system%20with%20only%20feedforward)%20is%20called%20nonrecursive%2C%20and%20its%20%24a%24%20coefficients%20are%20zero.%20Recursive%20and%20nonrecursive%20filters%20are%20also%20known%20as%20infinite%20impulse%20response%20(IIR)%20and%20finite%20impulse%20response%20(FIR)%20filters%2C%20respectively.%0A%0A%20%20%20%20A%20filter%20with%20only%20the%20terms%20on%20previous%20values%20of%20%24y%24%20is%20also%20known%20as%20an%20autoregressive%20(AR)%20filter%3B%20one%20with%20only%20the%20terms%20on%20current%20and%20previous%20values%20of%20%24x%24%2C%20as%20a%20moving-average%20(MA)%20filter%3B%20and%20one%20with%20both%2C%20as%20an%20autoregressive%20moving-average%20(ARMA)%20filter.%20The%20moving%20average%20of%20the%20previous%20section%20is%20an%20FIR%20filter%20with%20%24m%24%20coefficients%20%24b%24%2C%20each%20equal%20to%20%241%2Fm%24%2C%20and%20all%20%24a%24%20coefficients%20equal%20to%20zero.%0A%0A%20%20%20%20%23%23%23%20Transfer%20function%0A%0A%20%20%20%20Another%20way%20to%20characterize%20a%20filter%20is%20by%20its%20%5Btransfer%20function%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FTransfer_function).%20In%20simple%20terms%2C%20the%20transfer%20function%20is%20the%20ratio%2C%20in%20the%20frequency%20domain%2C%20between%20the%20output%20and%20the%20input%20signals%20of%20a%20filter.%20For%20a%20continuous-time%20input%20%24x(t)%24%20and%20output%20%24y(t)%24%2C%20the%20transfer%20function%20%24H(s)%24%20is%20the%20ratio%20between%20the%20%5BLaplace%20transforms%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FLaplace_transform)%20of%20the%20output%20and%20of%20the%20input%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20H(s)%20%3D%20%5Cfrac%7BY(s)%7D%7BX(s)%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24s%20%3D%20%5Csigma%20%2B%20j%5Comega%24%2C%20%24j%24%20is%20the%20imaginary%20unit%20and%20%24%5Comega%24%20is%20the%20angular%20frequency%2C%20%242%5Cpi%20f%24.%20In%20the%20steady-state%20response%20we%20can%20take%20%24%5Csigma%3D0%24%2C%20and%20the%20Laplace%20transforms%20with%20complex%20argument%20reduce%20to%20%5BFourier%20transforms%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FFourier_transform)%20with%20real%20argument%20%24%5Comega%24.%0A%0A%20%20%20%20For%20a%20discrete-time%20input%20and%20output%2C%20the%20transfer%20function%20%24H(z)%24%20is%20the%20ratio%20between%20their%20%5Bz-transforms%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FZ-transform)%2C%20and%20the%20formalism%20is%20similar.%20Taking%20the%20z-transform%20of%20the%20difference%20equation%20above%20gives%20the%20transfer%20function%20of%20a%20linear%2C%20time-invariant%2C%20causal%20digital%20filter%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20H(z)%20%3D%20%5Cfrac%7BY(z)%7D%7BX(z)%7D%20%3D%20%5Cfrac%7Bb_0%20%2B%20b_1%20z%5E%7B-1%7D%20%2B%20b_2%20z%5E%7B-2%7D%20%2B%20%5Ccdots%20%2B%20b_M%20z%5E%7B-M%7D%7D%7B1%20%2B%20a_1%20z%5E%7B-1%7D%20%2B%20a_2%20z%5E%7B-2%7D%20%2B%20%5Ccdots%20%2B%20a_N%20z%5E%7B-N%7D%7D%0A%20%20%20%20%3D%20%5Cfrac%7B%5Csum_%7Bk%3D0%7D%5EM%20b_kz%5E%7B-k%7D%7D%7B1%20%2B%20%5Csum_%7Bk%3D1%7D%5EN%20a_kz%5E%7B-k%7D%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%20again%20the%20order%20of%20the%20filter%20is%20the%20larger%20of%20%24N%24%20and%20%24M%24.%20As%20with%20the%20difference%20equation%2C%20this%20is%20the%20transfer%20function%20of%20a%20recursive%20(IIR)%20filter%3B%20if%20the%20%24a%24%20coefficients%20are%20zero%2C%20the%20denominator%20equals%20one%20and%20the%20filter%20is%20nonrecursive%20(FIR).%0A%0A%20%20%20%20%23%23%23%20The%20Fourier%20transform%0A%0A%20%20%20%20The%20%5BFourier%20transform%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FFourier_transform)%20is%20a%20mathematical%20operation%20that%20transforms%20a%20signal%20that%20is%20a%20function%20of%20time%2C%20%24g(t)%24%2C%20into%20a%20function%20of%20frequency%2C%20%24G(f)%24%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathcal%7BF%7D%5Bg(t)%5D%20%3D%20G(f)%20%3D%20%5Cint_%7B-%5Cinfty%7D%5E%7B%5Cinfty%7D%20g(t)%5C%2C%20e%5E%7B-j%202%5Cpi%20f%20t%7D%5C%2C%20dt%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Its%20inverse%20operation%20is%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathcal%7BF%7D%5E%7B-1%7D%5BG(f)%5D%20%3D%20g(t)%20%3D%20%5Cint_%7B-%5Cinfty%7D%5E%7B%5Cinfty%7D%20G(f)%5C%2C%20e%5E%7Bj%202%5Cpi%20f%20t%7D%5C%2C%20df%0A%20%20%20%20%24%24%0A%0A%20%20%20%20%24G(f)%24%20is%20the%20representation%20of%20the%20time-domain%20signal%20%24g(t)%24%20in%20the%20frequency%20domain%2C%20and%20vice%20versa%3B%20the%20two%20are%20called%20a%20Fourier%20transform%20pair.%20See%20the%20notebook%20%5BFourier%20transform%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Ffourier_transform.py)%2C%20or%20%5Bthis%20text%5D(https%3A%2F%2Fwww.thefouriertransform.com%2Ftransform%2Ffourier.php)%2C%20for%20an%20introduction.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Types%20of%20filters%0A%0A%20%20%20%20With%20respect%20to%20the%20frequencies%20that%20are%20*not*%20removed%20from%20the%20data%2C%20with%20a%20boundary%20given%20by%20the%20critical%20or%20cutoff%20frequency%2C%20a%20filter%20can%20be%20low-pass%2C%20high-pass%2C%20band-pass%20or%20band-stop.%20The%20frequency%20responses%20of%20these%20filters%20are%20illustrated%20in%20the%20next%20figure.%0A%0A%20%20%20%20%3Cfigure%20aria-describedby%3D%22figure-caption-1%22%3E%3Ccenter%3E%3Cimg%20src%3D%22https%3A%2F%2Fthumb.wikimedia.org%2Fwikipedia%2Fcommons%2Fthumb%2Fe%2Fec%2FBandform_template.svg%2F960px-Bandform_template.svg.png%22%20width%3D500%20alt%3D%22Filters%22%2F%3E%3C%2Fcenter%3E%3C%2Ffigure%3E%0A%0A%20%20%20%20%3Cspan%20id%3D%22figure-caption-1%22%20class%3D%22figure-caption%22%20style%3D%22display%3Ablock%3Btext-align%3Acenter%3Bfont-style%3Aitalic%22%3EFigure.%20Frequency%20response%20of%20filters%20(%3Ca%20href%3D%22https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FFilter_(signal_processing)%22%3Efrom%20Wikipedia%3C%2Fa%3E).%3C%2Fspan%3E%0A%0A%20%20%20%20The%20critical%20or%20cutoff%20frequency%20of%20a%20filter%20is%20the%20frequency%20at%20which%20the%20power%20(the%20amplitude%20squared)%20of%20the%20filtered%20signal%20is%20half%20the%20power%20of%20the%20input%20signal%2C%20or%2C%20equivalently%2C%20the%20output%20amplitude%20is%200.707%20times%20the%20input%20amplitude.%20For%20instance%2C%20if%20a%20low-pass%20filter%20has%20a%20cutoff%20frequency%20of%2010%20Hz%2C%20a%2010%20Hz%20component%20comes%20out%20with%2050%25%20of%20its%20power%20and%20about%2071%25%20of%20its%20amplitude.%20(The%20same%20%241%2F%5Csqrt%7B2%7D%24%20as%20the%20RMS%20of%20a%20sinusoid%2C%20in%20%5BBasic%20properties%20of%20signals%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fsignal_basic_properties.py).)%0A%0A%20%20%20%20The%20gain%20of%20a%20filter%20(the%20ratio%20between%20the%20output%20and%20input%20powers)%20is%20usually%20expressed%20in%20decibels%20(dB).%0A%0A%20%20%20%20%23%23%23%20Decibel%20(dB)%0A%0A%20%20%20%20The%20%3Ca%20href%3D%22https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FDecibel%22%3Edecibel%20(dB)%3C%2Fa%3E%20is%20a%20logarithmic%20unit%20used%20to%20express%20the%20ratio%20between%20two%20values.%20For%20the%20gain%20of%20a%20filter%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20Gain%3D10%5C%2C%5Clog_%7B10%7D%5Cleft(%5Cfrac%7BA_%7Bout%7D%5E2%7D%7BA_%7Bin%7D%5E2%7D%5Cright)%3D20%5C%2C%5Clog_%7B10%7D%5Cleft(%5Cfrac%7BA_%7Bout%7D%7D%7BA_%7Bin%7D%7D%5Cright)%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24A_%7Bout%7D%24%20and%20%24A_%7Bin%7D%24%20are%20the%20amplitudes%20of%20the%20output%20(filtered)%20and%20input%20(raw)%20signals.%20A%20decibel%20is%20one%20tenth%20of%20a%20bel%2C%20a%20unit%20named%20in%20honor%20of%20%3Ca%20href%3D%22https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FAlexander_Graham_Bell%22%3EAlexander%20Graham%20Bell%3C%2Fa%3E.%0A%0A%20%20%20%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20the%20gain%20in%20decibels%20at%20the%20cutoff%20frequency%2C%20where%20the%20power%20is%20halved%2C%20and%20for%20an%20output%20amplitude%20attenuated%20by%20a%20factor%20of%2010%20and%20of%201000.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20print(%22Amplitude%20ratio%20%20%20Power%20ratio%20%20%20Gain%20%5BdB%5D%22)%0A%20%20%20%20for%20_ratio%20in%20(10%2C%202%2C%201%2C%201%20%2F%20np.sqrt(2)%2C%200.5%2C%200.1%2C%200.001)%3A%0A%20%20%20%20%20%20%20%20print(f%22%7B_ratio%3A15.4f%7D%20%7B_ratio**2%3A13.6f%7D%20%7B20%20*%20np.log10(_ratio)%3A11.1f%7D%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20cutoff%20frequency%20is%20the%20**%24-3%24%20dB**%20point%3A%20%2410%5Clog_%7B10%7D(0.5)%20%5Capprox%20-3%24%20dB.%20Doubling%20the%20power%20is%20%24%2B3%24%20dB.%20Every%20factor%20of%20ten%20in%20amplitude%20is%2020%20dB%2C%20so%20an%20attenuation%20of%20the%20amplitude%20by%201000%20times%20is%20%24-60%24%20dB%2C%20and%20%24-120%24%20dB%20would%20be%20a%20million%20times.%20That%20compression%20is%20why%20the%20decibel%20is%20useful%20for%20describing%20filters%2C%20whose%20attenuation%20spans%20many%20orders%20of%20magnitude.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Butterworth%20filter%0A%0A%20%20%20%20A%20common%20filter%20in%20biomechanics%20and%20motor%20control%20is%20the%20%5BButterworth%20filter%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FButterworth_filter).%20It%20is%20popular%20because%20it%20is%20simple%20to%20design%20and%20to%20use%2C%20and%20because%20its%20frequency%20response%20is%20maximally%20flat%20in%20the%20pass%20band%3A%20it%20does%20not%20ripple%2C%20so%20the%20frequencies%20it%20keeps%20are%20kept%20with%20nearly%20the%20same%20gain.%20It%20is%20a%20recursive%20(IIR)%20filter%2C%20and%20both%20its%20%24a%24%20and%20%24b%24%20coefficients%20are%20used.%0A%0A%20%20%20%20In%20SciPy%2C%20the%20function%20%60butter%60%20calculates%20the%20filter%20coefficients%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20butter(N%2C%20Wn%2C%20btype%3D'low'%2C%20analog%3DFalse%2C%20output%3D'ba')%0A%20%20%20%20%60%60%60%0A%0A%20%20%20%20where%20%60N%60%20is%20the%20order%20of%20the%20filter%2C%20%60Wn%60%20is%20the%20cutoff%20frequency%20given%20as%20a%20fraction%20of%20the%20%5BNyquist%20frequency%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FNyquist_frequency)%20(half%20the%20sampling%20frequency)%2C%20and%20%60btype%60%20is%20the%20type%20of%20filter%20(one%20of%20%60'lowpass'%60%2C%20%60'highpass'%60%2C%20%60'bandpass'%60%2C%20%60'bandstop'%60%3B%20the%20default%20is%20%60'lowpass'%60).%20The%20filtering%20itself%20is%20done%20by%20the%20function%20%60lfilter%60%2C%20which%20implements%20the%20difference%20equation%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20lfilter(b%2C%20a%2C%20x%2C%20axis%3D-1%2C%20zi%3DNone)%0A%20%20%20%20%60%60%60%0A%0A%20%20%20%20where%20%60b%60%20and%20%60a%60%20are%20the%20coefficients%20calculated%20by%20%60butter%60%2C%20and%20%60x%60%20is%20the%20data%20to%20be%20filtered.%0A%0A%20%20%20%20Let's%20filter%20the%20signal-plus-noise%20function%20from%20the%20section%20on%20differentiation%2C%20sampled%20at%20100%20Hz%2C%20with%20a%20second-order%20Butterworth%20low-pass%20filter%20at%205%20Hz%2C%20between%20the%201%20Hz%20signal%20and%20the%2010%20Hz%20noise.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20signal)%3A%0A%20%20%20%20fs%20%3D%20100%20%20%23%20sampling%20frequency%20%5BHz%5D%0A%20%20%20%20bw_t%20%3D%20np.arange(0%2C%201%2C%201%20%2F%20fs)%0A%20%20%20%20_w%20%3D%202%20*%20np.pi%20*%201%20%20%23%201%20Hz%0A%20%20%20%20bw_y%20%3D%20np.sin(_w%20*%20bw_t)%20%2B%200.1%20*%20np.sin(10%20*%20_w%20*%20bw_t)%0A%0A%20%20%20%20%23%20second-order%20Butterworth%20low-pass%20filter%20at%205%20Hz%0A%20%20%20%20_b%2C%20_a%20%3D%20signal.butter(2%2C%205%20%2F%20(fs%20%2F%202)%2C%20btype%3D%22low%22)%0A%20%20%20%20bw_y2%20%3D%20signal.lfilter(_b%2C%20_a%2C%20bw_y)%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(9%2C%204))%0A%20%20%20%20_ax.plot(bw_t%2C%20bw_y%2C%20%22r.-%22%2C%20linewidth%3D2%2C%20label%3D%22raw%20data%22)%0A%20%20%20%20_ax.plot(bw_t%2C%20bw_y2%2C%20%22b.-%22%2C%20linewidth%3D2%2C%20label%3D%22filter%20%40%205%20Hz%22)%0A%20%20%20%20_ax.legend(frameon%3DFalse)%0A%20%20%20%20_ax.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22Amplitude%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20bw_t%2C%20bw_y%2C%20bw_y2%2C%20fs%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20ripple%20is%20gone.%20But%20the%20filtered%20signal%20is%20also%20late%3A%20the%20Butterworth%20filter%20introduces%20a%20phase%20shift%2C%20a%20delay%2C%20between%20the%20raw%20and%20filtered%20signals.%20We%20will%20measure%20it%2C%20and%20then%20remove%20it.%0A%0A%20%20%20%20First%2C%20a%20property%20of%20the%20coefficients.%20For%20a%20low-pass%20filter%20to%20leave%20a%20constant%20signal%20unchanged%20(a%20gain%20of%20one%20at%20zero%20frequency)%2C%20its%20transfer%20function%20at%20%24z%20%3D%201%24%20must%20equal%20one%2C%20and%20from%20the%20expression%20for%20%24H(z)%24%20this%20means%20that%20the%20sum%20of%20the%20%24b%24%20coefficients%20minus%20the%20sum%20of%20the%20%24a%24%20coefficients%2C%20excluding%20the%20first%2C%20is%20one%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20signal)%3A%0A%20%20%20%20print(%22Low-pass%20Butterworth%20filter%20coefficients%20(cutoff%20at%200.1%20of%20the%20Nyquist%20frequency)%22)%0A%20%20%20%20for%20_order%20in%20(1%2C%202)%3A%0A%20%20%20%20%20%20%20%20_b%2C%20_a%20%3D%20signal.butter(_order%2C%200.1%2C%20btype%3D%22low%22)%0A%20%20%20%20%20%20%20%20print(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22Order%20%7B_order%7D%3A%5Cn%20%20b%3A%20%7B_b%7D%5Cn%20%20a%3A%20%7B_a%7D%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22%20%20sum(b)%20-%20sum(a%5B1%3A%5D)%20%3D%20%7Bnp.sum(_b)%20-%20np.sum(_a%5B1%3A%5D)%3A.6f%7D%22%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Bode%20plot%0A%0A%20%20%20%20How%20much%20the%20amplitude%20of%20the%20filtered%20signal%20is%20attenuated%20relative%20to%20the%20raw%20signal%2C%20as%20a%20function%20of%20frequency%2C%20is%20given%20by%20the%20frequency%20response.%20The%20plots%20of%20the%20magnitude%20and%20phase%20responses%20form%20the%20%5BBode%20plot%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FBode_plot).%20Here%20it%20is%20for%20the%20filter%20we%20just%20used%20(Butterworth%2C%20low-pass%20at%205%20Hz%2C%20second%20order)%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(fs%2C%20np%2C%20plt%2C%20signal)%3A%0A%20%20%20%20_b%2C%20_a%20%3D%20signal.butter(2%2C%205%20%2F%20(fs%20%2F%202)%2C%20btype%3D%22low%22)%0A%20%20%20%20_f%2C%20_h%20%3D%20signal.freqz(_b%2C%20_a%2C%20worN%3D512%2C%20fs%3Dfs)%20%20%23%20frequency%20response%0A%20%20%20%20_angles%20%3D%20np.rad2deg(np.unwrap(np.angle(_h)))%20%20%23%20phase%20%5Bdegrees%5D%0A%20%20%20%20_h_db%20%3D%2020%20*%20np.log10(np.abs(_h))%20%20%23%20magnitude%20%5BdB%5D%0A%0A%20%20%20%20_fig%2C%20(_ax1%2C%20_ax2)%20%3D%20plt.subplots(2%2C%201%2C%20sharex%3DTrue%2C%20figsize%3D(9%2C%206.5))%0A%20%20%20%20_ax1.plot(_f%2C%20_h_db%2C%20linewidth%3D2)%0A%20%20%20%20_ax1.plot(5%2C%20-3.01%2C%20%22ro%22)%0A%20%20%20%20_ax1.set_ylim(-80%2C%203)%0A%20%20%20%20_ax1.set_title(%22Frequency%20response%22)%0A%20%20%20%20_ax1.set_ylabel(%22Magnitude%20%5BdB%5D%22)%0A%20%20%20%20_axi%20%3D%20_ax1.inset_axes(%5B0.06%2C%200.08%2C%200.33%2C%200.5%5D)%20%20%23%20zoom%20around%20the%20cutoff%0A%20%20%20%20_axi.plot(_f%2C%20_h_db%2C%20linewidth%3D2)%0A%20%20%20%20_axi.plot(5%2C%20-3.01%2C%20%22ro%22)%0A%20%20%20%20_axi.set_xlim(%5B0%2C%2010%5D)%0A%20%20%20%20_axi.set_ylim(%5B-6%2C%200.5%5D)%0A%20%20%20%20_axi.grid(True%2C%20linestyle%3D%22%3A%22)%0A%20%20%20%20_ax2.plot(_f%2C%20_angles%2C%20linewidth%3D2)%0A%20%20%20%20_ax2.plot(5%2C%20-90%2C%20%22ro%22)%0A%20%20%20%20_ax2.set_title(%22Phase%20response%22)%0A%20%20%20%20_ax2.set_xlabel(%22Frequency%20%5BHz%5D%22)%0A%20%20%20%20_ax2.set_ylabel(%22Phase%20%5Bdegrees%5D%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20inset%20shows%20that%20at%20the%20cutoff%20frequency%2C%205%20Hz%2C%20the%20filtered%20signal%20is%20indeed%20attenuated%20by%203%20dB.%20The%20phase%20response%20shows%20that%20at%20the%20cutoff%20frequency%2C%20the%20output%20lags%20the%20input%20by%20%2490%5Eo%24.%20A%205%20Hz%20component%20has%20a%20period%20of%200.2%20s%2C%20so%20%2490%5Eo%24%2C%20a%20quarter%20of%20a%20period%2C%20is%20a%20delay%20of%200.05%20s.%0A%0A%20%20%20%20But%20the%20delay%20we%20saw%20in%20the%20plot%20was%20of%20the%201%20Hz%20signal%2C%20not%20of%20a%205%20Hz%20component.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20the%20delay%20of%20the%201%20Hz%20component.%20Is%20it%200.05%20s%2C%20five%20times%20longer%2C%20or%20something%20else%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(fs%2C%20np%2C%20signal)%3A%0A%20%20%20%20_b%2C%20_a%20%3D%20signal.butter(2%2C%205%20%2F%20(fs%20%2F%202)%2C%20btype%3D%22low%22)%0A%20%20%20%20_freqs%20%3D%20np.array(%5B1%2C%202%2C%205%2C%2010%5D)%0A%20%20%20%20_%2C%20_h%20%3D%20signal.freqz(_b%2C%20_a%2C%20worN%3D_freqs%2C%20fs%3Dfs)%0A%20%20%20%20_phase%20%3D%20np.angle(_h)%0A%20%20%20%20print(%22Frequency%20%5BHz%5D%20%20%20Gain%20%5BdB%5D%20%20%20Phase%20%5Bdegrees%5D%20%20%20Delay%20%5Bs%5D%22)%0A%20%20%20%20for%20_f%2C%20_hh%2C%20_ph%20in%20zip(_freqs%2C%20_h%2C%20_phase)%3A%0A%20%20%20%20%20%20%20%20print(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22%7B_f%3A14d%7D%20%7B20%20*%20np.log10(np.abs(_hh))%3A11.2f%7D%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22%7Bnp.rad2deg(_ph)%3A17.1f%7D%20%7B-_ph%20%2F%20(2%20*%20np.pi%20*%20_f)%3A11.3f%7D%22%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%201%20Hz%20component%20lags%20by%20only%20%2416%5Eo%24%2C%20but%20%2416%5Eo%24%20of%20a%201%20s%20period%20is%200.045%20s%2C%20almost%20the%20same%20delay%20as%20the%205%20Hz%20component.%20In%20its%20pass%20band%2C%20the%20phase%20of%20the%20Butterworth%20filter%20grows%20nearly%20in%20proportion%20to%20the%20frequency%2C%20so%20every%20component%20the%20filter%20keeps%20is%20delayed%20by%20roughly%20the%20same%20time%2C%20about%20%24%5Csqrt%7B2%7D%2F(2%5Cpi%20f_c)%20%5Capprox%200.045%24%20s%20for%20this%20second-order%20filter.%20The%20filtered%20signal%20is%20close%20to%20a%20delayed%20copy%20of%20the%20slow%20part%20of%20the%20raw%20one.%0A%0A%20%20%20%20A%20delay%20of%2045%20ms%20is%20not%20a%20detail.%20If%20you%20filter%20a%20kinematic%20signal%20this%20way%20and%20compare%20it%20with%20an%20EMG%20or%20a%20force%20filtered%20differently%2C%20every%20timing%20between%20them%20is%20off%20by%20that%20much.%0A%0A%20%20%20%20%23%23%23%20Order%20of%20a%20filter%0A%0A%20%20%20%20The%20order%20of%20a%20filter%20determines%20how%20steep%20the%20%22wall%22%20of%20the%20frequency%20response%20is%20around%20the%20cutoff%20frequency.%20A%20vertical%20wall%20exactly%20at%20the%20cutoff%20would%20be%20ideal%2C%20but%20it%20is%20impossible%20to%20implement.%0A%0A%20%20%20%20A%20first-order%20Butterworth%20filter%20attenuates%206%20dB%20for%20each%20doubling%20of%20the%20frequency%20(per%20octave)%2C%20or%2C%20which%20is%20the%20same%2C%2020%20dB%20each%20time%20the%20frequency%20is%20multiplied%20by%2010%20(per%20decade).%20In%20technical%20terms%2C%20a%20first-order%20filter%20rolls%20off%20at%20%24-6%24%20dB%20per%20octave%2C%20or%20%24-20%24%20dB%20per%20decade.%20A%20second-order%20filter%20rolls%20off%20at%20%24-12%24%20dB%20per%20octave%20(%24-40%24%20dB%20per%20decade)%2C%20and%20so%20on%2C%20as%20the%20next%20figure%20shows.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.function%0Adef%20butterworth_plot()%3A%0A%20%20%20%20%22%22%22Plot%20the%20frequency%20response%20of%20Butterworth%20filters%20of%20different%20orders.%22%22%22%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20from%20scipy%20import%20signal%0A%0A%20%20%20%20fig%2C%20ax%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(11%2C%204.2))%0A%20%20%20%20colors%20%3D%20%7B1%3A%20%22b%22%2C%202%3A%20%22r%22%2C%204%3A%20%22g%22%2C%206%3A%20%22y%22%7D%0A%20%20%20%20for%20order%2C%20color%20in%20colors.items()%3A%0A%20%20%20%20%20%20%20%20b%2C%20a%20%3D%20signal.butter(order%2C%2010%2C%20%22low%22%2C%20analog%3DTrue)%0A%20%20%20%20%20%20%20%20w%2C%20h%20%3D%20signal.freqs(b%2C%20a%2C%20worN%3Dnp.logspace(0%2C%202%2C%20500))%0A%20%20%20%20%20%20%20%20ax%5B0%5D.plot(w%20%2F%2010%2C%2020%20*%20np.log10(abs(h))%2C%20color%2C%20linewidth%3D2)%0A%20%20%20%20%20%20%20%20ax%5B1%5D.plot(w%20%2F%2010%2C%20np.rad2deg(np.unwrap(np.angle(h)))%2C%20color%2C%20linewidth%3D2)%0A%20%20%20%20for%20axi%20in%20ax%3A%0A%20%20%20%20%20%20%20%20axi.axvline(1%2C%20color%3D%22black%22)%20%20%23%20cutoff%20frequency%0A%20%20%20%20%20%20%20%20axi.set_xscale(%22log%22)%0A%20%20%20%20%20%20%20%20axi.set_xlabel(%22Frequency%20%2F%20Critical%20frequency%22)%0A%20%20%20%20%20%20%20%20axi.grid(which%3D%22both%22%2C%20axis%3D%22both%22)%0A%20%20%20%20ax%5B0%5D.scatter(1%2C%20-3%2C%20marker%3D%22s%22%2C%20edgecolor%3D%220%22%2C%20facecolor%3D%221%22%2C%20s%3D400)%0A%20%20%20%20ax%5B0%5D.set_ylabel(%22Magnitude%20%5BdB%5D%22)%0A%20%20%20%20ax%5B0%5D.set_xlim(0.1%2C%2010)%0A%20%20%20%20ax%5B0%5D.set_ylim(-120%2C%2010)%0A%20%20%20%20ax%5B1%5D.legend(%5Bstr(order)%20for%20order%20in%20colors%5D%2C%20title%3D%22Filter%20order%22%2C%20loc%3D%22best%22)%0A%20%20%20%20ax%5B1%5D.set_ylabel(%22Phase%20%5Bdegrees%5D%22)%0A%20%20%20%20ax%5B1%5D.set_yticks(np.arange(0%2C%20-300%2C%20-45))%0A%20%20%20%20ax%5B1%5D.set_ylim(-300%2C%2010)%0A%20%20%20%20axi%20%3D%20ax%5B0%5D.inset_axes(%5B0.08%2C%200.08%2C%200.38%2C%200.45%5D)%20%20%23%20zoom%20around%20the%20cutoff%0A%20%20%20%20for%20order%2C%20color%20in%20colors.items()%3A%0A%20%20%20%20%20%20%20%20b%2C%20a%20%3D%20signal.butter(order%2C%2010%2C%20%22low%22%2C%20analog%3DTrue)%0A%20%20%20%20%20%20%20%20w%2C%20h%20%3D%20signal.freqs(b%2C%20a%2C%20worN%3Dnp.linspace(5%2C%2015%2C%20200))%0A%20%20%20%20%20%20%20%20axi.plot(w%20%2F%2010%2C%2020%20*%20np.log10(abs(h))%2C%20color%2C%20linewidth%3D2)%0A%20%20%20%20axi.set_xticks((0.6%2C%201%2C%201.4))%0A%20%20%20%20axi.set_yticks((-6%2C%20-3%2C%200))%0A%20%20%20%20axi.set_ylim(%5B-7%2C%201%5D)%0A%20%20%20%20axi.set_xlim(%5B0.5%2C%201.5%5D)%0A%20%20%20%20axi.grid(which%3D%22both%22%2C%20axis%3D%22both%22)%0A%20%20%20%20fig.suptitle(%0A%20%20%20%20%20%20%20%20%22Bode%20plot%20of%20low-pass%20Butterworth%20filters%20of%20different%20orders%22%2C%20fontsize%3D15%0A%20%20%20%20)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20butterworth_plot()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20All%20orders%20pass%20through%20%24-3%24%20dB%20at%20the%20cutoff%20frequency%20(the%20square).%20Above%20it%2C%20the%20higher%20orders%20fall%20much%20faster%3A%20at%20ten%20times%20the%20cutoff%20a%20first-order%20filter%20has%20attenuated%20by%2020%20dB%2C%20a%20sixth-order%20one%20by%20120%20dB.%20Note%20the%20cost%20in%20the%20right%20panel%2C%20though%3A%20the%20higher%20the%20order%2C%20the%20larger%20the%20phase%20lag.%0A%0A%20%20%20%20**Guiding%20questions%201.**%0A%0A%20%20%20%201.%20Read%20the%20left%20panel%20at%20twice%20the%20cutoff%20frequency.%20How%20many%20decibels%20does%20each%20order%20attenuate%20there%3F%20Does%20it%20match%20%24-6%24%20dB%20per%20octave%20per%20order%3F%0A%20%20%20%202.%20The%20noise%20in%20Pezzack's%20data%20extends%20up%20to%20the%20Nyquist%20frequency%2C%20about%2025%20Hz.%20With%20a%205%20Hz%20cutoff%2C%20roughly%20how%20much%20would%20a%20second-order%20filter%20attenuate%20noise%20at%2020%20Hz%3F%0A%20%20%20%203.%20Why%20not%20always%20use%20a%20very%20high%20order%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Exercise%203.**%20For%20an%20analog%20second-order%20Butterworth%20filter%2C%20estimate%20the%20attenuation%20at%20twice%20and%20four%20times%20its%20cutoff.%20Are%20%24-6%24%20dB%20per%20octave%20per%20order%20exact%20at%20the%20cutoff%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%203%0A%20%20%20%20The%20magnitude%20is%20%241%2F%5Csqrt%7B1%2B(f%2Ff_c)%5E4%7D%24.%20At%20ratios%202%20and%204%20the%20gains%20are%20about%20%24-12.30%24%20and%20%24-24.10%24%20dB.%20The%20%24-12%24%20dB%2Foctave%20slope%20is%20a%20high-frequency%20asymptote%2C%20not%20an%20exact%20rule%20at%20every%20frequency.%20Higher%20order%20attenuates%20more%20strongly%20but%20changes%20phase%20and%20can%20increase%20ringing%3B%20digital%20frequencies%20also%20undergo%20warping.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Butterworth%20filter%20with%20zero-phase%20shift%0A%0A%20%20%20%20The%20phase%20introduced%20by%20the%20Butterworth%20filter%20can%20be%20cancelled%20in%20a%20digital%20implementation%20by%20filtering%20the%20data%20twice%2C%20once%20forward%20and%20once%20backwards.%20The%20lag%20introduced%20in%20the%20first%20pass%20is%20undone%20by%20the%20same%20lag%20in%20the%20opposite%20direction%20in%20the%20second%20pass.%20The%20result%20is%20a%20zero-phase%20shift%20(or%20zero-phase%20lag)%20filter.%20Note%20that%20it%20needs%20the%20whole%20recording%20to%20run%20backwards%20from%20its%20end%2C%20so%20it%20cannot%20be%20used%20in%20real%20time.%0A%0A%20%20%20%20There%20is%20a%20catch.%20Each%20pass%20attenuates%20the%20power%20at%20the%20cutoff%20frequency%20by%20half%2C%20so%20two%20passes%20attenuate%20it%20by%20a%20quarter%2C%20and%20the%20cutoff%20frequency%20of%20the%20combined%20filter%20is%20no%20longer%20the%20one%20we%20asked%20for.%20We%20have%20to%20correct%20the%20cutoff%20given%20to%20%60butter%60%20so%20that%2C%20after%20both%20passes%2C%20the%20attenuation%20is%20again%20by%20half.%20For%20a%20second-order%20Butterworth%20filter%20applied%20%24n%24%20times%2C%20the%20correction%20factor%20is%20(Winter%2C%202009)%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20C%20%3D%20%5Csqrt%5B4%5D%7B2%5E%7B%5Cfrac%7B1%7D%7Bn%7D%7D%20-%201%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20For%20two%20passes%2C%20%24n%3D2%24%2C%20%24C%3D%5Csqrt%5B4%5D%7B2%5E%7B%5Cfrac%7B1%7D%7B2%7D%7D%20-%201%7D%20%5Capprox%200.802%24%2C%20and%20the%20low-frequency%20approximation%20to%20the%20design%20cutoff%20is%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20f_%7Bc%2C%5C%2Cdesign%7D%20%5Capprox%20%5Cfrac%7Bf_%7Bc%2C%5C%2Cdesired%7D%7D%7BC%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20For%20instance%2C%20this%20approximation%20suggests%2012.47%20Hz%20for%20a%20desired%2010%20Hz%20cutoff.%20A%20digital%20filter%20uses%20a%20bilinear%20transformation%2C%20so%20the%20exact%20correction%20must%20be%20made%20in%20the%20prewarped%20frequency%3A%0A%0A%20%20%20%20%24%24f_%7Bc%2C%5C%2Cdesign%7D%3D%5Cfrac%7Bf_s%7D%7B%5Cpi%7D%5Carctan%5C!%5Cleft(%5Cfrac%7B%5Ctan(%5Cpi%20f_%7Bc%2C%5C%2Cdesired%7D%2Ff_s)%7D%7BC%7D%5Cright).%24%24%0A%0A%20%20%20%20At%20%24f_s%3D100%24%20Hz%2C%20this%20gives%20about%2012.25%20Hz.%20The%20examples%20below%20use%20this%20exact%20correction.%20The%20effective%20cutoff%20is%20the%20frequency%20where%20the%20combined%20two-pass%20power%20gain%20is%20%241%2F2%24%2C%20or%20approximately%20%24-3.01%24%20dB.%20The%20combined%20response%20has%20fourth-order%20roll-off%20but%20is%20not%20an%20ordinary%20single-pass%20fourth-order%20Butterworth%20response.%0A%0A%20%20%20%20SciPy's%20function%20%60filtfilt%60%20does%20the%20forward%20and%20backward%20filtering.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20where%20the%20filtfilt%20curve%20will%20be%20relative%20to%20the%20raw%20data%20and%20to%20the%20single-pass%20filter.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(bw_t%2C%20bw_y%2C%20bw_y2%2C%20fs%2C%20np%2C%20plt%2C%20signal)%3A%0A%20%20%20%20%23%20correct%20the%20cutoff%20frequency%20for%20the%20two%20passes%20of%20the%20filter%0A%20%20%20%20C%20%3D%20(2**0.5%20-%201)**0.25%0A%20%20%20%20_b%2C%20_a%20%3D%20signal.butter(2%2C%20(2%20%2F%20np.pi)%20*%20np.arctan(np.tan(np.pi%20*%205%20%2F%20fs)%20%2F%20C)%2C%20btype%3D%22low%22)%0A%20%20%20%20_y3%20%3D%20signal.filtfilt(_b%2C%20_a%2C%20bw_y)%20%20%23%20zero-phase%20filtering%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(9%2C%204))%0A%20%20%20%20_ax.plot(bw_t%2C%20bw_y%2C%20%22r.-%22%2C%20linewidth%3D2%2C%20label%3D%22raw%20data%22)%0A%20%20%20%20_ax.plot(bw_t%2C%20bw_y2%2C%20%22b.-%22%2C%20linewidth%3D2%2C%20label%3D%22lfilter%20%40%205%20Hz%22)%0A%20%20%20%20_ax.plot(bw_t%2C%20_y3%2C%20%22g.-%22%2C%20linewidth%3D2%2C%20label%3D%22filtfilt%20%40%205%20Hz%22)%0A%20%20%20%20_ax.legend(frameon%3DFalse)%0A%20%20%20%20_ax.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22Amplitude%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20(C%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20filtfilt%20output%20sits%20right%20on%20the%20slow%20part%20of%20the%20raw%20data%2C%20with%20no%20delay%2C%20with%20much%20smaller%20edge%20transients%20than%20this%20zero-initialized%20single%20pass.%20Padding%20reduces%20boundary%20artifacts%3B%20it%20does%20not%20guarantee%20their%20absence.%20Let's%20also%20check%20the%20correction%3A%20at%205%20Hz%2C%20one%20pass%20of%20the%20corrected%20filter%20should%20attenuate%20by%20about%201.5%20dB%2C%20so%20that%20two%20passes%20attenuate%20by%203%20dB.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(C%2C%20fs%2C%20np%2C%20signal)%3A%0A%20%20%20%20_%2C%20_h%20%3D%20signal.freqz(*signal.butter(2%2C%20(2%20%2F%20np.pi)%20*%20np.arctan(np.tan(np.pi%20*%205%20%2F%20fs)%20%2F%20C))%2C%20worN%3D%5B5%5D%2C%20fs%3Dfs)%0A%20%20%20%20_%2C%20_h0%20%3D%20signal.freqz(*signal.butter(2%2C%205%20%2F%20(fs%20%2F%202))%2C%20worN%3D%5B5%5D%2C%20fs%3Dfs)%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22Gain%20at%205%20Hz%2C%20corrected%20cutoff%3A%20%20%20one%20pass%20%7B20%20*%20np.log10(abs(_h%5B0%5D))%3A.2f%7D%20dB%2C%20%22%0A%20%20%20%20%20%20%20%20f%22two%20passes%20%7B40%20*%20np.log10(abs(_h%5B0%5D))%3A.2f%7D%20dB%22%0A%20%20%20%20)%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22Gain%20at%205%20Hz%2C%20uncorrected%20cutoff%3A%20one%20pass%20%7B20%20*%20np.log10(abs(_h0%5B0%5D))%3A.2f%7D%20dB%2C%20%22%0A%20%20%20%20%20%20%20%20f%22two%20passes%20%7B40%20*%20np.log10(abs(_h0%5B0%5D))%3A.2f%7D%20dB%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Exercise%204.**%20A%20second-order%20Butterworth%20filter%20is%20applied%20forward%20and%20backward%20without%20correction.%20What%20is%20the%20combined%20gain%20at%20the%20single-pass%20cutoff%3F%20How%20can%20you%20check%20a%20corrected%20design%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%204%0A%20%20%20%20Each%20pass%20has%20amplitude%20%241%2F%5Csqrt2%24%2C%20so%20the%20combined%20amplitude%20is%20%241%2F2%24%20and%20the%20power%20is%20%241%2F4%24%3A%20%2420%5Clog_%7B10%7D(1%2F2)%3D-6.0206%24%20dB.%20Evaluate%20%60freqz%60%20at%20the%20desired%20frequency%3B%20for%20two%20passes%20require%20%60abs(h)**4%20%3D%200.5%60%2C%20or%20%6040*np.log10(abs(h))%20%3D%20-3.0103%60.%20The%20prewarped%20correction%20above%20satisfies%20this%20condition.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Without%20the%20correction%2C%20filtering%20twice%20at%20a%20nominal%205%20Hz%20attenuates%205%20Hz%20by%206%20dB%3A%20the%20real%20cutoff%20would%20be%20lower%20than%20the%20one%20you%20report%20in%20your%20methods%20section.%0A%0A%20%20%20%20%23%23%23%20Critically%20damped%20digital%20filter%0A%0A%20%20%20%20A%20problem%20with%20the%20low-pass%20Butterworth%20filter%20is%20that%20it%20tends%20to%20overshoot%20or%20undershoot%20data%20with%20rapid%20changes%20(see%2C%20for%20example%2C%20Winter%20(2009)%2C%20Robertson%20et%20al.%20(2013)%2C%20and%20Robertson%20and%20Dowling%20(2003)).%20The%20Butterworth%20filter%20behaves%20as%20an%20underdamped%20second-order%20system%3B%20the%20analog%20critically%20damped%20response%20avoids%20that%20overshoot.%20Its%20digital%20approximation%20must%20still%20be%20checked%2C%20especially%20close%20to%20Nyquist.%0A%0A%20%20%20%20The%20function%20%60critic_damp%60%20below%20calculates%20the%20coefficients%20(the%20%24b%24's%20and%20%24a%24's)%20of%20an%20IIR%20critically%20damped%20digital%20filter%2C%20and%20corrects%20the%20cutoff%20frequency%20for%20the%20number%20of%20passes%20of%20the%20filter.%20The%20calculation%20is%20very%20similar%20to%20that%20of%20the%20Butterworth%20filter%2C%20so%20the%20function%20can%20also%20calculate%20the%20Butterworth%20coefficients%20if%20asked%20to%3B%20only%20the%20damping%20term%20and%20the%20correction%20factor%20differ.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.function%0Adef%20critic_damp(fcut%2C%20freq%2C%20npass%3D2%2C%20fcorr%3DTrue%2C%20filt%3D%22critic%22)%3A%0A%20%20%20%20%22%22%22Coefficients%20of%20a%20critically%20damped%20or%20Butterworth%20digital%20low-pass%20filter.%0A%0A%20%20%20%20A%20problem%20with%20a%20low-pass%20Butterworth%20filter%20is%20that%20it%20tends%20to%20overshoot%0A%20%20%20%20or%20undershoot%20data%20with%20rapid%20changes%20(see%20for%20example%2C%20Winter%20(2009)%2C%0A%20%20%20%20Robertson%20et%20al.%20(2013)%2C%20and%20Robertson%20%26%20Dowling%20(2003)).%0A%20%20%20%20The%20analog%20critically%20damped%20response%20avoids%20step%20overshoot.%20Check%20the%0A%20%20%20%20actual%20digital%20response%2C%20particularly%20for%20cutoffs%20near%20Nyquist.%0A%0A%20%20%20%20Parameters%0A%20%20%20%20----------%0A%20%20%20%20fcut%20%3A%20number%0A%20%20%20%20%20%20%20%20desired%20cutoff%20frequency%20for%20the%20low-pass%20digital%20filter%20(Hz).%0A%20%20%20%20freq%20%3A%20number%0A%20%20%20%20%20%20%20%20sampling%20frequency%20(Hz).%0A%20%20%20%20npass%20%3A%20number%2C%20optional%20(default%20%3D%202)%0A%20%20%20%20%20%20%20%20number%20of%20passes%20the%20filter%20will%20be%20applied.%0A%20%20%20%20%20%20%20%20choose%202%20for%20two%20passes%20of%20the%20second-order%20section%20(fourth-order%20roll-off).%0A%20%20%20%20fcorr%20%3A%20bool%2C%20optional%20(default%20%3D%20True)%0A%20%20%20%20%20%20%20%20correct%20(True)%20or%20not%20the%20cutoff%20frequency%20for%20the%20number%20of%20passes.%0A%20%20%20%20filt%20%3A%20string%20('critic'%2C%20'butter')%2C%20optional%20(default%20%3D%20'critic')%0A%20%20%20%20%20%20%20%20'critic'%20to%20calculate%20coefficients%20for%20a%20critically%20damped%20filter%2C%0A%20%20%20%20%20%20%20%20'butter'%20to%20calculate%20coefficients%20for%20a%20Butterworth%20filter.%0A%0A%20%20%20%20Returns%0A%20%20%20%20-------%0A%20%20%20%20b%20%3A%201D%20array%0A%20%20%20%20%20%20%20%20b%20coefficients%20for%20the%20filter%0A%20%20%20%20a%20%3A%201D%20array%0A%20%20%20%20%20%20%20%20a%20coefficients%20for%20the%20filter%0A%20%20%20%20fc%20%3A%20number%0A%20%20%20%20%20%20%20%20digital%20design-frequency%20parameter%20after%20exact%20prewarped%20correction%0A%0A%20%20%20%20References%0A%20%20%20%20----------%0A%20%20%20%20Robertson%20DG%2C%20Dowling%20JJ%20(2003)%20Design%20and%20responses%20of%20Butterworth%20and%0A%20%20%20%20critically%20damped%20digital%20filters.%20J%20Electromyogr%20Kinesiol%2013%2C%20569-573.%0A%20%20%20%20%22%22%22%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20filt%20%3D%20filt.lower()%0A%20%20%20%20if%20filt%20not%20in%20(%22critic%22%2C%20%22butter%22)%3A%0A%20%20%20%20%20%20%20%20raise%20ValueError(f%22Invalid%20option%20for%20parameter%20filt%3A%20%7Bfilt!r%7D%22)%0A%20%20%20%20if%20not%20np.isfinite(freq)%20or%20freq%20%3C%3D%200%20or%20not%200%20%3C%20fcut%20%3C%20freq%20%2F%202%3A%0A%20%20%20%20%20%20%20%20raise%20ValueError(%22Require%20a%20positive%20sampling%20rate%20and%200%20%3C%20fcut%20%3C%20freq%2F2.%22)%0A%20%20%20%20if%20not%20np.isfinite(npass)%20or%20npass%20%3C%201%20or%20int(npass)%20!%3D%20npass%3A%0A%20%20%20%20%20%20%20%20raise%20ValueError(%22npass%20must%20be%20a%20positive%20integer.%22)%0A%0A%20%20%20%20%23%20Correct%20the%20prewarped%20frequency%2C%20not%20its%20digital%20frequency%20in%20hertz.%0A%20%20%20%20corr%20%3D%201.0%0A%20%20%20%20if%20fcorr%3A%0A%20%20%20%20%20%20%20%20if%20filt%20%3D%3D%20%22critic%22%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20corr%20%3D%201%20%2F%20np.sqrt(2%20**%20(1%20%2F%20(2%20*%20npass))%20-%201)%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20corr%20%3D%201%20%2F%20(2%20**%20(1%20%2F%20npass)%20-%201)%20**%200.25%0A%20%20%20%20wc%20%3D%20np.tan(np.pi%20*%20fcut%20%2F%20freq)%20*%20corr%0A%20%20%20%20fc%20%3D%20freq%20%2F%20np.pi%20*%20np.arctan(wc)%0A%20%20%20%20%23%20low-pass%20coefficients%0A%20%20%20%20k1%20%3D%20np.sqrt(2)%20*%20wc%20if%20filt%20%3D%3D%20%22butter%22%20else%202%20*%20wc%0A%20%20%20%20k2%20%3D%20wc%20*%20wc%0A%20%20%20%20a0%20%3D%20k2%20%2F%20(1%20%2B%20k1%20%2B%20k2)%0A%20%20%20%20a1%20%3D%202%20*%20a0%0A%20%20%20%20a2%20%3D%20k2%20%2F%20(1%20%2B%20k1%20%2B%20k2)%0A%20%20%20%20b1%20%3D%202%20*%20a0%20*%20(1%20%2F%20k2%20-%201)%0A%20%20%20%20b2%20%3D%201%20-%20(a0%20%2B%20a1%20%2B%20a2%20%2B%20b1)%0A%20%20%20%20%23%20transform%20parameters%20to%20be%20consistent%20with%20SciPy%20butter%20output%0A%20%20%20%20b%20%3D%20np.array(%5Ba0%2C%20a1%2C%20a2%5D)%0A%20%20%20%20a%20%3D%20np.array(%5B1%2C%20-b1%2C%20-b2%5D)%0A%0A%20%20%20%20return%20b%2C%20a%2C%20fc%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20First%2C%20a%20sanity%20check%3A%20with%20%60filt%3D'butter'%60%2C%20the%20function%20should%20give%20the%20same%20coefficients%20as%20SciPy's%20%60butter%60%20at%20the%20corrected%20cutoff%20frequency.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20signal)%3A%0A%20%20%20%20_b_bw%2C%20_a_bw%2C%20_fc_bw%20%3D%20critic_damp(fcut%3D10%2C%20freq%3D100%2C%20npass%3D2%2C%20filt%3D%22butter%22)%0A%20%20%20%20_b_sp%2C%20_a_sp%20%3D%20signal.butter(2%2C%20_fc_bw%20%2F%20(100%20%2F%202))%0A%20%20%20%20print(f%22Butterworth%2C%20corrected%20cutoff%3A%20%7B_fc_bw%3A.2f%7D%20Hz%22)%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20%22Same%20coefficients%20as%20scipy.signal.butter%3A%22%2C%0A%20%20%20%20%20%20%20%20np.allclose(_b_bw%2C%20_b_sp)%20and%20np.allclose(_a_bw%2C%20_a_sp)%2C%0A%20%20%20%20)%0A%0A%20%20%20%20_b_cd%2C%20_a_cd%2C%20_fc_cd%20%3D%20critic_damp(fcut%3D10%2C%20freq%3D100%2C%20npass%3D2%2C%20filt%3D%22critic%22)%0A%20%20%20%20print(f%22Critically%20damped%2C%20corrected%20cutoff%3A%20%7B_fc_cd%3A.2f%7D%20Hz%22)%0A%20%20%20%20print(%22%20%20b%3A%22%2C%20_b_cd%2C%20%22%5Cn%20%20a%3A%22%2C%20_a_cd)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20two%20filters%20need%20very%20different%20corrections%3A%20the%20critically%20damped%20response%20needs%20a%20larger%20design-frequency%20parameter%20because%20it%20falls%20off%20more%20gently.%20The%20factor%20of%20about%202.30%20applies%20in%20the%20low-frequency%20limit%3B%20the%20exact%20digital%20value%20depends%20on%20the%20sampling%20rate.%0A%0A%20%20%20%20Now%20the%20reason%20for%20its%20existence.%20Let's%20filter%20a%20step%2C%20the%20most%20rapid%20change%20possible%2C%20sampled%20at%20100%20Hz%2C%20with%20both%20filters%20at%20a%2010%20Hz%20cutoff%2C%20two%20passes%20each.%20We%20add%20a%20third%2C%20nonrecursive%20option%20for%20comparison%3A%20an%20FIR%20low-pass%20filter%20with%2011%20coefficients%2C%20designed%20with%20SciPy's%20%60firwin%60%20(and%2C%20unlike%20the%20other%20two%2C%20without%20correcting%20its%20cutoff%20for%20the%20two%20passes).%0A%0A%20%20%20%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20which%20filters%20will%20overshoot%20the%20step%2C%20and%20by%20how%20much.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20signal)%3A%0A%20%20%20%20_y%20%3D%20np.hstack((np.zeros(20)%2C%20np.ones(20)))%0A%20%20%20%20_t%20%3D%20np.linspace(0%2C%200.39%2C%2040)%20-%200.19%0A%0A%20%20%20%20_b_cd%2C%20_a_cd%2C%20_%20%3D%20critic_damp(fcut%3D10%2C%20freq%3D100%2C%20npass%3D2%2C%20filt%3D%22critic%22)%0A%20%20%20%20_b_bw%2C%20_a_bw%2C%20_%20%3D%20critic_damp(fcut%3D10%2C%20freq%3D100%2C%20npass%3D2%2C%20filt%3D%22butter%22)%0A%20%20%20%20_b_fir%20%3D%20signal.firwin(numtaps%3D11%2C%20cutoff%3D10%2C%20fs%3D100)%0A%20%20%20%20_filtered%20%3D%20%7B%0A%20%20%20%20%20%20%20%20%22Butterworth%22%3A%20(%22b%22%2C%20signal.filtfilt(_b_bw%2C%20_a_bw%2C%20_y))%2C%0A%20%20%20%20%20%20%20%20%22Critically%20damped%22%3A%20(%22r%22%2C%20signal.filtfilt(_b_cd%2C%20_a_cd%2C%20_y))%2C%0A%20%20%20%20%20%20%20%20%22FIR%2C%2011%20coefficients%22%3A%20(%22g%22%2C%20signal.filtfilt(_b_fir%2C%201%2C%20_y))%2C%0A%20%20%20%20%7D%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(9%2C%204))%0A%20%20%20%20_ax.plot(_t%2C%20_y%2C%20%22k%22%2C%20linewidth%3D2%2C%20drawstyle%3D%22steps-mid%22%2C%20label%3D%22raw%20data%22)%0A%20%20%20%20for%20_name%2C%20(_color%2C%20_yf)%20in%20_filtered.items()%3A%0A%20%20%20%20%20%20%20%20_ax.plot(_t%2C%20_yf%2C%20_color%20%2B%20%22.-%22%2C%20linewidth%3D2%2C%20label%3D_name)%0A%20%20%20%20%20%20%20%20print(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22%7B_name%3A21s%7D%20overshoot%3A%20%7B100%20*%20(_yf.max()%20-%201)%3A5.1f%7D%25%2C%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22undershoot%3A%20%7Bmax(0.0%2C%20-100%20*%20_yf.min())%3A5.1f%7D%25%22%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20_ax.legend(frameon%3DFalse%2C%20loc%3D%22upper%20left%22)%0A%20%20%20%20_ax.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22Amplitude%22)%0A%20%20%20%20_ax.set_title(%22Step%20response%3A%20100%20Hz%20sampling%2C%2010%20Hz%20cutoff%2C%20zero-phase%20filters%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20Butterworth%20filter%20overshoots%20the%20step%20by%20almost%204%25%2C%20and%20undershoots%20it%20by%20the%20same%20amount%20*before*%20the%20step%20happens%3A%20this%20is%20a%20zero-phase%20filter%2C%20so%20the%20ringing%20appears%20on%20both%20sides.%20The%20critically%20damped%20filter%20does%20neither%2C%20and%20neither%20does%20this%20FIR%20filter.%20For%20a%20step%2C%20an%20overshoot%20of%204%25%20may%20not%20matter%3B%20in%20the%20second%20derivative%20of%20a%20step%2C%20such%20as%20a%20foot%20striking%20the%20ground%2C%20the%20ringing%20becomes%20a%20spurious%20oscillation%20in%20the%20acceleration.%0A%0A%20%20%20%20**Challenge%201.**%20Change%20the%20cutoff%20frequency%20in%20the%20cell%20above%20to%205%20Hz%20and%20to%2020%20Hz.%20Does%20the%20overshoot%20of%20the%20Butterworth%20filter%20change%3F%20What%20about%20the%20width%20of%20the%20transition%3F%20What%20property%20of%20the%20filter%20determines%20the%20overshoot%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Exercise%201.**%20Challenge%200%20needs%20more%20than%20an%20error%20amplitude.%20If%20independent%20position%20errors%20have%20standard%20deviation%201%20mm%20at%20100%20Hz%2C%20what%20is%20the%20standard%20deviation%20of%20the%20central%20second-difference%20acceleration%20error%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%201%0A%20%20%20%20The%20error%20is%20%24(e_%7Bi%2B1%7D-2e_i%2Be_%7Bi-1%7D)%2F%5CDelta%20t%5E2%24.%20Independence%20gives%20variance%20%24(1%2B4%2B1)%5Csigma%5E2%2F%5CDelta%20t%5E4%24%2C%20so%20the%20standard%20deviation%20is%20%24%5Csqrt6(0.001)%2F(0.01)%5E2%3D24.5%24%20m%2Fs%C2%B2.%20Correlated%20errors%20give%20a%20different%20value.%20A%20position-error%20amplitude%20alone%20cannot%20determine%20the%20acceleration%20error.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Exercise%205.**%20Answer%20Challenge%201%3A%20how%20does%20changing%20cutoff%20alter%20a%20step%20response%2C%20and%20why%20can%20a%20zero-phase%20filter%20show%20ringing%20before%20the%20step%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%205%0A%20%20%20%20A%20lower%20cutoff%20spreads%20the%20transition%20over%20more%20samples%3B%20a%20higher%20cutoff%20narrows%20it.%20Damping%20and%20pole%20locations%20govern%20overshoot%2C%20not%20zero%20phase%20alone.%20Forward%E2%80%93backward%20filtering%20is%20noncausal%2C%20so%20its%20symmetric%20impulse%20response%20can%20spread%20ringing%20to%20both%20sides%20of%20the%20event.%20At%20cutoffs%20approaching%20Nyquist%2C%20digital%20pole%20locations%20can%20also%20change%20the%20overshoot%3B%20inspect%20the%20actual%20response.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Moving%20filters%0A%0A%20%20%20%20%23%23%23%20Moving-average%20filter%2C%20four%20ways%0A%0A%20%20%20%20Back%20to%20the%20moving%20average.%20Here%20are%20four%20implementations%20of%20it%3A%20the%20naive%20loop%20from%20the%20beginning%20of%20the%20notebook%2C%20one%20with%20a%20cumulative%20sum%2C%20one%20with%20a%20convolution%2C%20and%20one%20with%20the%20general%20filter%20%60lfilter%60%2C%20using%20the%20coefficients%20%24b_k%20%3D%201%2Fm%24%20and%20no%20%24a%24%20coefficients.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.function%0Adef%20moving_average_cumsum(x%2C%20window)%3A%0A%20%20%20%20%22%22%22Moving%20average%20of%20%60x%60%20with%20window%20size%20%60window%60%2C%20by%20a%20cumulative%20sum.%22%22%22%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20xsum%20%3D%20np.cumsum(x)%0A%20%20%20%20xsum%5Bwindow%3A%5D%20%3D%20xsum%5Bwindow%3A%5D%20-%20xsum%5B%3A-window%5D%0A%20%20%20%20return%20xsum%5Bwindow%20-%201%20%3A%5D%20%2F%20window%0A%0A%0A%40app.function%0Adef%20moving_average_convolve(x%2C%20window)%3A%0A%20%20%20%20%22%22%22Moving%20average%20of%20%60x%60%20with%20window%20size%20%60window%60%2C%20by%20convolution.%22%22%22%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20return%20np.convolve(x%2C%20np.ones(window)%20%2F%20window%2C%20%22same%22)%0A%0A%0A%40app.function%0Adef%20moving_average_lfilter(x%2C%20window)%3A%0A%20%20%20%20%22%22%22Moving%20average%20of%20%60x%60%20with%20window%20size%20%60window%60%2C%20as%20an%20FIR%20filter.%22%22%22%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20from%20scipy.signal%20import%20lfilter%0A%0A%20%20%20%20return%20lfilter(np.ones(window)%20%2F%20window%2C%201%2C%20x)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Let's%20test%20these%20versions%20on%20the%20same%20kind%20of%20step%20signal%20as%20before.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20which%20of%20the%20four%20curves%20will%20coincide%2C%20and%20which%20will%20be%20shifted%20in%20time%20relative%20to%20the%20raw%20data.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt)%3A%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(seed%3D1)%0A%20%20%20%20ma_x%20%3D%20_rng.standard_normal(300)%20%2F%2010%0A%20%20%20%20ma_x%5B100%3A200%5D%20%2B%3D%201%0A%20%20%20%20ma_window%20%3D%2010%0A%0A%20%20%20%20_versions%20%3D%20%7B%0A%20%20%20%20%20%20%20%20%22loop%22%3A%20(moving_average%2C%20%22y-%22)%2C%0A%20%20%20%20%20%20%20%20%22cumsum%22%3A%20(moving_average_cumsum%2C%20%22m--%22)%2C%0A%20%20%20%20%20%20%20%20%22convolve%22%3A%20(moving_average_convolve%2C%20%22r-%22)%2C%0A%20%20%20%20%20%20%20%20%22lfilter%22%3A%20(moving_average_lfilter%2C%20%22g-%22)%2C%0A%20%20%20%20%7D%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(10%2C%205))%0A%20%20%20%20_ax.plot(ma_x%2C%20%22b-%22%2C%20linewidth%3D1%2C%20label%3D%22raw%20data%22)%0A%20%20%20%20for%20_name%2C%20(_func%2C%20_style)%20in%20_versions.items()%3A%0A%20%20%20%20%20%20%20%20_y%20%3D%20_func(ma_x%2C%20ma_window)%0A%20%20%20%20%20%20%20%20_ax.plot(_y%2C%20_style%2C%20linewidth%3D2%2C%20label%3Df%22moving%20average%2C%20%7B_name%7D%22)%0A%20%20%20%20%20%20%20%20print(f%22%7B_name%3A9s%7D%20output%20length%3A%20%7B_y.size%7D%22)%0A%20%20%20%20_ax.legend(frameon%3DFalse%2C%20loc%3D%22upper%20right%22)%0A%20%20%20%20_ax.set_xlabel(%22Sample%22)%0A%20%20%20%20_ax.set_ylabel(%22Amplitude%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20ma_window%2C%20ma_x%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20loop%20and%20the%20cumulative-sum%20versions%20give%20identical%20results%2C%20%24m-1%24%20samples%20shorter%20than%20the%20input.%20The%20convolution%20and%20%60lfilter%60%20versions%20keep%20the%20length%20of%20the%20input.%20But%20only%20the%20convolution%20version%20is%20aligned%20with%20the%20raw%20data.%20The%20other%20three%20compute%20the%20same%20averages%20and%20attribute%20them%20to%20different%20instants%3A%0A%0A%20%20%20%20-%20the%20loop%20and%20cumsum%20versions%20put%20the%20mean%20of%20%60x%5Bi%5D%60%20to%20%60x%5Bi%2Bm-1%5D%60%20at%20sample%20%60i%60%2C%20the%20*start*%20of%20the%20window%2C%20so%20they%20lead%20the%20data%20by%20%24(m-1)%2F2%24%20samples%3B%0A%20%20%20%20-%20%60lfilter%60%20puts%20the%20mean%20of%20%60x%5Bi-m%2B1%5D%60%20to%20%60x%5Bi%5D%60%20at%20sample%20%60i%60%2C%20the%20*end*%20of%20the%20window%2C%20so%20it%20lags%20by%20%24(m-1)%2F2%24%20samples%2C%20as%20a%20causal%20filter%20must%3B%0A%20%20%20%20-%20%60np.convolve%60%20with%20%60'same'%60%20puts%20it%20at%20the%20*centre*%20of%20the%20window.%0A%0A%20%20%20%20It%20is%20the%20same%20question%20as%20where%20to%20place%20a%20finite-difference%20velocity%20in%20%5BKinematics%20of%20a%20particle%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fkinematics_particle.py)%3A%20an%20average%20over%20an%20interval%20belongs%20to%20its%20middle.%20The%20first%20three%20could%20be%20fixed%20by%20shifting%20their%20output%2C%20or%2C%20for%20%60lfilter%60%2C%20by%20using%20%60filtfilt%60%20instead.%0A%0A%20%20%20%20Now%20the%20speed.%20**Before%20you%20run%20the%20next%20cell**%2C%20guess%20how%20much%20slower%20the%20loop%20is%20than%20the%20others.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(ma_window%2C%20ma_x%2C%20timeit)%3A%0A%20%20%20%20for%20_func%20in%20(%0A%20%20%20%20%20%20%20%20moving_average%2C%0A%20%20%20%20%20%20%20%20moving_average_cumsum%2C%0A%20%20%20%20%20%20%20%20moving_average_convolve%2C%0A%20%20%20%20%20%20%20%20moving_average_lfilter%2C%0A%20%20%20%20)%3A%0A%20%20%20%20%20%20%20%20_time%20%3D%20timeit.timeit(lambda%3A%20_func(ma_x%2C%20ma_window)%2C%20number%3D200)%20%2F%20200%0A%20%20%20%20%20%20%20%20print(f%22%7B_func.__name__%3A24s%7D%20%7B1e6%20*%20_time%3A9.1f%7D%20microseconds%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20exact%20numbers%20depend%20on%20your%20computer%2C%20but%20the%20pattern%20should%20not%3A%20the%20version%20with%20the%20Python%20loop%20is%20around%20a%20hundred%20times%20slower%20than%20the%20others%2C%20which%20hand%20the%20loop%20to%20compiled%20code%20inside%20NumPy%20or%20SciPy.%20Avoid%20explicit%20loops%20over%20the%20samples%20of%20a%20signal%20whenever%20a%20vectorized%20operation%20exists.%0A%0A%20%20%20%20%23%23%23%20Moving-RMS%20filter%0A%0A%20%20%20%20The%20root%20mean%20square%20(RMS)%20is%20a%20measure%20of%20the%20absolute%20amplitude%20of%20the%20data%2C%20useful%20when%20the%20data%20have%20positive%20and%20negative%20values%2C%20as%20an%20EMG%20does.%20The%20RMS%20is%20defined%20as%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20RMS%20%3D%20%5Csqrt%7B%5Cfrac%7B1%7D%7BN%7D%5Csum_%7Bi%3D1%7D%5E%7BN%7D%20x_i%5E2%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%2C%20like%20the%20moving%20average%2C%20the%20moving%20RMS%20applies%20it%20to%20a%20sliding%20window%20of%20%24m%24%20samples%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20y%5Bi%5D%20%3D%20%5Csqrt%7B%5Cfrac%7B1%7D%7Bm%7D%5Csum_%7Bj%3D0%7D%5E%7Bm-1%7D%20(x%5Bi%2Bj%5D)%5E2%7D%20%5Cquad%20%5Ctext%7Bfor%7D%20%5Cquad%20i%3D1%2C%20%5C%3B%20%5Cdots%2C%20%5C%3B%20n-m%2B1%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Here%20are%20two%20implementations%20of%20a%20moving-RMS%20filter%2C%20very%20similar%20to%20the%20moving%20average%3A%20one%20with%20a%20convolution%20over%20a%20centred%20window%20of%20%242m%2B1%24%20samples%2C%20and%20one%20with%20%60filtfilt%60%2C%20which%20applies%20a%20window%20of%20%24m%24%20samples%20forward%20and%20then%20backward.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.function%0Adef%20moving_rms_convolve(x%2C%20window)%3A%0A%20%20%20%20%22%22%22Moving%20RMS%20of%20%60x%60%20over%20a%20centred%20window%20of%202*%60window%60%20%2B%201%20samples.%22%22%22%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20window%20%3D%202%20*%20window%20%2B%201%0A%20%20%20%20return%20np.sqrt(np.convolve(x%20*%20x%2C%20np.ones(window)%20%2F%20window%2C%20%22same%22))%0A%0A%0A%40app.function%0Adef%20moving_rms_filtfilt(x%2C%20window)%3A%0A%20%20%20%20%22%22%22Moving%20RMS%20of%20%60x%60%2C%20a%20%60window%60-sample%20average%20applied%20forward%20and%20backward.%22%22%22%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20from%20scipy.signal%20import%20filtfilt%0A%0A%20%20%20%20return%20np.sqrt(filtfilt(np.ones(window)%20%2F%20window%2C%20%5B1%5D%2C%20x%20*%20x))%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Let's%20apply%20them%20to%20the%20electromyographic%20data%20of%20%5BBasic%20properties%20of%20signals%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fsignal_basic_properties.py)%2C%20sampled%20at%201000%20Hz%2C%20with%20the%20mean%20removed%20first%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20pd%2C%20plt%2C%20repo_file)%3A%0A%20%20%20%20_data%20%3D%20pd.read_csv(%0A%20%20%20%20%20%20%20%20repo_file('data%2Femg.csv')%2C%20header%3DNone%0A%20%20%20%20).to_numpy()%5B300%3A1000%5D%0A%20%20%20%20emg_time%20%3D%20_data%5B%3A%2C%200%5D%0A%20%20%20%20emg%20%3D%20_data%5B%3A%2C%201%5D%20-%20np.mean(_data%5B%3A%2C%201%5D)%0A%20%20%20%20emg_window%20%3D%2050%0A%0A%20%20%20%20_y1%20%3D%20moving_rms_convolve(emg%2C%20emg_window)%0A%20%20%20%20_y2%20%3D%20moving_rms_filtfilt(emg%2C%20emg_window)%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(9%2C%205))%0A%20%20%20%20_ax.plot(emg_time%2C%20emg%2C%20%22k-%22%2C%20linewidth%3D1%2C%20label%3D%22raw%20data%22)%0A%20%20%20%20_ax.plot(emg_time%2C%20_y1%2C%20%22r-%22%2C%20linewidth%3D2%2C%20label%3D%22moving%20RMS%2C%20convolve%22)%0A%20%20%20%20_ax.plot(emg_time%2C%20_y2%2C%20%22b-%22%2C%20linewidth%3D2%2C%20label%3D%22moving%20RMS%2C%20filtfilt%22)%0A%20%20%20%20_ax.legend(frameon%3DFalse%2C%20loc%3D%22upper%20right%22)%0A%20%20%20%20_ax.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22Amplitude%22)%0A%20%20%20%20_ax.set_ylim(-0.1%2C%200.1)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20emg%2C%20emg_window%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Similar%2C%20but%20not%20the%20same.%20Both%20span%20about%20100%20samples%20(0.1%20s)%2C%20but%20they%20weight%20them%20differently.%20The%20convolution%20version%20averages%20101%20samples%20with%20equal%20weights.%20The%20filtfilt%20version%20averages%2050%20samples%20twice%2C%20and%20a%20box%20averaged%20with%20itself%20is%20a%20triangle%2099%20samples%20wide%2C%20which%20weights%20the%20centre%20more%20than%20the%20edges%3B%20its%20result%20is%20a%20little%20smoother.%20Choosing%20a%20window%2C%20and%20its%20shape%2C%20is%20part%20of%20the%20method%2C%20and%20should%20be%20reported.%0A%0A%20%20%20%20The%20convolution%20version%20is%20also%20faster%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(emg%2C%20emg_window%2C%20timeit)%3A%0A%20%20%20%20for%20_func%20in%20(moving_rms_convolve%2C%20moving_rms_filtfilt)%3A%0A%20%20%20%20%20%20%20%20_time%20%3D%20timeit.timeit(lambda%3A%20_func(emg%2C%20emg_window)%2C%20number%3D100)%20%2F%20100%0A%20%20%20%20%20%20%20%20print(f%22%7B_func.__name__%3A20s%7D%20%7B1e6%20*%20_time%3A9.1f%7D%20microseconds%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Exercise%206.**%20Why%20do%20the%20two%20moving-RMS%20functions%20give%20different%20results%20even%20with%20the%20same%20%60window%60%20argument%3F%20Can%20RMS%20be%20computed%20after%20a%20moving%20average%20instead%20of%20before%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%206%0A%20%20%20%20The%20convolution%20version%20averages%20squared%20samples%20over%20a%20rectangle%20of%20length%20%242m%2B1%24.%20The%20forward%E2%80%93backward%20version%20applies%20an%20%24m%24-sample%20rectangle%20twice%2C%20giving%20a%20triangular%20interior%20weighting%20of%20length%20%242m-1%24%2C%20with%20different%20edge%20handling.%20These%20are%20different%20estimators.%20Averaging%20a%20signed%20EMG%20first%20cancels%20positive%20and%20negative%20values%3B%20square%20first%2C%20average%2C%20then%20take%20the%20square%20root%20to%20measure%20amplitude.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Moving-median%20filter%0A%0A%20%20%20%20The%20moving-median%20filter%20is%20similar%20in%20concept%20to%20the%20other%20moving%20filters%20but%20uses%20the%20median%20instead%20of%20the%20mean.%20Because%20the%20median%20ignores%20the%20extreme%20values%20in%20the%20window%2C%20it%20preserves%20abrupt%20changes%20better%20than%20the%20moving%20average%2C%20and%20it%20removes%20isolated%20spikes%20almost%20completely.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20how%20each%20filter%20will%20treat%20the%20two%20jumps.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20signal)%3A%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(seed%3D2)%0A%20%20%20%20_x%20%3D%20_rng.standard_normal(300)%20%2F%2010%0A%20%20%20%20_x%5B100%3A200%5D%20%2B%3D%201%0A%20%20%20%20_window%20%3D%2011%0A%0A%20%20%20%20_y%20%3D%20np.convolve(_x%2C%20np.ones(_window)%20%2F%20_window%2C%20%22same%22)%0A%20%20%20%20_y2%20%3D%20signal.medfilt(_x%2C%20_window)%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(10%2C%204))%0A%20%20%20%20_ax.plot(_x%2C%20%22b-%22%2C%20linewidth%3D1%2C%20label%3D%22raw%20data%22)%0A%20%20%20%20_ax.plot(_y%2C%20%22r-%22%2C%20linewidth%3D2%2C%20label%3D%22moving%20average%22)%0A%20%20%20%20_ax.plot(_y2%2C%20%22g-%22%2C%20linewidth%3D2%2C%20label%3D%22moving%20median%22)%0A%20%20%20%20_ax.legend(frameon%3DFalse%2C%20loc%3D%22upper%20right%22)%0A%20%20%20%20_ax.set_xlabel(%22Sample%22)%0A%20%20%20%20_ax.set_ylabel(%22Amplitude%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20moving%20median%20keeps%20the%20jumps%20nearly%20vertical%2C%20where%20the%20moving%20average%20turns%20them%20into%20ramps.%20Both%20filters%20also%20show%20what%20happens%20at%20the%20ends%20of%20the%20data%3A%20%60np.convolve%60%20with%20%60'same'%60%20and%20%60medfilt%60%20pad%20the%20data%20with%20zeros%20beyond%20its%20edges%2C%20which%20here%2C%20where%20the%20signal%20is%20near%20zero%2C%20does%20little%20harm.%20Every%20filter%20has%20to%20make%20some%20assumption%20about%20the%20data%20it%20does%20not%20have%2C%20and%20the%20edges%20are%20where%20that%20shows.%0A%0A%20%20%20%20%23%23%23%20More%20moving%20filters%0A%0A%20%20%20%20The%20library%20%5Bpandas%5D(https%3A%2F%2Fpandas.pydata.org%2F)%20has%20many%20types%20of%20%5Bmoving%20(rolling)%20window%20functions%5D(https%3A%2F%2Fpandas.pydata.org%2Fdocs%2Fuser_guide%2Fwindow.html)%3A%20mean%2C%20median%2C%20standard%20deviation%2C%20quantiles%2C%20weighted%20windows%20and%20others%2C%20as%20in%20%60pd.Series(x).rolling(11%2C%20center%3DTrue).mean()%60.%20They%20mark%20the%20edges%2C%20where%20the%20window%20is%20incomplete%2C%20with%20%60NaN%60%20instead%20of%20padding%20the%20data.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Filtering%20before%20differentiating%0A%0A%20%20%20%20Back%20to%20the%20function%20%24f%20%3D%20%5Csin(%5Comega%20t)%20%2B%200.1%5Csin(10%5Comega%20t)%24.%20Let's%20filter%20it%20with%20a%20zero-phase%20Butterworth%20low-pass%20at%205%20Hz%2C%20differentiate%20the%20raw%20and%20the%20filtered%20data%20twice%2C%20and%20look%20at%20the%20results%20in%20the%20time%20domain%20and%20in%20the%20frequency%20domain%2C%20using%20the%20%5BFourier%20transform%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FFourier_transform)%20to%20see%20their%20frequency%20content.%0A%0A%20%20%20%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20what%20the%20amplitude%20spectrum%20of%20the%20raw%20%24f''%24%20will%20look%20like%3A%20which%20peak%20will%20be%20larger%2C%20the%20one%20at%201%20Hz%20or%20the%20one%20at%2010%20Hz%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(C%2C%20bw_t%2C%20bw_y%2C%20fs%2C%20np%2C%20plt%2C%20signal)%3A%0A%20%20%20%20_b%2C%20_a%20%3D%20signal.butter(2%2C%20(2%20%2F%20np.pi)%20*%20np.arctan(np.tan(np.pi%20*%205%20%2F%20fs)%20%2F%20C)%2C%20btype%3D%22low%22)%0A%20%20%20%20_yf%20%3D%20signal.filtfilt(_b%2C%20_a%2C%20bw_y)%0A%20%20%20%20%23%20second%20derivatives%0A%20%20%20%20_ydd%20%3D%20np.diff(bw_y%2C%202)%20*%20fs%20*%20fs%20%20%23%20raw%20data%0A%20%20%20%20_yfdd%20%3D%20np.diff(_yf%2C%202)%20*%20fs%20*%20fs%20%20%23%20filtered%20data%0A%0A%0A%20%20%20%20def%20_amplitude_spectrum(x)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Amplitude%20spectrum%20of%20%60x%60%2C%20up%20to%2025%20Hz.%22%22%22%0A%20%20%20%20%20%20%20%20_freqs%20%3D%20np.fft.rfftfreq(x.size%2C%201%20%2F%20fs)%0A%20%20%20%20%20%20%20%20_amp%20%3D%20np.abs(np.fft.rfft(x))%20%2F%20(x.size%20%2F%202)%0A%20%20%20%20%20%20%20%20return%20_freqs%5B_freqs%20%3C%3D%2025%5D%2C%20_amp%5B_freqs%20%3C%3D%2025%5D%0A%0A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(2%2C%202%2C%20figsize%3D(11%2C%205.5))%0A%20%20%20%20_axs%5B0%2C%200%5D.set_title(%22Time%20domain%22)%0A%20%20%20%20_axs%5B0%2C%200%5D.plot(bw_t%2C%20bw_y%2C%20%22r%22%2C%20linewidth%3D2%2C%20label%3D%22raw%20data%22)%0A%20%20%20%20_axs%5B0%2C%200%5D.plot(bw_t%2C%20_yf%2C%20%22b%22%2C%20linewidth%3D2%2C%20label%3D%22filtered%20%40%205%20Hz%22)%0A%20%20%20%20_axs%5B0%2C%200%5D.set_ylabel(%22f%22)%0A%20%20%20%20_axs%5B0%2C%200%5D.legend(frameon%3DFalse)%0A%20%20%20%20_axs%5B0%2C%201%5D.set_title(%22Frequency%20domain%22)%0A%20%20%20%20_axs%5B0%2C%201%5D.plot(*_amplitude_spectrum(bw_y)%2C%20%22r%22%2C%20linewidth%3D2%2C%20label%3D%22raw%20data%22)%0A%20%20%20%20_axs%5B0%2C%201%5D.plot(%0A%20%20%20%20%20%20%20%20*_amplitude_spectrum(_yf)%2C%20%22b--%22%2C%20linewidth%3D2%2C%20label%3D%22filtered%20%40%205%20Hz%22%0A%20%20%20%20)%0A%20%20%20%20_axs%5B0%2C%201%5D.set_ylabel(%22Amplitude%20of%20f%22)%0A%20%20%20%20_axs%5B0%2C%201%5D.legend(frameon%3DFalse)%0A%20%20%20%20_axs%5B1%2C%200%5D.plot(bw_t%5B1%3A-1%5D%2C%20_ydd%2C%20%22r%22%2C%20linewidth%3D2)%0A%20%20%20%20_axs%5B1%2C%200%5D.plot(bw_t%5B1%3A-1%5D%2C%20_yfdd%2C%20%22b%22%2C%20linewidth%3D2)%0A%20%20%20%20_axs%5B1%2C%200%5D.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_axs%5B1%2C%200%5D.set_ylabel(%22f%20''%22)%0A%20%20%20%20_axs%5B1%2C%201%5D.plot(*_amplitude_spectrum(_ydd)%2C%20%22r%22%2C%20linewidth%3D2)%0A%20%20%20%20_axs%5B1%2C%201%5D.plot(*_amplitude_spectrum(_yfdd)%2C%20%22b--%22%2C%20linewidth%3D2)%0A%20%20%20%20_axs%5B1%2C%201%5D.set_xlabel(%22Frequency%20%5BHz%5D%22)%0A%20%20%20%20_axs%5B1%2C%201%5D.set_ylabel(%22Amplitude%20of%20f%20''%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Exercise%207.**%20Compare%20lowering%20the%20cutoff%20to%203%20Hz%20with%20applying%20the%205%20Hz%20forward%E2%80%93backward%20filter%20twice.%20Which%20retains%20more%20of%20a%201%20Hz%20component%20while%20suppressing%20the%2010%20Hz%20interference%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%207%0A%20%20%20%20For%20the%20100%20Hz%20example%2C%20use%20the%20exact%20response%20to%20compare%3A%20the%20corrected%20two-pass%203%20Hz%20filter%20retains%20about%2099.5%25%20of%20the%201%20Hz%20amplitude%20and%20about%201.7%25%20of%20the%2010%20Hz%20amplitude.%20Applying%20the%20corrected%205%20Hz%20two-pass%20filter%20twice%20retains%20about%2099.9%25%20of%201%20Hz%20and%20about%201.4%25%20of%2010%20Hz.%20The%20latter%20keeps%20the%20slow%20component%20better%20here%2C%20but%20has%20higher%20total%20order%20and%20different%20boundary%2Fringing%20behavior.%20Neither%20choice%20is%20universally%20better.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20In%20the%20raw%20data%20the%2010%20Hz%20peak%20is%20a%20tenth%20of%20the%201%20Hz%20peak%3B%20in%20the%20raw%20second%20derivative%20it%20is%20about%20nine%20times%20*larger*%2C%20close%20to%20the%20factor%20of%2010%20in%20amplitude%20computed%20by%20hand%20(a%20little%20less%2C%20because%20a%20finite%20difference%20underestimates%20the%20derivative%20of%20a%20fast%20oscillation).%0A%0A%20%20%20%20After%20filtering%2C%20the%2010%20Hz%20peak%20of%20%24f%24%20has%20almost%20disappeared%3A%20it%20is%20about%201%25%20of%20the%20signal.%20And%20yet%20in%20the%20second%20derivative%20of%20the%20filtered%20data%20it%20is%20back%2C%20about%20as%20large%20as%20the%20signal%20itself.%20Nothing%20went%20wrong.%20Two%20passes%20of%20a%20second-order%20filter%20at%205%20Hz%20reduce%20the%20amplitude%20at%2010%20Hz%20only%20about%20eight%20times%2C%20and%20the%20second%20derivative%20multiplies%20it%20by%20100%20again.%20**A%20filter%20that%20looks%20perfect%20on%20the%20data%20can%20be%20far%20too%20weak%20for%20their%20derivatives**%3B%20how%20much%20attenuation%20you%20need%20depends%20on%20how%20many%20times%20you%20will%20differentiate.%0A%0A%20%20%20%20**Challenge%202.**%20Change%20the%20cutoff%20frequency%20in%20the%20cell%20above%20to%203%20Hz%2C%20then%20go%20back%20to%205%20Hz%20and%20filter%20twice%20with%20%60filtfilt%60%20(four%20passes%20in%20all).%20Which%20works%20better%20for%20%24f''%24%3F%20What%20does%20each%20choice%20cost%20the%201%20Hz%20signal%3F%0A%0A%20%20%20%20%23%23%20Pezzack's%20benchmark%20data%2C%20again%0A%0A%20%20%20%20Now%20we%20can%20attack%20the%20problem%20from%20the%20beginning%20of%20the%20notebook.%20The%20noise%20in%20Pezzack's%20angle%20comes%20from%20small%20random%20errors%20in%20the%20digitization%20of%20each%20frame%2C%20so%20it%20is%20spread%20up%20to%20the%20Nyquist%20frequency%2C%20above%20the%20frequency%20content%20of%20the%20movement.%20Let's%20try%20three%20different%20methods%20to%20attenuate%20it%3A%20a%20%5BButterworth%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FButterworth_filter)%20filter%2C%20a%20%5BSavitzky-Golay%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FSavitzky%25E2%2580%2593Golay_filter)%20filter%2C%20and%20a%20smoothing%20%5Bspline%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FSpline_function).%0A%0A%20%20%20%20The%20Savitzky-Golay%20filter%20and%20the%20spline%20both%20fit%20polynomials%20to%20the%20data%2C%20and%20both%20can%20differentiate%20the%20fitted%20polynomials%20to%20get%20the%20derivatives%20directly%2C%20instead%20of%20differentiating%20the%20data%20numerically.%20Their%20signatures%20in%20SciPy%20are%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20savgol_filter(x%2C%20window_length%2C%20polyorder%2C%20deriv%3D0%2C%20delta%3D1.0%2C%20axis%3D-1%2C%20mode%3D'interp'%2C%20cval%3D0.0)%0A%20%20%20%20splrep(x%2C%20y%2C%20w%3DNone%2C%20xb%3DNone%2C%20xe%3DNone%2C%20k%3D3%2C%20task%3D0%2C%20s%3DNone%2C%20t%3DNone%2C%20full_output%3D0%2C%20per%3D0%2C%20quiet%3D1)%0A%20%20%20%20%60%60%60%0A%0A%20%20%20%20and%20the%20spline%20derivatives%20are%20evaluated%20with%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20splev(x%2C%20tck%2C%20der%3D0%2C%20ext%3D0)%0A%20%20%20%20%60%60%60%0A%0A%20%20%20%20Both%20fits%20behave%20poorly%20at%20the%20ends%20of%20the%20data%2C%20where%20there%20are%20points%20on%20only%20one%20side.%20A%20common%20workaround%20is%20to%20pad%20the%20data%20before%20fitting%2C%20with%20an%20*odd%20extension*%3A%20the%20data%20reflected%20about%20their%20first%20and%20last%20values%2C%20which%20continues%20the%20trend%20of%20the%20data%20instead%20of%20flattening%20it.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.function%0Adef%20odd_extension(x%2C%20n)%3A%0A%20%20%20%20%22%22%22Extend%20%60x%60%20by%20%60n%60%20samples%20at%20each%20end%2C%20reflected%20about%20its%20end%20values.%22%22%22%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20x%20%3D%20np.asarray(x)%0A%20%20%20%20return%20np.concatenate((2%20*%20x%5B0%5D%20-%20x%5Bn%3A0%3A-1%5D%2C%20x%2C%202%20*%20x%5B-1%5D%20-%20x%5B-2%20%3A%20-n%20-%202%20%3A%20-1%5D))%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20We%20will%20compare%20the%20methods%20with%20the%20%5Broot-mean-square%20error%20(RMSE)%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FRoot-mean-square_deviation)%20of%20each%20estimated%20acceleration%20relative%20to%20the%20measured%20reference%2C%20using%20the%20same%20interior%20samples%20for%20every%20method.%20The%20plain%20double%20difference%20had%20an%20RMSE%20of%20about%2011%20rad%2Fs%C2%B2.%20**Before%20you%20run%20the%20next%20cell**%2C%20guess%20how%20much%20each%20method%20will%20reduce%20it%3A%20by%20half%2C%20by%20a%20factor%20of%20ten%2C%20more%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(C%2C%20interpolate%2C%20np%2C%20plt%2C%20pz_aacc%2C%20pz_disp%2C%20pz_dt%2C%20pz_time%2C%20signal)%3A%0A%20%20%20%20%23%20Butterworth%20filter%20at%209%20Hz%2C%20zero%20phase%0A%20%20%20%20_b%2C%20_a%20%3D%20signal.butter(2%2C%20(2%20%2F%20np.pi)%20*%20np.arctan(np.tan(np.pi%20*%209%20*%20pz_dt)%20%2F%20C))%0A%20%20%20%20_aacc_bw%20%3D%20np.diff(signal.filtfilt(_b%2C%20_a%2C%20pz_disp)%2C%202)%20%2F%20pz_dt**2%20%20%23%202%20samples%20shorter%0A%0A%20%20%20%20%23%20pad%20the%20data%20at%20the%20extremities%20to%20reduce%20edge%20effects%0A%20%20%20%20_n%20%3D%2011%0A%20%20%20%20_disp_pad%20%3D%20odd_extension(pz_disp%2C%20_n)%0A%20%20%20%20_time_pad%20%3D%20odd_extension(pz_time%2C%20_n)%0A%0A%20%20%20%20%23%20Savitzky-Golay%20filter%0A%20%20%20%20_aacc_sg%20%3D%20signal.savgol_filter(%0A%20%20%20%20%20%20%20%20_disp_pad%2C%20window_length%3D5%2C%20polyorder%3D3%2C%20deriv%3D2%2C%20delta%3Dpz_dt%0A%20%20%20%20)%5B_n%3A-_n%5D%0A%0A%20%20%20%20%23%20quintic%20smoothing%20spline%0A%20%20%20%20_s%20%3D%200.15%20*%20np.var(_disp_pad)%20%2F%20np.size(_disp_pad)%0A%20%20%20%20_tck%20%3D%20interpolate.splrep(_time_pad%2C%20_disp_pad%2C%20k%3D5%2C%20s%3D_s)%0A%20%20%20%20_aacc_sp%20%3D%20interpolate.splev(_time_pad%2C%20_tck%2C%20der%3D2)%5B_n%3A-_n%5D%0A%0A%20%20%20%20_rmse_bw%20%3D%20np.sqrt(np.mean((_aacc_bw%20-%20pz_aacc%5B1%3A-1%5D)%20**%202))%0A%20%20%20%20_rmse_sg%20%3D%20np.sqrt(np.mean((_aacc_sg%5B1%3A-1%5D%20-%20pz_aacc%5B1%3A-1%5D)%20**%202))%0A%20%20%20%20_rmse_sp%20%3D%20np.sqrt(np.mean((_aacc_sp%5B1%3A-1%5D%20-%20pz_aacc%5B1%3A-1%5D)%20**%202))%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(11%2C%204.5))%0A%20%20%20%20_ax.plot(pz_time%2C%20pz_aacc%2C%20%22g%22%2C%20linewidth%3D3%2C%20label%3D%22Accelerometer%20(reference)%22)%0A%20%20%20%20_ax.plot(pz_time%5B1%3A-1%5D%2C%20_aacc_bw%2C%20%22r%22%2C%20label%3Df%22Butterworth%209%20Hz%3A%20RMSE%20%3D%20%7B_rmse_bw%3A.2f%7D%22)%0A%20%20%20%20_ax.plot(pz_time%2C%20_aacc_sg%2C%20%22b%22%2C%20label%3Df%22Savitzky-Golay%205%20points%3A%20RMSE%20%3D%20%7B_rmse_sg%3A.2f%7D%22)%0A%20%20%20%20_ax.plot(%0A%20%20%20%20%20%20%20%20pz_time%2C%20_aacc_sp%2C%20%22m%22%2C%20label%3Df%22Quintic%20spline%2C%20s%20%3D%20%7B_s%3A.5f%7D%3A%20RMSE%20%3D%20%7B_rmse_sp%3A.2f%7D%22%0A%20%20%20%20)%0A%20%20%20%20_ax.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22Angular%20acceleration%20%5Brad%2Fs%24%5E2%24%5D%22)%0A%20%20%20%20_ax.set_title(%22Pezzack's%20benchmark%20data%22)%0A%20%20%20%20_ax.legend(frameon%3DFalse%2C%20loc%3D%22upper%20left%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20All%20three%20methods%20cut%20the%20error%20by%20more%20than%20half%2C%20to%20about%204%E2%80%935%20rad%2Fs%C2%B2%2C%20and%20all%20three%20follow%20the%20reference%20acceleration%20closely.%20None%20of%20them%20is%20clearly%20better%20than%20the%20others%20here.%20With%20all%20of%20them%2C%20and%20particularly%20with%20the%20spline%2C%20the%20parameters%20had%20to%20be%20tuned%3A%20the%209%20Hz%20cutoff%2C%20the%205-point%20window%2C%20the%20smoothing%20factor%20%60s%60.%20The%20Butterworth%20filter%20is%20often%20the%20easiest%20to%20tune%2C%20because%20a%20cutoff%20frequency%20is%20a%20parameter%20with%20a%20physical%20meaning%20for%20human%20movement.%0A%0A%20%20%20%20But%20where%20did%209%20Hz%20come%20from%3F%20It%20was%20chosen%20by%20hand%2C%20looking%20at%20the%20result%2C%20and%20here%20we%20had%20the%20luxury%20of%20knowing%20the%20reference%20acceleration.%20In%20a%20real%20experiment%20you%20would%20not.%0A%0A%20%20%20%20**Challenge%203.**%20The%20data%20file%20also%20has%20a%20second%2C%20noisier%20version%20of%20the%20angle%2C%20loaded%20above%20as%20%60pz_disp_noisy%60.%0A%0A%20%20%20%201.%20Differentiate%20it%20twice%20without%20filtering%20and%20compare%20the%20RMSE%20with%20that%20of%20%60pz_disp%60.%0A%20%20%20%202.%20Filter%20it%20with%20the%20same%20Butterworth%20filter%20at%209%20Hz.%20Is%209%20Hz%20still%20a%20good%20choice%3F%20Try%20a%20few%20other%20cutoff%20frequencies%20and%20find%20the%20one%20that%20gives%20the%20smallest%20RMSE.%0A%20%20%20%203.%20Without%20the%20accelerometer%2C%20how%20would%20you%20have%20chosen%20the%20cutoff%3F%20Keep%20your%20answer%20for%20%5BResidual%20analysis%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fresidual_analysis.py).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Kinematics%20of%20a%20ball%20toss%0A%0A%20%20%20%20Filtering%20is%20not%20the%20only%20way%20to%20deal%20with%20noise.%20Let's%20analyse%20the%20kinematic%20data%20of%20a%20ball%20tossed%20in%20the%20air.%20They%20were%20obtained%20with%20%5BTracker%5D(https%3A%2F%2Fopensourcephysics.github.io%2Ftracker-website%2F)%2C%20a%20free%20video%20analysis%20and%20modeling%20tool%20built%20on%20the%20%5BOpen%20Source%20Physics%5D(https%3A%2F%2Fwww.compadre.org%2Fosp%2F)%20(OSP)%20Java%20framework%2C%20from%20the%20video%20*balltossout.mov*%20in%20the%20mechanics%20video%20collection%20on%20the%20Tracker%20website.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20pd%2C%20plt%2C%20repo_file)%3A%0A%20%20%20%20_data%20%3D%20pd.read_csv(%0A%20%20%20%20%20%20%20%20repo_file('data%2Fballtoss.txt')%2C%0A%20%20%20%20%20%20%20%20sep%3D%22%5Ct%22%2C%0A%20%20%20%20%20%20%20%20header%3DNone%2C%0A%20%20%20%20%20%20%20%20skiprows%3D2%2C%0A%20%20%20%20).to_numpy()%0A%20%20%20%20ball_t%2C%20ball_x%2C%20ball_y%20%3D%20_data.T%0A%20%20%20%20ball_dt%20%3D%20np.mean(np.diff(ball_t))%0A%20%20%20%20print(f%22Time%20interval%3A%20%7Bball_dt%3A.4f%7D%20s%20(%7B1%20%2F%20ball_dt%3A.0f%7D%20Hz)%2C%20%7Bball_t.size%7D%20frames%22)%0A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(12%2C%203.2))%0A%20%20%20%20_axs%5B0%5D.plot(ball_x%2C%20ball_y%2C%20%22go%22)%0A%20%20%20%20_axs%5B0%5D.set_xlabel(%22x%20%5Bm%5D%22)%0A%20%20%20%20_axs%5B0%5D.set_ylabel(%22y%20%5Bm%5D%22)%0A%20%20%20%20_axs%5B1%5D.plot(ball_t%2C%20ball_x%2C%20%22bo%22)%0A%20%20%20%20_axs%5B1%5D.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_axs%5B1%5D.set_ylabel(%22x%20%5Bm%5D%22)%0A%20%20%20%20_axs%5B2%5D.plot(ball_t%2C%20ball_y%2C%20%22ro%22)%0A%20%20%20%20_axs%5B2%5D.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_axs%5B2%5D.set_ylabel(%22y%20%5Bm%5D%22)%0A%20%20%20%20plt.suptitle(%22Kinematics%20of%20a%20ball%20toss%22%2C%20fontsize%3D14)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20ball_dt%2C%20ball_t%2C%20ball_x%2C%20ball_y%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Twenty-two%20frames%20at%2030%20Hz.%20Let's%20compute%20the%20velocity%20and%20acceleration%20numerically%2C%20with%20two%20algorithms%3A%20the%20forward%20difference%2C%20%24(x_%7Bi%2B1%7D%20-%20x_i)%2F%5CDelta%20t%24%2C%20and%20the%20central%20difference%2C%20%24(x_%7Bi%2B1%7D%20-%20x_%7Bi-1%7D)%2F(2%5CDelta%20t)%24.%0A%0A%20%20%20%20**Before%20you%20run%20the%20next%20cells**%2C%20predict%20the%20horizontal%20and%20vertical%20accelerations%20of%20the%20ball.%20Which%20of%20the%20two%20algorithms%20do%20you%20expect%20to%20be%20noisier%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(ball_dt%2C%20ball_x%2C%20ball_y%2C%20np)%3A%0A%20%20%20%20%23%20forward%20difference%20algorithm%0A%20%20%20%20vx%2C%20vy%20%3D%20np.diff(ball_x)%20%2F%20ball_dt%2C%20np.diff(ball_y)%20%2F%20ball_dt%0A%20%20%20%20ax%2C%20ay%20%3D%20np.diff(vx)%20%2F%20ball_dt%2C%20np.diff(vy)%20%2F%20ball_dt%0A%20%20%20%20%23%20central%20difference%20algorithm%0A%20%20%20%20vx2%20%3D%20(ball_x%5B2%3A%5D%20-%20ball_x%5B%3A-2%5D)%20%2F%20(2%20*%20ball_dt)%0A%20%20%20%20vy2%20%3D%20(ball_y%5B2%3A%5D%20-%20ball_y%5B%3A-2%5D)%20%2F%20(2%20*%20ball_dt)%0A%20%20%20%20ax2%2C%20ay2%20%3D%20(vx2%5B2%3A%5D%20-%20vx2%5B%3A-2%5D)%20%2F%20(2%20*%20ball_dt)%2C%20(vy2%5B2%3A%5D%20-%20vy2%5B%3A-2%5D)%20%2F%20(2%20*%20ball_dt)%0A%20%20%20%20return%20ax%2C%20ax2%2C%20ay%2C%20ay2%2C%20vx%2C%20vx2%2C%20vy%2C%20vy2%0A%0A%0A%40app.cell%0Adef%20_(ax%2C%20ax2%2C%20ay%2C%20ay2%2C%20ball_t%2C%20ball_x%2C%20ball_y%2C%20plt%2C%20vx%2C%20vx2%2C%20vy%2C%20vy2)%3A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(2%2C%203%2C%20sharex%3DTrue%2C%20figsize%3D(11%2C%206))%0A%20%20%20%20_axs%5B0%2C%200%5D.plot(ball_t%2C%20ball_x%2C%20%22bo%22)%0A%20%20%20%20_axs%5B0%2C%200%5D.set_ylabel(%22x%20%5Bm%5D%22)%0A%20%20%20%20_axs%5B0%2C%201%5D.plot(ball_t%5B%3A-1%5D%2C%20vx%2C%20%22bo%22%2C%20label%3D%22forward%20difference%22)%0A%20%20%20%20_axs%5B0%2C%201%5D.plot(ball_t%5B1%3A-1%5D%2C%20vx2%2C%20%22m%2B%22%2C%20markersize%3D10%2C%20label%3D%22central%20difference%22)%0A%20%20%20%20_axs%5B0%2C%201%5D.set_ylabel(%22vx%20%5Bm%2Fs%5D%22)%0A%20%20%20%20_axs%5B0%2C%201%5D.legend(frameon%3DFalse%2C%20fontsize%3D10%2C%20loc%3D%22upper%20left%22%2C%20numpoints%3D1)%0A%20%20%20%20_axs%5B0%2C%202%5D.plot(ball_t%5B%3A-2%5D%2C%20ax%2C%20%22bo%22)%0A%20%20%20%20_axs%5B0%2C%202%5D.plot(ball_t%5B2%3A-2%5D%2C%20ax2%2C%20%22m%2B%22%2C%20markersize%3D10)%0A%20%20%20%20_axs%5B0%2C%202%5D.set_ylabel(%22ax%20%5Bm%2Fs%24%5E2%24%5D%22)%0A%20%20%20%20_axs%5B1%2C%200%5D.plot(ball_t%2C%20ball_y%2C%20%22ro%22)%0A%20%20%20%20_axs%5B1%2C%200%5D.set_ylabel(%22y%20%5Bm%5D%22)%0A%20%20%20%20_axs%5B1%2C%201%5D.plot(ball_t%5B%3A-1%5D%2C%20vy%2C%20%22ro%22)%0A%20%20%20%20_axs%5B1%2C%201%5D.plot(ball_t%5B1%3A-1%5D%2C%20vy2%2C%20%22m%2B%22%2C%20markersize%3D10)%0A%20%20%20%20_axs%5B1%2C%201%5D.set_ylabel(%22vy%20%5Bm%2Fs%5D%22)%0A%20%20%20%20_axs%5B1%2C%202%5D.plot(ball_t%5B%3A-2%5D%2C%20ay%2C%20%22ro%22)%0A%20%20%20%20_axs%5B1%2C%202%5D.plot(ball_t%5B2%3A-2%5D%2C%20ay2%2C%20%22m%2B%22%2C%20markersize%3D10)%0A%20%20%20%20_axs%5B1%2C%202%5D.set_ylabel(%22ay%20%5Bm%2Fs%24%5E2%24%5D%22)%0A%20%20%20%20for%20_ax%20in%20_axs%5B1%5D%3A%0A%20%20%20%20%20%20%20%20_ax.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20plt.suptitle(%22Kinematics%20of%20a%20ball%20toss%22%2C%20fontsize%3D14)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20noise%20shows%20up%20in%20the%20derivatives%2C%20as%20always.%20Once%20in%20the%20air%2C%20the%20ball%20has%20only%20gravity%20acting%20on%20it%20(neglecting%20air%20resistance)%2C%20so%20its%20horizontal%20acceleration%20should%20be%20zero%20and%20its%20vertical%20acceleration%20constant%2C%20about%20%24-9.8%24%20m%2Fs%C2%B2.%20The%20forward-difference%20accelerations%20scatter%20around%20those%20values%20by%20about%201%20m%2Fs%C2%B2%3B%20the%20central%20differences%2C%20which%20here%20average%20over%20a%20wider%20interval%2C%20scatter%20less.%0A%0A%20%20%20%20To%20estimate%20the%20acceleration%20we%20could%20filter%20the%20data%2C%20but%20with%2022%20frames%20there%20is%20little%20room%20for%20a%20filter%2C%20and%20we%20know%20the%20physics%20of%20the%20phenomenon.%20So%20we%20can%20instead%20fit%20a%20*model*%20to%20the%20data.%20For%20the%20vertical%20position%2C%20the%20model%20is%20that%20of%20a%20particle%20at%20constant%20acceleration%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20y(t)%20%3D%20y_0%20%2B%20v_0%20t%20%2B%20%5Cfrac%7B1%7D%7B2%7D%20g%20t%5E2%0A%20%20%20%20%24%24%0A%0A%20%20%20%20a%20second-order%20polynomial%20whose%20leading%20coefficient%20is%20%24g%2F2%24.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20how%20close%20to%20%24-9.8%24%20m%2Fs%C2%B2%20the%20fitted%20%24g%24%20will%20be.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(ball_t%2C%20ball_y%2C%20np)%3A%0A%20%20%20%20_p%20%3D%20np.polyfit(ball_t%2C%20ball_y%2C%202)%0A%20%20%20%20print(f%22g%20%3D%20%7B2%20*%20_p%5B0%5D%3A.2f%7D%20m%2Fs2%20(all%20%7Bball_t.size%7D%20frames)%22)%0A%20%20%20%20_p%20%3D%20np.polyfit(ball_t%5B%3A-4%5D%2C%20ball_y%5B%3A-4%5D%2C%202)%0A%20%20%20%20print(f%22g%20%3D%20%7B2%20*%20_p%5B0%5D%3A.2f%7D%20m%2Fs2%20(without%20the%20last%204%20frames)%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20A%20good%20estimate%2C%20from%20data%20whose%20point-by-point%20accelerations%20scatter%20by%20a%20whole%20m%2Fs%C2%B2.%20Fitting%20a%20model%20uses%20all%20the%20frames%20at%20once%20to%20estimate%20one%20parameter%2C%20which%20is%20why%20it%20is%20so%20much%20more%20precise%20than%20any%20derivative.%20The%20price%20is%20that%20you%20must%20trust%20the%20model.%0A%0A%20%20%20%20Look%20again%20at%20the%20central-difference%20vertical%20acceleration%3A%20it%20drifts%20from%20about%20%24-9.7%24%20to%20almost%20%24-11%24%20m%2Fs%C2%B2%20over%20the%20last%20frames%2C%20which%20no%20ball%20in%20free%20flight%20does.%20Something%20is%20wrong%20with%20the%20last%20frames%2C%20perhaps%20a%20distortion%20of%20the%20lens%20near%20the%20edge%20of%20the%20video%20image%2C%20or%20the%20ball%20moving%20out%20of%20the%20plane%20of%20the%20camera.%20Leaving%20out%20the%20last%20four%20frames%20changes%20%24g%24%20by%20about%200.1%20m%2Fs%C2%B2.%20To%20read%20more%20about%20fitting%20a%20model%20to%20data%2C%20see%20%5BCurve%20fitting%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fcurve_fitting.py).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20The%20optimal%20cutoff%20frequency%0A%0A%20%20%20%20You%20are%20probably%20now%20wondering%20how%20to%20determine%20automatically%20the%20optimal%20cutoff%20frequency%20of%20a%20low-pass%20filter%2C%20the%20one%20that%20attenuates%20as%20much%20of%20the%20noise%20as%20possible%20without%20compromising%20the%20signal.%20This%20is%20an%20important%20topic%20in%20signal%20processing%2C%20particularly%20in%20movement%20science%2C%20and%20one%20method%20for%20it%20is%20the%20subject%20of%20the%20notebook%20%5BResidual%20analysis%20to%20determine%20the%20optimal%20cutoff%20frequency%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fresidual_analysis.py).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Exercise%208.**%20Why%20can%20a%20quadratic%20fit%20estimate%20gravity%20from%20a%20short%20ball%20flight%20more%20reliably%20than%20two%20numerical%20differences%3F%20What%20would%20invalidate%20that%20fit%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%208%0A%20%20%20%20A%20quadratic%20uses%20the%20constant-acceleration%20model%20to%20combine%20all%20frames%20into%20three%20coefficients%3B%20its%20acceleration%20is%20the%20constant%20%242p_2%24.%20Differencing%20instead%20magnifies%20local%20tracking%20errors.%20The%20model%20fails%20if%20drag%2C%20camera%20perspective%2C%20out-of-plane%20motion%20or%20contact%20forces%20matter.%20A%20small%20fit%20residual%20alone%20cannot%20establish%20the%20model%E2%80%99s%20physical%20validity.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Checkpoint%20questions%0A%0A%20%20%20%20Pause%20here%20before%20the%20problems.%0A%0A%20%20%20%201.%20Go%20back%20to%20your%20guess%20in%20Challenge%200.%20By%20what%20factor%20does%20each%20differentiation%20multiply%20a%20noise%20component%20at%20frequency%20%24f%24%3F%20What%20does%20that%20mean%20for%20noise%20near%20the%20Nyquist%20frequency%20of%20a%20100%20Hz%20recording%3F%0A%20%20%20%202.%20A%20single-pass%20Butterworth%20filter%20at%206%20Hz%20is%20applied%20to%20a%20marker%20trajectory%20but%20not%20to%20the%20EMG%20recorded%20with%20it.%20What%20happens%20to%20the%20timing%20between%20muscle%20activation%20and%20movement%3F%0A%20%20%20%203.%20Why%20must%20the%20cutoff%20frequency%20be%20corrected%20when%20filtering%20with%20%60filtfilt%60%3F%20What%20do%20you%20report%20in%20your%20methods%20if%20you%20forget%3F%0A%20%20%20%204.%20When%20would%20you%20prefer%20a%20critically%20damped%20filter%20to%20a%20Butterworth%20filter%3F%0A%20%20%20%205.%20Why%20can%20%60filtfilt%60%20not%20be%20used%20to%20filter%20data%20in%20real%20time%2C%20for%20instance%20in%20a%20biofeedback%20system%3F%0A%20%20%20%206.%20The%20ball%20toss%20was%20solved%20by%20fitting%20a%20model%2C%20the%20Pezzack%20data%20by%20filtering.%20When%20is%20each%20approach%20appropriate%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Checkpoint%20answers%0A%20%20%20%201.%20Continuous%20differentiation%20multiplies%20amplitude%20by%20%242%5Cpi%20f%24.%20For%20the%20discrete%20second%20difference%2C%20the%20exact%20magnitude%20is%20%244f_s%5E2%5Csin%5E2(%5Cpi%20f%2Ff_s)%24%3B%20it%20is%20large%20near%20Nyquist%20and%20only%20approaches%20%24(2%5Cpi%20f)%5E2%24%20at%20low%20frequency.%0A%20%20%20%202.%20The%20filtered%20marker%20trajectory%20lags%20by%20a%20frequency-dependent%20amount%2C%20biasing%20comparison%20with%20the%20unfiltered%20EMG.%0A%20%20%20%203.%20Two%20passes%20square%20the%20amplitude%20response.%20Report%20the%20design%20cutoff%2C%20order%20per%20pass%2C%20number%20of%20passes%20and%20effective%20response%3B%20do%20not%20label%20the%20uncorrected%20design%20frequency%20as%20the%20combined%20%24-3%24%20dB%20cutoff.%0A%20%20%20%204.%20When%20avoiding%20step%20ringing%20matters%20more%20than%20a%20sharp%20frequency%20transition%2C%20after%20checking%20the%20actual%20digital%20response.%0A%20%20%20%205.%20The%20backward%20pass%20needs%20future%20samples%20from%20the%20whole%20record.%0A%20%20%20%206.%20Fit%20a%20model%20when%20its%20physical%20assumptions%20are%20justified%3B%20use%20filtering%20when%20an%20appropriate%20signal%2Fnoise%20separation%20is%20defensible%20without%20prescribing%20the%20whole%20trajectory.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Problems%0A%0A%20%20%20%201.%20Show%20that%20the%20coefficients%20of%20the%20moving-average%20filter%20satisfy%20%24%5Csum%20b_k%20-%20%5Csum_%7Bk%20%5Cgeq%201%7D%20a_k%20%3D%201%24.%20What%20does%20this%20property%20mean%20for%20a%20constant%20input%3F%0A%0A%20%20%20%202.%20Use%20%60signal.freqz%60%20to%20plot%20the%20frequency%20response%20of%20moving-average%20filters%20of%205%20and%2011%20samples%2C%20for%20data%20sampled%20at%20100%20Hz.%20At%20what%20frequency%20is%20the%20first%20zero%20of%20each%20response%3F%20Is%20the%20moving%20average%20a%20good%20low-pass%20filter%3F%20Compare%20it%20with%20a%20second-order%20Butterworth%20filter%20whose%20%24-3%24%20dB%20point%20matches%20the%20moving%20average%E2%80%99s%20%24-3%24%20dB%20point%20(not%20its%20first%20zero).%0A%0A%20%20%20%203.%20A%20high-pass%20filter%20attenuates%20low%20frequencies%20and%20passes%20high%20frequencies.%20Remove%20the%20DC%20component%20and%20any%20slow%20drift%20of%20the%20EMG%20above%20with%20a%20second-order%20Butterworth%20high-pass%20filter%20at%2020%20Hz.%20Compare%20the%20result%20with%20simply%20subtracting%20the%20mean.%0A%0A%20%20%20%204.%20The%20linear%20envelope%20of%20an%20EMG%20is%20obtained%20by%20full-wave%20rectifying%20it%20(taking%20its%20absolute%20value)%20and%20filtering%20it%20with%20a%20low-pass%20filter%20at%20a%20few%20hertz.%20Compute%20the%20linear%20envelope%20of%20the%20EMG%20above%20with%20a%20zero-phase%20Butterworth%20filter%20at%205%20Hz%2C%20and%20compare%20it%20with%20the%20moving%20RMS.%20See%20%5BElectromyography%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Felectromyography.py).%0A%0A%20%20%20%205.%20Filter%20the%20vertical%20position%20of%20the%20ball%20toss%20with%20a%20zero-phase%20Butterworth%20low-pass%20filter%2C%20differentiate%20it%20twice%2C%20and%20compare%20the%20acceleration%20with%20the%20%24g%24%20estimated%20by%20the%20polynomial%20fit.%20Which%20cutoff%20frequency%20did%20you%20use%2C%20and%20why%20is%20filtering%20harder%20here%20than%20with%20the%20Pezzack%20data%3F%0A%0A%20%20%20%206.%20Repeat%20the%20comparison%20of%20the%20Butterworth%2C%20Savitzky-Golay%20and%20spline%20methods%20with%20the%20noisy%20version%20of%20Pezzack's%20angle%2C%20%60pz_disp_noisy%60%2C%20adjusting%20the%20parameters%20of%20each%20method.%20Does%20the%20ranking%20of%20the%20methods%20change%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Final%20answers%0A%20%20%20%201.%20%24b_k%3D1%2Fm%24%20for%20%24k%3D0%2C%5Cldots%2Cm-1%24%20and%20all%20feedback%20coefficients%20are%20zero%2C%20so%20the%20expression%20is%201.%20A%20constant%20input%20is%20preserved%20away%20from%20initialization%2Fpadding%20boundaries.%0A%0A%20%20%20%202.%20First%20zeros%3A%2020%20Hz%20for%205%20samples%20and%20%24100%2F11%5Capprox9.0909%24%20Hz%20for%2011%20samples.%20Their%20single-pass%20%24-3%24%20dB%20points%20are%20approximately%209.0159%20and%204.0412%20Hz.%20A%20matched%20Butterworth%20response%20decreases%20monotonically%3B%20the%20moving%20average%20has%20sidelobes%20and%20passes%20some%20higher-frequency%20components%2C%20so%20it%20is%20a%20limited%20low-pass%20filter.%0A%0A%20%20%20%203.%20Mean%20subtraction%20removes%20only%20DC.%20A%2020%20Hz%20high-pass%20also%20attenuates%20slow%20fluctuations%20and%20movement%20artifacts%3B%20a%20single%20pass%20changes%20phase%2C%20and%20two%20uncorrected%20passes%20give%20%24-6.02%24%20dB%20at%2020%20Hz.%20High-frequency%20EMG%20remains%2C%20subject%20to%20the%20chosen%20response.%0A%0A%20%20%20%204.%20Both%20estimates%20follow%20the%20burst%20timing.%20The%20rectified-and-smoothed%20envelope%20estimates%20local%20mean%20absolute%20amplitude%3B%20RMS%20estimates%20the%20square%20root%20of%20local%20mean%20square%2C%20so%20their%20amplitudes%20are%20not%20interchangeable.%20For%20locally%20zero-mean%20Gaussian%20amplitude%20distributions%2C%20mean%20absolute%20amplitude%20is%20%24%5Csqrt%7B2%2F%5Cpi%7D%5Capprox0.798%24%20times%20RMS.%20Different%20windows%20and%20edge%20handling%20give%20further%20differences.%0A%0A%20%20%20%205.%20One%20defensible%20illustrative%20choice%20is%20a%20corrected%203%20Hz%20two-pass%20Butterworth%20filter.%20Its%20acceleration%20is%20not%20constant%20and%20its%20endpoints%20are%20especially%20distorted.%20The%20quadratic%20estimates%20are%20%24g%3D-9.9781%24%20m%2Fs%C2%B2%20for%20all%2022%20frames%20and%20%24g%3D-9.8720%24%20m%2Fs%C2%B2%20without%20the%20last%20four.%20This%20record%20is%20only%20about%200.7%20s%20long%2C%20so%20boundaries%20and%20model%20violations%20matter%20much%20more%20than%20in%20the%20longer%20Pezzack%20record%3B%20no%20unique%20data-only%20cutoff%20is%20justified.%0A%0A%20%20%20%206.%20A%20reproducible%20tuned%20comparison%20using%20RMSE%20on%20the%20common%20interior%20samples%2C%20odd%20padding%20of%2011%20samples%20for%20polynomial%20methods%2C%20and%20the%20grids%20specified%20here%20gives%20approximately%3A%20Butterworth%20cutoff%205.1%20Hz%2C%20RMSE%204.51%20rad%2Fs%C2%B2%3B%20Savitzky%E2%80%93Golay%20window%2015%2Forder%204%2C%20RMSE%204.59%3B%20quintic%20spline%20%24s%5Capprox0.00545247%24%2C%20RMSE%204.24.%20For%20the%20cleaner%20angle%2C%20a%20tuned%20spline%20is%20also%20slightly%20best%20(about%204.05)%2C%20followed%20by%20Savitzky%E2%80%93Golay%20(4.14)%20and%20Butterworth%20(4.16).%20Thus%20this%20tuning%20does%20not%20reverse%20the%20ranking.%20These%20are%20benchmark-specific%20comparisons%2C%20not%20general%20superiority%20claims.%20The%20search%20grids%20are%20cutoff%202%E2%80%9318%20Hz%20in%200.1%20Hz%20steps%2C%20odd%20Savitzky%E2%80%93Golay%20windows%205%E2%80%9323%2Forders%202%E2%80%934%2C%20and%20spline%20%24s%3Dk%5C%2C%5Cmathrm%7Bvar%7D(x_%7Bpad%7D)%2FN_%7Bpad%7D%24%20with%2081%20logarithmically%20spaced%20%24k%24%20values%20from%200.01%20to%20100.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Go%20deeper%0A%0A%20%20%20%20-%20%5BResidual%20analysis%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fresidual_analysis.py)%20%E2%80%94%20choosing%20the%20cutoff%20frequency%20from%20the%20data.%0A%20%20%20%20-%20%5BBasic%20properties%20of%20signals%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fsignal_basic_properties.py)%20%E2%80%94%20sampling%2C%20aliasing%2C%20quantization%20and%20the%20signal-to-noise%20ratio.%0A%20%20%20%20-%20%5BFourier%20transform%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Ffourier_transform.py)%20%E2%80%94%20the%20frequency%20content%20of%20a%20signal.%0A%20%20%20%20-%20%5BCurve%20fitting%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fcurve_fitting.py)%20%E2%80%94%20fitting%20a%20model%20instead%20of%20filtering.%0A%20%20%20%20-%20%5BElectromyography%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Felectromyography.py)%20%E2%80%94%20filtering%2C%20rectifying%20and%20enveloping%20EMG.%0A%20%20%20%20-%20%5BThe%20Scientist%20and%20Engineer%E2%80%99s%20Guide%20to%20Digital%20Signal%20Processing%5D(https%3A%2F%2Fwww.dspguide.com%2F).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Video%20lectures%20on%20the%20Internet%0A%0A%20%20%20%20-%20MIT%20OpenCourseWare%2C%20Alan%20Oppenheim%3A%20%5BDigital%20Butterworth%20Filters%5D(https%3A%2F%2Focw.mit.edu%2Fcourses%2Fres-6-008-digital-signal-processing-spring-2011%2Fresources%2Flecture-16-digital-butterworth-filters%2F).%20Follow%20the%20bilinear-transform%20section%20to%20connect%20the%20analog%20response%20to%20the%20exact%20digital%20cutoff%20correction.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20References%0A%0A%20%20%20%20-%20Lyons%20RG%20(2010)%20%5BUnderstanding%20Digital%20Signal%20Processing%5D(https%3A%2F%2Fbooks.google.com.br%2Fbooks%3Fid%3DUBU7Y2tpwWUC%26hl).%203rd%20edition.%20Prentice%20Hall.%0A%0A%20%20%20%20-%20Pezzack%20JC%2C%20Norman%20RW%2C%20Winter%20DA%20(1977)%20%5BAn%20assessment%20of%20derivative%20determining%20techniques%20used%20for%20motion%20analysis%5D(https%3A%2F%2Fpubmed.ncbi.nlm.nih.gov%2F893476%2F).%20Journal%20of%20Biomechanics%2C%2010%2C%20377-382.%20%5BPubMed%5D(https%3A%2F%2Fpubmed.ncbi.nlm.nih.gov%2F893476%2F).%0A%0A%20%20%20%20-%20Robertson%20DG%2C%20Dowling%20JJ%20(2003)%20%5BDesign%20and%20responses%20of%20Butterworth%20and%20critically%20damped%20digital%20filters%5D(https%3A%2F%2Fpubmed.ncbi.nlm.nih.gov%2F14573371%2F).%20Journal%20of%20Electromyography%20and%20Kinesiology%2C%2013(6)%2C%20569-573.%0A%0A%20%20%20%20-%20Robertson%20G%2C%20Caldwell%20G%2C%20Hamill%20J%2C%20Kamen%20G%20(2013)%20%5BResearch%20Methods%20in%20Biomechanics%5D(https%3A%2F%2Fbooks.google.com.br%2Fbooks%3Fid%3DgRn8AAAAQBAJ).%202nd%20edition.%20Human%20Kinetics.%0A%0A%20%20%20%20-%20Vint%20PF%2C%20Hinrichs%20RN%20(1996)%20Endpoint%20error%20in%20smoothing%20and%20differentiating%20raw%20kinematic%20data%3A%20an%20evaluation%20of%20four%20popular%20methods.%20Journal%20of%20Biomechanics%2C%2029%2C%201637-1642.%20(The%20source%20of%20the%20Pezzack%20data%20file%20used%20here.)%0A%0A%20%20%20%20-%20Winter%20DA%20(2009)%20%5BBiomechanics%20and%20Motor%20Control%20of%20Human%20Movement%5D(https%3A%2F%2Fbooks.google.com.br%2Fbooks%3Fid%3D_bFHL08IWfwC).%204th%20edition.%20Hoboken%2C%20USA%3A%20Wiley.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A
5a2628a13ea2f5daa9fc19457c8bec51