residual analysis

%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%22optcutfreq%3E%3D0.0.10%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%20Residual%20analysis%20to%20determine%20the%20optimal%20cutoff%20frequency%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%20Every%20low-pass%20filter%20needs%20a%20cutoff%20frequency%2C%20and%20the%20choice%20matters%3A%20too%20high%20and%20the%20noise%20passes%20through%2C%20too%20low%20and%20the%20signal%20is%20distorted.%20A%20common%20problem%20in%20signal%20processing%20is%20to%20determine%20automatically%20the%20cutoff%20frequency%20that%20attenuates%20as%20much%20of%20the%20noise%20as%20possible%20without%20compromising%20the%20signal.%20There%20is%20no%20definitive%20solution%20to%20it%2C%20but%20there%20are%20techniques%2C%20with%20different%20degrees%20of%20success.%0A%0A%20%20%20%20This%20notebook%20builds%20one%20of%20them%2C%20proposed%20by%20David%20Winter%2C%20step%20by%20step%3A%20first%20the%20idea%2C%20then%20the%20computation%20by%20hand%2C%20then%20an%20automatic%20implementation.%20It%20then%20does%20what%20is%20rarely%20possible%20in%20practice%3A%20it%20checks%20the%20answer%20against%20an%20independently%20measured%20acceleration%20reference.%20That%20measurement%20is%20a%20benchmark%2C%20not%20an%20uncertainty-free%20truth.%0A%0A%20%20%20%20You%20should%20know%20the%20material%20of%20%5BData%20filtering%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fdata_filtering.py)%20first%2C%20in%20particular%20the%20zero-phase%20Butterworth%20filter%20and%20the%20correction%20of%20its%20cutoff%20frequency%20for%20two%20passes.%0A%0A%20%20%20%20Read%20it%20in%20order%20and%20run%20each%20cell%20as%20you%20reach%20it.%20Try%20the%20distributed%20Exercises%20before%20opening%20their%20worked%20solutions.%20Use%20the%20collapsed%20final%20answers%20to%20check%20the%20independent%20Problems%20at%20the%20end.%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.**%20In%20%5BData%20filtering%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fdata_filtering.py)%2C%20Pezzack's%20angle%20was%20filtered%20at%209%20Hz%20before%20being%20differentiated%2C%20and%20that%20value%20was%20picked%20by%20looking%20at%20the%20result.%20Suppose%20you%20had%20no%20accelerometer%20to%20compare%20with.%20Write%20down%2C%20in%20two%20or%20three%20sentences%2C%20how%20you%20would%20choose%20the%20cutoff%20frequency%20from%20the%20angle%20alone.%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%20Winter's%20idea%0A%0A%20%20%20%20David%20Winter%2C%20in%20his%20classic%20book%20*Biomechanics%20and%20Motor%20Control%20of%20Human%20Movement*%2C%20proposed%20finding%20the%20cutoff%20frequency%20from%20a%20residual%20analysis%3A%20filter%20the%20data%20with%20a%20range%20of%20cutoff%20frequencies%2C%20and%20for%20each%20one%20compute%20the%20residual%2C%20the%20RMS%20difference%20between%20the%20filtered%20and%20the%20unfiltered%20signals%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20R(f_c)%20%3D%20%5Csqrt%7B%5Cfrac%7B1%7D%7BN%7D%5Csum_%7Bi%3D1%7D%5E%7BN%7D%20%5Cleft(x_i%20-%20%5Chat%7Bx%7D_i(f_c)%5Cright)%5E2%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24x%24%20is%20the%20raw%20signal%20and%20%24%5Chat%7Bx%7D(f_c)%24%20is%20the%20signal%20filtered%20with%20the%20cutoff%20frequency%20%24f_c%24.%0A%0A%20%20%20%20The%20residual%20is%20whatever%20the%20filter%20removed.%20If%20the%20movement%20is%20concentrated%20at%20low%20frequencies%2C%20a%20sufficiently%20high%20cutoff%20may%20remove%20mostly%20noise.%20With%20broadband%20noise%2C%20a%20segment%20of%20the%20residual%20curve%20can%20look%20approximately%20straight%3B%20that%20straightness%20is%20an%20empirical%20local%20approximation%2C%20not%20a%20mathematical%20consequence%20of%20whiteness.%20At%20low%20cutoff%20frequencies%20the%20filter%20starts%20removing%20signal%20too%2C%20and%20the%20residual%20grows%20quickly.%20Winter's%20choice%20is%20the%20cutoff%20at%20the%20transition%2C%20where%20the%20residual%20starts%20to%20change%20very%20little%2C%20because%20from%20that%20point%20on%20the%20filter%20is%2C%20ideally%2C%20removing%20mostly%20noise%20and%20little%20signal.%0A%0A%20%20%20%20The%20concept%20is%20straightforward%20to%20implement.%20Let's%20do%20it%20on%20data%20where%20we%20can%20check%20the%20answer.%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.**%20A%20filter%20removes%20%24%5B0.002%2C-0.002%2C0.002%2C-0.002%5D%24%20radians%20from%20four%20samples.%20Compute%20the%20residual.%20Does%20it%20establish%20that%20the%20removed%20component%20was%20noise%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%201%0A%20%20%20%20The%20mean%20squared%20difference%20is%20%244%5Ctimes10%5E%7B-6%7D%24%20rad%C2%B2%2C%20giving%20%24R%3D0.002%24%20rad.%20The%20residual%20measures%20the%20size%20of%20what%20was%20removed%2C%20not%20its%20origin.%20A%20real%20movement%20component%20of%20the%20same%20amplitude%20gives%20exactly%20the%20same%20residual.%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%20Python%20setup%0A%0A%20%20%20%20NumPy%20and%20SciPy%20do%20the%20computation%2C%20pandas%20reads%20the%20data%20file%2C%20and%20Matplotlib%20draws%20the%20plots.%20%60optcutfreq%60%20implements%20the%20automatic%20residual%20criterion.%20All%20five%20packages%20are%20declared%20in%20this%20notebook%E2%80%99s%20dependency%20metadata.%20Open%20the%20file%20in%20marimo%3B%20the%20Pezzack%20benchmark%20requires%20access%20to%20the%20repository%20data.%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%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%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%20np%2C%20pd%2C%20plt%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%20Pezzack's%20benchmark%20data%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%3A%20the%20angle%20of%20a%20bar%20rotated%20by%20hand%2C%20digitized%20from%20film%2C%20and%20its%20angular%20acceleration%2C%20measured%20directly%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%20its%20own%20measurement%20uncertainty.%20The%20file%20has%20a%20second%2C%20noisier%20version%20of%20the%20angle%2C%20which%20we%20will%20use%20later.%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%20time%2C%20disp%2C%20disp_noisy%2C%20aacc%20%3D%20_data.T%0A%20%20%20%20freq%20%3D%201%20%2F%20np.mean(np.diff(time))%20%20%23%20sampling%20frequency%20%5BHz%5D%0A%0A%20%20%20%20_fig%2C%20(_ax1%2C%20_ax2)%20%3D%20plt.subplots(1%2C%202%2C%20sharex%3DTrue%2C%20figsize%3D(11%2C%204))%0A%20%20%20%20_ax1.plot(time%2C%20disp%2C%20%22b%22)%0A%20%20%20%20_ax1.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax1.set_ylabel(%22Angular%20displacement%20%5Brad%5D%22)%0A%20%20%20%20_ax2.plot(time%2C%20aacc%2C%20%22g%22)%0A%20%20%20%20_ax2.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax2.set_ylabel(%22Angular%20acceleration%20%5Brad%2Fs%24%5E2%24%5D%22)%0A%20%20%20%20plt.suptitle(%22Pezzack's%20benchmark%20data%22%2C%20fontsize%3D14)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20print(f%22Sampling%20frequency%3A%20%7Bfreq%3A.2f%7D%20Hz%20(Nyquist%20frequency%3A%20%7Bfreq%20%2F%202%3A.2f%7D%20Hz)%22)%0A%20%20%20%20return%20aacc%2C%20disp%2C%20disp_noisy%2C%20freq%2C%20time%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%20Step%201%3A%20the%20residuals%0A%0A%20%20%20%20Let's%20filter%20the%20angle%20with%20a%20second-order%20zero-phase%20Butterworth%20filter%20at%20100%20cutoff%20frequencies%2C%20from%20about%200.25%20Hz%20to%20just%20below%2020%20Hz%2C%20and%20compute%20the%20residual%20for%20each.%20The%20function%20below%20does%20it%3B%20the%20cutoff%20given%20to%20%60butter%60%20is%20corrected%20for%20the%20two%20passes%20of%20%60filtfilt%60%2C%20with%20%24C%20%3D%200.802%24.%20This%20is%20the%20historical%20low-frequency%20approximation%20used%20by%20%60optcutfreq%60%3B%20we%20retain%20it%20here%20to%20reproduce%20that%20package%E2%80%99s%20criterion.%20Its%20reported%20cutoffs%20are%20nominal%20values%20under%20this%20convention.%20The%20exact%20digital%20correction%2C%20explained%20in%20%5BData%20filtering%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fdata_filtering.py)%2C%20is%20%24f_%7Bdesign%7D%3D(f_s%2F%5Cpi)%5Carctan%5B%5Ctan(%5Cpi%20f_c%2Ff_s)%2FC%5D%24.%20The%20two%20conventions%20increasingly%20differ%20toward%20Nyquist%2C%20so%20report%20which%20one%20you%20use.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.function%0Adef%20residuals(y%2C%20freq%2C%20n%3D101)%3A%0A%20%20%20%20%22%22%22Residual%20RMS%20between%20%60y%60%20and%20%60y%60%20filtered%20at%20%60n%60%20cutoff%20frequencies.%0A%0A%20%20%20%20The%20filter%20is%20a%20second-order%20zero-phase%20Butterworth%20low-pass%20filter%2C%20with%0A%20%20%20%20the%20historical%20approximate%20cutoff%20correction%20used%20by%20optcutfreq.%20The%20cutoff%20frequencies%0A%20%20%20%20go%20from%201%25%20of%20the%20Nyquist%20frequency%20to%20just%20below%200.802%20times%20it.%0A%0A%20%20%20%20Returns%20(cutoff%20frequencies%20%5BHz%5D%2C%20residuals).%0A%20%20%20%20%22%22%22%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20from%20scipy.signal%20import%20butter%2C%20filtfilt%0A%0A%20%20%20%20C%20%3D%200.802%20%20%23%20historical%20approximate%20two-pass%20correction%3A%20C%20%3D%20(2**(1%2F2)%20-%201)**0.25%0A%20%20%20%20freqs%20%3D%20np.linspace((freq%20%2F%202)%20%2F%20100%2C%20(freq%20%2F%202)%20*%20C%2C%20n%2C%20endpoint%3DFalse)%0A%20%20%20%20res%20%3D%20np.empty(freqs.size)%0A%20%20%20%20for%20i%2C%20fc%20in%20enumerate(freqs)%3A%0A%20%20%20%20%20%20%20%20b%2C%20a%20%3D%20butter(2%2C%20(fc%20%2F%20C)%20%2F%20(freq%20%2F%202))%0A%20%20%20%20%20%20%20%20res%5Bi%5D%20%3D%20np.sqrt(np.mean((filtfilt(b%2C%20a%2C%20y)%20-%20y)%20**%202))%0A%20%20%20%20return%20freqs%2C%20res%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**Before%20you%20run%20the%20next%20cell**%2C%20sketch%20the%20residual%20against%20the%20cutoff%20frequency.%20Where%20is%20it%20largest%2C%20where%20is%20it%20smallest%2C%20and%20what%20does%20it%20look%20like%20in%20between%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(disp%2C%20freq%2C%20np%2C%20plt)%3A%0A%20%20%20%20_freqs%2C%20_res%20%3D%20residuals(disp%2C%20freq)%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(_freqs%2C%20_res%2C%20%22b.%22%2C%20markersize%3D8)%0A%20%20%20%20_ax.set_xlabel(%22Cutoff%20frequency%20%5BHz%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22Residual%20RMS%20%5Brad%5D%22)%0A%20%20%20%20_ax.set_title(%22Residual%20analysis%20of%20Pezzack's%20angle%22)%0A%20%20%20%20_ax.set_ylim(0%2C%200.01)%0A%20%20%20%20_ax.grid(True%2C%20linestyle%3D%22%3A%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20print(f%22Residual%20at%20%7B_freqs%5B0%5D%3A.2f%7D%20Hz%3A%20%7B_res%5B0%5D%3A.4f%7D%20rad%3B%20at%20%7B_freqs%5B-1%5D%3A.2f%7D%20Hz%3A%20%7B_res%5B-1%5D%3A.6f%7D%20rad%22)%0A%20%20%20%20print(f%22Residual%20at%20%7B_freqs%5Bnp.argmin(np.abs(_freqs%20-%2010))%5D%3A.2f%7D%20Hz%3A%20%7B_res%5Bnp.argmin(np.abs(_freqs%20-%2010))%5D%3A.4f%7D%20rad%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%20plot%20is%20zoomed%20in%3A%20at%20the%20lowest%20cutoff%20frequencies%20the%20residual%20is%20far%20above%20the%20top%20of%20the%20axis%2C%20because%20there%20the%20filter%20removes%20most%20of%20the%20movement.%20Below%20about%205%20Hz%20the%20residual%20falls%20steeply%3B%20above%20about%208%20Hz%20it%20falls%20slowly%20and%20almost%20linearly%20towards%20zero%20at%20the%20highest%20cutoff%2C%20where%20the%20filter%20removes%20almost%20nothing.%0A%0A%20%20%20%20Interpret%20that%20approximately%20straight%20tail%20as%20mostly%20removed%20noise%2C%20subject%20to%20the%20bandwidth%20assumptions.%20Extrapolate%20it%20back%20to%20a%20cutoff%20of%200%20Hz%2C%20and%20its%20intercept%2C%20call%20it%20%24a%24%2C%20estimates%20the%20residual%20that%20would%20be%20left%20if%20the%20filter%20removed%20only%20noise%2C%20all%20of%20it%3A%20approximately%20the%20RMS%20of%20the%20noise.%20Winter's%20choice%20of%20cutoff%20frequency%20is%20the%20one%20at%20which%20the%20residual%20curve%20reaches%20that%20value%2C%20%24R(f_c)%20%3D%20a%24.%20It%20is%20a%20compromise%20between%20the%20signal%20distortion%20the%20filter%20introduces%20and%20the%20noise%20it%20lets%20through.%0A%0A%20%20%20%20**Guiding%20questions%201.**%0A%0A%20%20%20%201.%20By%20eye%2C%20where%20would%20you%20draw%20the%20straight%20line%20through%20the%20noisy%20part%2C%20and%20where%20does%20it%20cross%20the%20vertical%20axis%3F%0A%20%20%20%202.%20Where%20does%20the%20curve%20reach%20the%20height%20of%20that%20intercept%3F%20That%20is%20your%20residual-analysis%20cutoff%20frequency.%0A%20%20%20%203.%20The%20intercept%20is%20in%20radians.%20Convert%20it%20to%20degrees.%20Is%20that%20a%20plausible%20error%20for%20an%20angle%20digitized%20from%20film%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%202.**%20Answer%20Guiding%20questions%201%3A%20use%20the%20nearly%20straight%20tail%20around%208%E2%80%9318%20Hz%20to%20estimate%20the%20noise%20level%20and%20the%20crossing.%20Convert%20the%20noise%20estimate%20to%20degrees.%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%202%0A%20%20%20%20The%20automatic%20fit%20gives%20%24a%5Capprox0.001829%24%20rad%20and%20a%20crossing%20near%205.59%20Hz.%20Multiplying%20by%20%24180%2F%5Cpi%24%20gives%20about%20%240.105%5E%5Ccirc%24.%20A%20visual%20fit%20can%20differ%20modestly.%20Such%20a%20small%20angular%20uncertainty%20is%20plausible%20for%20digitized%20film%2C%20but%20plausibility%20does%20not%20independently%20verify%20the%20noise%20model.%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%20Step%202%3A%20an%20automatic%20search%0A%0A%20%20%20%20Doing%20this%20by%20eye%20is%20easy%20for%20one%20signal%20and%20impossible%20for%20a%20thousand.%20The%20function%20%60optcutfreq%60%2C%20from%20the%20Python%20package%20%5Boptcutfreq%5D(https%3A%2F%2Fpypi.org%2Fproject%2Foptcutfreq%2F)%20(included%20in%20the%20notebook%20dependencies)%2C%20automates%20it%20in%20three%20parts%3A%0A%0A%20%20%20%201.%20it%20computes%20the%20residuals%20over%20a%20range%20of%20cutoff%20frequencies%2C%20as%20%60residuals%60%20does%20above%3B%0A%20%20%20%202.%20it%20finds%20the%20noisy%20region%20of%20the%20residual%20curve.%20For%20that%2C%20it%20fits%20an%20exponential%20decay%20to%20the%20residuals%20and%20considers%20that%20the%20noisy%2C%20linear%20tail%20starts%20after%20three%20lifetimes%20of%20the%20decay%20(a%2095%25%20drop)%3B%20a%20straight%20line%20is%20then%20fitted%20to%20the%20residuals%20from%20that%20point%20on%2C%20and%20its%20intercept%20%24a%24%20is%20the%20noise%20level%3B%0A%20%20%20%203.%20it%20finds%20the%20cutoff%20frequency%20at%20which%20the%20residual%20curve%20equals%20%24a%24%2C%20by%20interpolating%20the%20residuals%20with%20a%20spline%2C%20and%20optionally%20plots%20the%20results.%0A%0A%20%20%20%20Winter%20should%20not%20be%20blamed%20for%20the%20automatic%20search.%20It%20only%20follows%2C%20as%20closely%20as%20possible%2C%20his%20suggestion%20of%20fitting%20a%20regression%20line%20to%20the%20noisy%20part%20of%20the%20residuals.%20If%20the%20search%20fails%2C%20the%20frequency%20limits%20of%20the%20noisy%20part%20can%20be%20given%20by%20hand%2C%20with%20the%20parameter%20%60fclim%60.%0A%0A%20%20%20%20Its%20signature%20is%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20fc_opt%20%3D%20optcutfreq(y%2C%20freq%3D1%2C%20fclim%3DNone%2C%20show%3DFalse%2C%20ax%3DNone)%0A%20%20%20%20%60%60%60%0A%0A%20%20%20%20Run%20%60help(optcutfreq)%60%20to%20read%20its%20documentation%2C%20or%20read%20its%20code%20on%20%5BGitHub%5D(https%3A%2F%2Fgithub.com%2Fdemotu%2Foptcutfreq).%20The%20code%20is%20long%20relative%20to%20the%20simplicity%20of%20the%20idea%2C%20because%20of%20the%20help%20section%2C%20the%20automatic%20search%2C%20and%20a%20rich%20plot.%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%20from%20optcutfreq%20import%20optcutfreq%0A%0A%20%20%20%20return%20(optcutfreq%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%20%23%23%20Step%203%3A%20the%20optimal%20cutoff%20frequency%20of%20Pezzack's%20angle%0A%0A%20%20%20%20**Before%20you%20run%20the%20next%20cell**%2C%20compare%20your%20answer%20to%20Guiding%20question%201.2%20with%20the%209%20Hz%20used%20in%20%5BData%20filtering%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fdata_filtering.py).%20Will%20the%20automatic%20search%20give%20a%20higher%20or%20a%20lower%20cutoff%20frequency%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(disp%2C%20freq%2C%20optcutfreq)%3A%0A%20%20%20%20fc_opt%20%3D%20optcutfreq(disp%2C%20freq%3Dfreq%2C%20show%3DTrue)%0A%20%20%20%20print(f%22Optimal%20cutoff%20frequency%3A%20%7Bfc_opt%3A.2f%7D%20Hz%22)%0A%20%20%20%20return%20(fc_opt%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%20optimal%20cutoff%20frequency%20found%20is%205.6%20Hz%2C%20lower%20than%20the%209%20Hz%20chosen%20by%20hand.%20The%20noise%20level%20%24a%24%2C%20the%20RMSE%20shown%20in%20the%20title%20of%20the%20top%20right%20panel%2C%20is%20about%200.0018%20rad%2C%20a%20tenth%20of%20a%20degree.%0A%0A%20%20%20%20Note%20that%20the%20filtering%20matters%20only%20for%20the%20derivatives%3A%20in%20the%20top%20right%20panel%20the%20unfiltered%20and%20filtered%20angles%20cannot%20be%20told%20apart.%20In%20the%20bottom%20right%20panel%2C%20their%20second%20derivatives%20can.%0A%0A%20%20%20%20Let's%20filter%20the%20angle%20at%20this%20cutoff%20frequency%2C%20differentiate%20it%20twice%2C%20and%20compare%20the%20result%20with%20the%20reference%20acceleration%2C%20as%20in%20%5BData%20filtering%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fdata_filtering.py).%20The%20function%20below%20does%20the%20filtering%20and%20the%20double%20differentiation%3B%20its%20output%20is%20two%20samples%20shorter%20than%20its%20input.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.function%0Adef%20filtered_acceleration(y%2C%20freq%2C%20fc)%3A%0A%20%20%20%20%22%22%22Second%20derivative%20of%20%60y%60%20after%20a%20zero-phase%20Butterworth%20filter%20at%20%60fc%60%20Hz.%0A%0A%20%20%20%20The%20filter%20is%20second%20order%2C%20applied%20forward%20and%20backward%20with%20its%20cutoff%0A%20%20%20%20corrected%20for%20the%20two%20passes%2C%20and%20the%20derivative%20is%20the%20second-order%0A%20%20%20%20finite%20difference%2C%20so%20the%20output%20is%20two%20samples%20shorter%20than%20%60y%60.%0A%20%20%20%20%22%22%22%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20from%20scipy.signal%20import%20butter%2C%20filtfilt%0A%0A%20%20%20%20C%20%3D%200.802%20%20%23%20historical%20approximate%20two-pass%20correction%0A%20%20%20%20b%2C%20a%20%3D%20butter(2%2C%20(fc%20%2F%20C)%20%2F%20(freq%20%2F%202))%0A%20%20%20%20return%20np.diff(filtfilt(b%2C%20a%2C%20y)%2C%202)%20*%20freq**2%0A%0A%0A%40app.cell%0Adef%20_(aacc%2C%20disp%2C%20fc_opt%2C%20freq%2C%20np%2C%20plt%2C%20time)%3A%0A%20%20%20%20_aacc_bw%20%3D%20filtered_acceleration(disp%2C%20freq%2C%20fc_opt)%0A%20%20%20%20_rmse%20%3D%20np.sqrt(np.mean((_aacc_bw%20-%20aacc%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))%0A%20%20%20%20_ax.plot(time%5B1%3A-1%5D%2C%20aacc%5B1%3A-1%5D%2C%20%22g%22%2C%20linewidth%3D3%2C%20label%3D%22Accelerometer%20(true%20value)%22)%0A%20%20%20%20_ax.plot(time%5B1%3A-1%5D%2C%20_aacc_bw%2C%20%22r%22%2C%20label%3Df%22Butterworth%20%7Bfc_opt%3A.3g%7D%20Hz%3A%20RMSE%20%3D%20%7B_rmse%3A.2f%7D%22)%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%20An%20RMSE%20of%20about%204.2%20rad%2Fs%C2%B2%2C%20slightly%20better%20than%204.45%20rad%2Fs%C2%B2%20at%20a%20nominal%209%20Hz%20under%20the%20same%20historical%20correction%20convention%2C%20and%20found%20without%20looking%20at%20the%20reference%20acceleration.%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%20Is%20the%20optimal%20cutoff%20frequency%20optimal%3F%0A%0A%20%20%20%20%22Optimal%22%20here%20means%20optimal%20by%20Winter's%20criterion%2C%20a%20compromise%20judged%20on%20the%20angle.%20What%20we%20actually%20care%20about%20is%20the%20acceleration.%20Because%20we%20have%20a%20reference%20acceleration%2C%20we%20can%20sweep%20a%20grid%20of%20cutoffs%20and%20compare%20the%20estimated%20accelerations%20with%20that%20reference.%20This%20check%20is%20unavailable%20when%20the%20experiment%20records%20only%20position.%0A%0A%20%20%20%20The%20residual%20criterion%20can%20oversmooth%20for%20a%20derivative-based%20task%2C%20but%20that%20is%20not%20guaranteed%20for%20every%20signal.%20The%20measured%20benchmark%20lets%20us%20check%20rather%20than%20assume%20the%20direction%20of%20the%20difference.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20where%20the%20best%20cutoff%20frequency%20for%20the%20acceleration%20will%20be%3A%20below%205.6%20Hz%2C%20above%20it%2C%20or%20above%209%20Hz%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(aacc%2C%20disp%2C%20fc_opt%2C%20freq%2C%20np%2C%20plt)%3A%0A%20%20%20%20fc_sweep%20%3D%20np.arange(2%2C%2018.01%2C%200.25)%0A%0A%0A%20%20%20%20def%20rmse_sweep(y)%3A%0A%20%20%20%20%20%20%20%20%22%22%22RMSE%20of%20the%20acceleration%20of%20%60y%60%20filtered%20at%20each%20of%20%60fc_sweep%60.%22%22%22%0A%20%20%20%20%20%20%20%20return%20np.array(%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Bnp.sqrt(np.mean((filtered_acceleration(y%2C%20freq%2C%20_fc)%20-%20aacc%5B1%3A-1%5D)%20**%202))%20for%20_fc%20in%20fc_sweep%5D%0A%20%20%20%20%20%20%20%20)%0A%0A%0A%20%20%20%20_rmse%20%3D%20rmse_sweep(disp)%0A%20%20%20%20_best%20%3D%20fc_sweep%5Bnp.argmin(_rmse)%5D%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(fc_sweep%2C%20_rmse%2C%20%22b.-%22)%0A%20%20%20%20_ax.axvline(fc_opt%2C%20color%3D%22r%22%2C%20linestyle%3D%22--%22%2C%20label%3Df%22residual%20analysis%3A%20%7Bfc_opt%3A.1f%7D%20Hz%22)%0A%20%20%20%20_ax.axvline(_best%2C%20color%3D%22g%22%2C%20linestyle%3D%22--%22%2C%20label%3Df%22best%20for%20the%20acceleration%3A%20%7B_best%3A.1f%7D%20Hz%22)%0A%20%20%20%20_ax.set_xlabel(%22Cutoff%20frequency%20%5BHz%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22RMSE%20of%20the%20acceleration%20%5Brad%2Fs%24%5E2%24%5D%22)%0A%20%20%20%20_ax.set_ylim(0%2C%2012)%0A%20%20%20%20_ax.legend(frameon%3DFalse%2C%20loc%3D%22upper%20center%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20for%20_fc%20in%20(fc_opt%2C%20_best%2C%209)%3A%0A%20%20%20%20%20%20%20%20print(f%22Cutoff%20%7B_fc%3A5.2f%7D%20Hz%3A%20RMSE%20%3D%20%7Bnp.sqrt(np.mean((filtered_acceleration(disp%2C%20freq%2C%20_fc)%20-%20aacc%5B1%3A-1%5D)%20**%202))%3A.2f%7D%20rad%2Fs2%22)%0A%20%20%20%20return%20fc_sweep%2C%20rmse_sweep%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%20best%20cutoff%20frequency%20for%20the%20acceleration%20is%20about%206.5%20Hz%2C%20a%20little%20above%20the%205.6%20Hz%20of%20the%20residual%20analysis%3A%20it%20oversmoothed%20slightly%20for%20this%20particular%20acceleration%20benchmark.%20But%20look%20at%20the%20shape%20of%20the%20curve.%20It%20is%20flat%20around%20its%20minimum%2C%20and%20the%20price%20of%20the%20oversmoothing%20is%20tiny%2C%20an%20RMSE%20of%204.22%20instead%20of%204.16%20rad%2Fs%C2%B2.%20Cutoff%20frequencies%20that%20are%20too%20low%20are%20punished%20hard%20(at%202%20Hz%20the%20error%20is%20more%20than%20twice%20the%20minimum)%3B%20too%20high%2C%20the%20error%20grows%20slowly.%0A%0A%20%20%20%20**Guiding%20questions%202.**%0A%0A%20%20%20%201.%20Why%20is%20the%20curve%20so%20asymmetric%3F%20Think%20of%20what%20is%20removed%20below%20and%20above%20the%20best%20cutoff%20frequency.%0A%20%20%20%202.%20The%20residual%20analysis%20judged%20the%20cutoff%20on%20the%20angle%2C%20but%20we%20evaluated%20it%20on%20the%20acceleration.%20Would%20the%20best%20cutoff%20frequency%20for%20the%20angular%20*velocity*%20be%20higher%20or%20lower%20than%20for%20the%20acceleration%3F%0A%20%20%20%203.%20In%20a%20real%20experiment%20you%20cannot%20draw%20this%20curve.%20What%20does%20its%20flatness%20near%20the%20minimum%20tell%20you%20about%20how%20precise%20a%20cutoff%20frequency%20needs%20to%20be%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.**%20The%20residual%20criterion%20selects%205.59%20Hz%2C%20while%20the%20acceleration%20RMSE%20minimum%20is%206.5%20Hz.%20Is%20either%20answer%20%E2%80%9Cwrong%E2%80%9D%3F%20Would%20a%20velocity%20optimum%20necessarily%20lie%20between%20them%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%203%0A%20%20%20%20They%20optimize%20different%20criteria%3A%20one%20uses%20the%20angle%E2%80%99s%20residual%20curve%3B%20the%20other%20uses%20a%20measured%20acceleration%20reference.%20The%20latter%20has%20RMSE%20about%204.16%20rather%20than%204.22%20rad%2Fs%C2%B2%2C%20a%20small%20difference%20here.%20Differentiation%20changes%20the%20weighting%20of%20noise%20and%20signal%20across%20frequency.%20There%20is%20no%20general%20theorem%20that%20the%20velocity%20optimum%20lies%20between%20those%20two%20numbers%3B%20calculate%20it%20if%20a%20suitable%20reference%20is%20available.%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%20Noisier%20data%0A%0A%20%20%20%20The%20data%20file%20has%20a%20second%20version%20of%20the%20angle%20with%20more%20noise%20added%20to%20it.%20A%20fixed%20cutoff%20frequency%20ignores%20the%20noise%20level%3B%20the%20residual%20analysis%20does%20not.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20whether%20the%20optimal%20cutoff%20frequency%20for%20the%20noisier%20angle%20will%20be%20higher%20or%20lower%20than%205.6%20Hz%2C%20and%20whether%20the%209%20Hz%20chosen%20by%20hand%20for%20the%20clean%20angle%20will%20still%20work.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(aacc%2C%20disp_noisy%2C%20fc_sweep%2C%20freq%2C%20np%2C%20optcutfreq%2C%20rmse_sweep)%3A%0A%20%20%20%20_fc_opt%20%3D%20optcutfreq(disp_noisy%2C%20freq%3Dfreq%2C%20show%3DTrue)%0A%20%20%20%20_rmse%20%3D%20rmse_sweep(disp_noisy)%0A%20%20%20%20print(f%22Optimal%20cutoff%20frequency%2C%20noisy%20angle%3A%20%7B_fc_opt%3A.2f%7D%20Hz%22)%0A%20%20%20%20print(f%22Best%20cutoff%20for%20the%20acceleration%3A%20%20%20%20%20%20%7Bfc_sweep%5Bnp.argmin(_rmse)%5D%3A.2f%7D%20Hz%22)%0A%20%20%20%20for%20_fc%20in%20(_fc_opt%2C%20fc_sweep%5Bnp.argmin(_rmse)%5D%2C%209)%3A%0A%20%20%20%20%20%20%20%20print(f%22Cutoff%20%7B_fc%3A5.2f%7D%20Hz%3A%20RMSE%20%3D%20%7Bnp.sqrt(np.mean((filtered_acceleration(disp_noisy%2C%20freq%2C%20_fc)%20-%20aacc%5B1%3A-1%5D)%20**%202))%3A.2f%7D%20rad%2Fs2%22)%0A%20%20%20%20print(f%22No%20filter%20at%20all%3A%20%20%20%20%20RMSE%20%3D%20%7Bnp.sqrt(np.mean((np.diff(disp_noisy%2C%202)%20*%20freq**2%20-%20aacc%5B1%3A-1%5D)%20**%202))%3A.2f%7D%20rad%2Fs2%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%20With%20more%20noise%2C%20the%20residual%20analysis%20lowers%20its%20cutoff%20to%20about%203.9%20Hz%2C%20and%20its%20noise%20level%20%24a%24%20rises%20to%20about%200.007%20rad%2C%20four%20times%20higher.%20The%20best%20cutoff%20for%20the%20acceleration%20is%20now%205%20Hz%2C%20so%20the%20residual%20analysis%20oversmooths%20more%20than%20before%2C%20at%20a%20cost%20of%205.0%20instead%20of%204.5%20rad%2Fs%C2%B2.%20But%20the%209%20Hz%20that%20worked%20well%20on%20the%20clean%20angle%20now%20gives%20almost%20twice%20that%20error%2C%20and%20no%20filter%20at%20all%20gives%2037%20rad%2Fs%C2%B2.%0A%0A%20%20%20%20That%20is%20the%20real%20value%20of%20the%20method.%20It%20is%20not%20that%20it%20finds%20*the*%20optimal%20cutoff%3B%20it%20is%20that%20it%20adapts%20the%20cutoff%20to%20the%20noise%20in%20each%20recording%2C%20without%20guaranteeing%20an%20optimum%20for%20the%20quantity%20you%20ultimately%20analyze.%20For%20this%20benchmark%2C%20mild%20oversmoothing%20costs%20little%20near%20the%20RMSE%20minimum%3B%20rapid%20events%20or%20other%20noise%20spectra%20can%20behave%20differently.%0A%0A%20%20%20%20**Challenge%201.**%20The%20residual%20analysis%20assumes%20that%20the%20noise%20is%20white%20and%20that%20the%20signal%20has%20frequencies%20well%20below%20the%20Nyquist%20frequency.%0A%0A%20%20%20%201.%20Run%20%60optcutfreq(rng.standard_normal(100)%2C%20freq%3D100%2C%20show%3DTrue)%60%2C%20with%20%60rng%20%3D%20np.random.default_rng(seed%3D42)%60%3A%20a%20signal%20that%20is%20only%20noise.%20What%20happens%2C%20and%20why%3F%0A%20%20%20%202.%20Run%20it%20on%20a%20random%20walk%2C%20%60np.cumsum(rng.standard_normal(1000))%60%2C%20sampled%20at%201000%20Hz.%20Is%20the%20result%20plausible%3F%20Look%20at%20the%20residual%20plot%20and%20choose%20the%20limits%20of%20the%20noisy%20region%20yourself%20with%20%60fclim%60.%0A%20%20%20%203.%20Add%20a%2060%20Hz%20sinusoid%20evaluated%20at%20the%20angle%E2%80%99s%20sample%20times%20and%20run%20the%20residual%20analysis%20again.%20First%20compute%20its%20alias%20at%20this%20approximately%2049.75%20Hz%20sampling%20rate.%20Does%20the%20residual%20plot%20diagnose%20a%20narrowband%20contaminant%20uniquely%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%204.**%20Answer%20Challenge%201%3A%20what%20happens%20for%20pure%20noise%2C%20a%20random%20walk%2C%20and%2060%20Hz%20interference%20sampled%20with%20Pezzack%E2%80%99s%20time%20vector%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%204%0A%20%20%20%20For%20the%20specified%20100-sample%20white-noise%20record%2C%20%60optcutfreq%60%20returns%20%60None%60%3A%20it%20cannot%20identify%20a%20suitable%20signal%2Fnoise%20transition.%20A%20random%20walk%20has%20strongly%20colored%20low-frequency%20content%3B%20obtaining%20a%20numerical%20crossing%20or%20choosing%20%60fclim%60%20manually%20does%20not%20establish%20a%20physically%20meaningful%20signal%20bandwidth.%20At%20%24f_s%5Capprox49.7512%24%20Hz%2C%2060%20Hz%20aliases%20to%20%24%7C60-f_s%7C%5Capprox10.2488%24%20Hz.%20A%20residual%20curve%20may%20change%2C%20but%20cannot%20label%20this%20uniquely%20as%20interference%20or%20undo%20aliasing.%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%20answer%20to%20Challenge%200.%20Which%20parts%20of%20it%20are%20in%20Winter's%20method%2C%20and%20which%20are%20not%3F%0A%20%20%20%202.%20What%20does%20the%20residual%20measure%2C%20and%20why%20does%20it%20approach%20zero%20as%20the%20filter%20tends%20to%20the%20identity%3F%20(A%20digital%20Butterworth%20cutoff%20cannot%20be%20set%20exactly%20at%20Nyquist.)%0A%20%20%20%203.%20Why%20does%20the%20intercept%20of%20the%20straight%20part%20of%20the%20residual%20curve%20estimate%20the%20RMS%20of%20the%20noise%3F%20Which%20assumption%20about%20the%20noise%20does%20that%20require%3F%0A%20%20%20%204.%20The%20residual%20analysis%20chooses%20the%20cutoff%20from%20the%20angle%2C%20but%20you%20will%20use%20the%20acceleration.%20Why%20might%20the%20best%20cutoff%20for%20the%20two%20be%20different%3F%0A%20%20%20%205.%20Two%20recordings%20of%20the%20same%20movement%2C%20one%20noisier%20than%20the%20other%2C%20get%20different%20cutoff%20frequencies%20from%20the%20residual%20analysis.%20Is%20that%20a%20problem%20when%20you%20compare%20them%3F%20What%20would%20you%20report%3F%0A%20%20%20%206.%20A%20colleague%20uses%206%20Hz%20for%20every%20gait%20study%20%22because%20it%20is%20standard%22.%20Give%20one%20argument%20for%20and%20one%20against.%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.**%20Two%20trials%20have%20residual%20cutoffs%204%20and%206%20Hz.%20How%20would%20you%20compare%20their%20peak%20accelerations%20responsibly%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%205%0A%20%20%20%20First%20inspect%20the%20residual%20tails%20and%20document%20the%20filter%20family%2C%20pass%20count%2C%20cutoff%20convention%20and%20boundary%20handling.%20Since%20different%20responses%20can%20change%20peaks%20differently%2C%20repeat%20the%20comparison%20with%20defensible%20common%20cutoffs%20and%20assess%20sensitivity.%20Trial-specific%20filtering%20is%20not%20automatically%20invalid%2C%20but%20a%20peak%20difference%20that%20disappears%20under%20reasonable%20processing%20choices%20needs%20qualified%20interpretation.%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%2F%2F%2F%20details%20%7C%20Checkpoint%20answers%0A%20%20%20%201.%20Winter%E2%80%99s%20method%20supplies%20a%20residual%2Fnoise%20criterion%3B%20it%20does%20not%20supply%20a%20task-specific%20model%2C%20prove%20white%20noise%20or%20optimize%20every%20derivative.%0A%20%20%20%202.%20RMS%20of%20raw%20minus%20filtered%20data.%20It%20tends%20to%20zero%20when%20filtering%20changes%20nothing%3B%20the%20nominal%20grid%20endpoint%20in%20this%20implementation%20is%20about%20%240.802f_s%2F2%24%2C%20not%20%24f_s%2F2%24.%0A%20%20%20%203.%20The%20tail%20extrapolation%20is%20an%20empirical%20estimate%20under%20a%20low-bandwidth%20movement%2Fbroadband-noise%20separation.%20Correlated%20noise%20or%20a%20poor%20tail%20fit%20can%20invalidate%20it.%0A%20%20%20%204.%20Differentiation%20weights%20high%20frequencies%20more%20strongly%2C%20changing%20the%20relevant%20distortion%2Fnoise%20trade-off.%0A%20%20%20%205.%20Different%20noise%20levels%20can%20justify%20different%20cutoffs%2C%20but%20the%20resulting%20comparisons%20need%20sensitivity%20checks%20and%20a%20report%20of%20the%20per-trial%20choices%20and%20method.%0A%20%20%20%206.%20A%20fixed%20value%20makes%20processing%20consistent%20and%20may%20be%20validated%20for%20a%20particular%20protocol.%20It%20cannot%20automatically%20accommodate%20different%20noise%20levels%2C%20movement%20bandwidths%20or%20outcomes.%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.%20Generate%20a%20known%20signal%2C%20such%20as%20%24x(t)%20%3D%20%5Csin(2%5Cpi%20t)%20%2B%200.3%5Csin(6%5Cpi%20t)%24%20sampled%20at%20100%20Hz%20for%202%20s%2C%20and%20add%20white%20noise%20with%20standard%20deviations%20of%200.01%20and%200.05%20(use%20%60t%3Dnp.arange(200)%2F100%60%20and%20reset%20%60rng%3Dnp.random.default_rng(42)%60%20for%20each%20noise%20level).%20Use%20the%20historical%20approximate%20correction%20above.%20For%20each%3A%3Cbr%3E%0A%20%20%20%20%20%20%20a.%20Find%20the%20optimal%20cutoff%20frequency%20with%20%60optcutfreq%60%2C%20and%20compare%20the%20noise%20level%20%24a%24%20it%20estimates%20with%20the%20true%20standard%20deviation%20of%20the%20noise.%3Cbr%3E%0A%20%20%20%20%20%20%20b.%20Find%20the%20cutoff%20that%20minimizes%20the%20RMSE%20of%20the%20filtered%20signal%20relative%20to%20the%20noise-free%20one%2C%20and%20the%20one%20that%20minimizes%20the%20RMSE%20of%20its%20second%20derivative.%20Use%20cutoffs%20from%201%20to%2039.5%20Hz%20in%200.1%20Hz%20steps%2C%20all%20samples%20for%20signal%20RMSE%2C%20and%20the%20interior%20samples%20for%20the%20central%20second%20difference%20compared%20with%20the%20analytic%20second%20derivative.%20How%20do%20they%20compare%20with%20the%20residual%20analysis%3F%0A%0A%20%20%20%202.%20Repeat%20the%20analysis%20of%20this%20notebook%20with%20a%20critically%20damped%20filter%20instead%20of%20a%20Butterworth%20filter%20(see%20%60critic_damp%60%20in%20%5BData%20filtering%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fdata_filtering.py)).%20To%20isolate%20filter%20shape%2C%20use%20the%20exact%20digital%20correction%20for%20both%20filter%20families%2C%20the%20101-point%20frequency%20grid%20from%20%60residuals%60%2C%20and%20fit%20the%20same%20tail%20interval%20%248%5Cle%20f_c%3C18%24%20Hz.%20Does%20the%20residual-analysis%20cutoff%20frequency%20change%3F%0A%0A%20%20%20%203.%20Apply%20the%20residual%20analysis%20to%20each%20coordinate%20of%20the%20markers%20in%20Winter's%20Table%20A.1%2C%20used%20in%20%5BAngular%20kinematics%20in%20a%20plane%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fkinematics_angular_2d.py).%20Are%20the%20optimal%20cutoff%20frequencies%20similar%20for%20all%20markers%20and%20coordinates%3F%20Which%20would%20you%20use%2C%20and%20why%3F%0A%0A%20%20%20%204.%20Write%20a%20paragraph%20assessing%20residual%20analysis%20using%20this%20notebook%E2%80%99s%20assumptions%20and%20benchmark%20results%3A%20discuss%20tail%20selection%2C%20colored%20noise%2C%20differentiation%20and%20endpoint%20effects.%20Does%20the%20evidence%20support%20using%20the%20returned%20value%20without%20inspection%3F%20An%20optional%20historical%20gait-analysis%20discussion%20is%20linked%20under%20Go%20deeper%3B%20this%20problem%20does%20not%20depend%20on%20its%20availability.%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.%20With%20the%20specified%20seed%2C%20grids%20and%20historical%20%60optcutfreq%60%20convention%3A%0A%0A%20%20%20%20%7C%20Noise%20SD%20%7C%20Residual%20cutoff%20%7C%20Estimated%20%24a%24%20%7C%20Minimum%20signal-RMSE%20cutoff%20%7C%20Minimum%20acceleration-RMSE%20cutoff%20%7C%0A%20%20%20%20%7C---%7C---%3A%7C---%3A%7C---%3A%7C---%3A%7C%0A%20%20%20%20%7C%200.01%20%7C%206.2781%20Hz%20%7C%200.008895%20%7C%208.1%20Hz%20%7C%205.7%20Hz%20%7C%0A%20%20%20%20%7C%200.05%20%7C%204.0501%20Hz%20%7C%200.047831%20%7C%205.8%20Hz%20%7C%204.3%20Hz%20%7C%0A%0A%20%20%20%20The%20minimum%20signal%20RMSEs%20are%20approximately%200.003882%20and%200.014889%3B%20the%20minimum%20second-derivative%20RMSEs%20are%204.8043%20and%2012.5421%20in%20the%20signal%E2%80%99s%20units%20per%20second%20squared.%20The%20residual%20estimate%20is%20a%20distinct%20criterion%3B%20none%20of%20these%20cutoffs%20is%20universally%20optimal.%0A%0A%20%20%20%202.%20With%20the%20stated%20exact%20digital%20correction%20and%20common%208%E2%80%9318%20Hz%20tail%20fit%2C%20the%20cleaner%20Pezzack%20angle%20gives%20approximately%205.9842%20Hz%20and%20%24a%3D0.0016992%24%20rad%20for%20Butterworth%2C%20versus%207.0070%20Hz%20and%20%24a%3D0.0040652%24%20rad%20for%20the%20critically%20damped%20filter.%20The%20cutoff%20and%20inferred%20noise%20estimate%20depend%20on%20filter%20shape%20and%20tail%20choice.%20These%20values%20intentionally%20differ%20from%20the%20package%E2%80%99s%20automatic%20historical-convention%20result.%0A%0A%20%20%20%203.%20At%20the%20sampling%20rate%20inferred%20from%20the%20timestamps%2C%20about%2069.9534%20Hz%2C%20%60optcutfreq%60%20gives%20the%20following%20nominal%20cutoffs%20in%20Hz%3A%0A%0A%20%20%20%20%7C%20Marker%20%7C%20x%20%7C%20y%20%7C%0A%20%20%20%20%7C---%7C---%3A%7C---%3A%7C%0A%20%20%20%20%7C%20Rib%20cage%20%7C%203.64%20%7C%202.97%20%7C%0A%20%20%20%20%7C%20Hip%20%7C%205.17%20%7C%202.87%20%7C%0A%20%20%20%20%7C%20Knee%20%7C%206.28%20%7C%203.43%20%7C%0A%20%20%20%20%7C%20Fibula%20%7C%205.91%20%7C%203.08%20%7C%0A%20%20%20%20%7C%20Ankle%20%7C%206.55%20%7C%204.28%20%7C%0A%20%20%20%20%7C%20Heel%20%7C%205.56%20%7C%205.35%20%7C%0A%20%20%20%20%7C%20MT5%20%7C%205.71%20%7C%205.64%20%7C%0A%20%20%20%20%7C%20Toe%20%7C%207.28%20%7C%205.44%20%7C%0A%0A%20%20%20%20The%20values%20are%20not%20uniform.%20There%20is%20no%20unique%20justified%20common%20cutoff%20from%20this%20table%20alone.%20A%20defensible%20answer%20inspects%20the%20residual%20fits%2C%20identifies%20the%20movement%20features%20to%20retain%2C%20and%20reports%20the%20sensitivity%20of%20the%20intended%20outcome%20to%20common%20versus%20coordinate-specific%20filtering.%0A%0A%20%20%20%204.%20A%20suitable%20conclusion%20is%3A%20residual%20analysis%20is%20a%20useful%20diagnostic%2C%20but%20the%20tail%20may%20not%20isolate%20noise%3B%20colored%20noise%2C%20endpoint%20behavior%20and%20differentiation%20can%20alter%20the%20appropriate%20cutoff.%20The%20Pezzack%20examples%20show%20modest%20disagreement%20with%20an%20acceleration-based%20optimum%20and%20increased%20disagreement%20after%20adding%20noise.%20They%20support%20inspecting%20the%20curve%20and%20testing%20sensitivity%2C%20not%20accepting%20every%20returned%20frequency%20automatically.%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%5BData%20filtering%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fdata_filtering.py)%20%E2%80%94%20the%20filters%20used%20here%2C%20and%20other%20ways%20of%20dealing%20with%20noise.%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%20noise%2C%20the%20signal-to-noise%20ratio%20and%20the%20Nyquist%20frequency.%0A%0A%20%20%20%20The%20%5Bhistorical%20clinical%20gait-analysis%20discussion%5D(https%3A%2F%2Fwww.clinicalgaitanalysis.com%2Ffaq%2Fcutoff.html)%20is%20optional%20background%3B%20the%20self-contained%20critique%20in%20Problem%204%20uses%20the%20evidence%20in%20this%20guide.%0A%0A%20%20%20%20To%20read%20more%20about%20the%20determination%20of%20the%20optimal%20cutoff%20frequency%2C%20see%20the%20following%20papers%3A%0A%0A%20%20%20%20-%20Pezzack%20JC%2C%20Norman%20RW%2C%20Winter%20DA%20(1977)%20An%20assessment%20of%20derivative%20determining%20techniques%20used%20for%20motion%20analysis.%20Journal%20of%20Biomechanics%2C%2010%2C%20377-382.%0A%20%20%20%20-%20Giakas%20G%2C%20Baltzopoulos%20V%20(1997)%20A%20comparison%20of%20automatic%20filtering%20techniques%20applied%20to%20biomechanical%20walking%20data.%20Journal%20of%20Biomechanics%2C%2030%2C%20847-850.%0A%20%20%20%20-%20Alonso%20FJ%2C%20Salgado%20DR%2C%20Cuadrado%20J%2C%20Pintado%20P%20(2009)%20%5BAutomatic%20smoothing%20of%20raw%20kinematic%20signals%20using%20SSA%20and%20cluster%20analysis%5D(https%3A%2F%2Flim.ii.udc.es%2Fdocs%2Fproceedings%2F2009_09_EUROMECH_Automatic.pdf).%207th%20EUROMECH%20Solid%20Mechanics%20Conference.%0A%20%20%20%20-%20Kristianslund%20E%2C%20Krosshaug%20T%2C%20van%20den%20Bogert%20AJ%20(2012)%20Effect%20of%20low%20pass%20filtering%20on%20joint%20moments%20from%20inverse%20dynamics%3A%20Implications%20for%20injury%20prevention.%20Journal%20of%20Biomechanics%2C%2045%2C%20666-671.%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).%20Use%20the%20frequency-response%20explanation%20to%20understand%20why%20a%20residual%20curve%20depends%20on%20the%20chosen%20filter%20family.%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-%20Pezzack%20JC%2C%20Norman%20RW%2C%20Winter%20DA%20(1977)%20An%20assessment%20of%20derivative%20determining%20techniques%20used%20for%20motion%20analysis.%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-%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%2C%20with%20its%20noisier%20version%20of%20the%20angle.)%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