electromyography

%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%22detecta%3E%3D0.0.5%22%2C%20%22psd2%3E%3D0.0.5%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()%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%20Introduction%20to%20data%20analysis%20in%20electromyography%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%20Follow%20the%20recorded%20signal%20through%20filtering%2C%20rectification%2C%20envelope%20extraction%2C%20onset%20detection%2C%20amplitude%20and%20frequency%20summaries.%20Keep%20each%20operation%E2%80%99s%20units%20and%20time%20scale%20in%20view.%20An%20EMG%20amplitude%20is%20an%20electrical%20measurement%20whose%20relationship%20to%20force%20depends%20on%20the%20task%20and%20measurement%20setup.%0A%0A%20%20%20%20Try%20the%20distributed%20**Exercises**%20before%20opening%20their%20collapsed%20worked%20solutions.%20Then%20use%20the%20**Checkpoint%20questions**%20and%20attempt%20the%20**Problems**%20independently%3B%20their%20panels%20contain%20final%20answers%20only.%20Predict%20the%20effect%20of%20a%20parameter%20change%20before%20running%20it%2C%20and%20check%20units%20and%20limiting%20cases.%0A%0A%20%20%20%20**Opening%20challenge.**%20An%20EMG%20trace%20oscillates%20around%20zero%20during%20a%20strong%20contraction.%20Does%20a%20near-zero%20average%20mean%20that%20the%20muscle%20is%20inactive%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%20Python%20setup%0A%0A%20%20%20%20The%20dependency%20declaration%20at%20the%20top%20of%20this%20file%20lists%20the%20packages%20used%20by%20this%20guide.%20Marimo%20runs%20cells%20in%20dependency%20order.%20Run%20the%20notebook%20in%20an%20editor%20to%20use%20its%20displays%20and%20controls%3B%20the%20computational%20examples%20below%20explain%20their%20inputs%20before%20using%20them.%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%20Electromyography%20(EMG)%2C%20from%20Latin%2C%20literally%20means%20recording%20of%20muscle%20electricity.%20Nowadays%2C%20a%20more%20accurate%20definition%20is%20that%20EMG%20is%20the%20technique%20employed%20to%20register%20the%20electrical%20signal%20related%20to%20muscle%20activation.%0A%0A%20%20%20%20This%20electrical%20signal%20represents%20the%20electrical%20current%20generated%20by%20the%20ionic%20flow%20across%20the%20membranes%20of%20muscle%20fibers%20when%20the%20motor%20unit%20is%20fired.%20A%20motor%20unit%20is%20the%20functional%20unit%20of%20a%20muscle%20and%20consists%20of%20an%20alpha%20motor%20neuron%20and%20the%20muscle%20fibers%20innervated%20by%20this%20motor%20neuron.%20This%20electrical%20signal%20can%20be%20captured%20by%20invasive%20needle%20or%20wire%20electrodes%20placed%20directly%20near%20the%20muscle%20fibers%20(intramuscular%20EMG)%20of%20by%20electrodes%20placed%20on%20the%20skin%20where%20it%20is%20captured%20the%20electrical%20signal%20that%20propagates%20through%20the%20tissues%20of%20the%20human%20body%20(surface%20EMG).%20The%20use%20of%20surface%20EMG%20is%20much%20more%20common%20in%20human%20motion%20studies%20and%20here%20will%20see%20an%20introduction%20to%20the%20analysis%20of%20the%20surface%20EMG%20signal.%0A%0A%20%20%20%20For%20a%20typical%20muscle%20activation%20during%20voluntary%20human%20movement%2C%20the%20measured%20surface%20EMG%20signal%20is%20the%20result%20of%20an%20algebraic%20summation%20of%20many%20motor%20unit%20action%20potentials.%20A%20typical%20surface%20EMG%20signal%20has%20a%20random%20characteristic%20with%20zero%20mean%2C%20about%201%20mV%20of%20peak-to-peak%20amplitude%2C%20and%20frequencies%20between%205-10%20and%20400-450%20Hz%2C%20but%20these%20values%20can%20be%20affected%20by%20many%20factors%20including%20intensity%20of%20the%20muscle%20activation%2C%20the%20distance%20from%20the%20motor%20units%20to%20the%20skin%20surface%2C%20the%20geometry%20of%20the%20electrodes%2C%20etc.%0A%0A%20%20%20%20Let's%20load%20a%20surface%20EMG%20signal%20from%20a%20voluntary%20contraction%20and%20analyze%20it.%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%20numpy%20as%20np%0A%20%20%20%20%23%20'%25matplotlib%20inline'%20command%20supported%20automatically%20in%20marimo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%0A%20%20%20%20return%20np%2C%20plt%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20repo_file)%3A%0A%20%20%20%20import%20pandas%20as%20pd%0A%0A%20%20%20%20%23%20load%20data%20file%0A%20%20%20%20time%2C%20data%20%3D%20pd.read_csv(repo_file('data%2Femg.csv')%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20header%3DNone).to_numpy().T%0A%20%20%20%20freq%20%3D%201%2Fnp.mean(np.diff(time))%0A%20%20%20%20return%20data%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%20And%20let's%20plot%20the%20data%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20plt%2C%20time)%3A%0A%20%20%20%20(_fig%2C%20_ax)%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(8%2C%203))%0A%20%20%20%20_ax.plot(time%2C%20data%2C%20'b')%0A%20%20%20%20_ax.set_xlabel('Time%20%5Bs%5D'%2C%20fontsize%3D12)%0A%20%20%20%20_ax.set_ylabel('EMG%20amplitude%20%5BV%5D'%2C%20fontsize%3D12)%0A%20%20%20%20_ax.set_xlim(time%5B0%5D%2C%20time%5B-1%5D)%0A%20%20%20%20_ax.set_title('Raw%20data'%2C%20fontsize%3D16)%0A%20%20%20%20plt.locator_params(axis%3D'both'%2C%20nbins%3D4)%0A%20%20%20%20plt.grid()%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%23%23%20Filtering%0A%0A%20%20%20%20Any%20bioelectrical%20signal%20is%20contaminated%20by%20noise%2C%20be%20it%20from%20other%20electrical%20sources%20in%20the%20human%20body%2C%20from%20external%20sources%2C%20or%20by%20the%20own%20process%20of%20measurement.%0A%20%20%20%20To%20get%20rid%20of%20part%20of%20this%20noise%2C%20we%20can%20apply%20a%20band-pass%20filter%20to%20only%20pass%20signals%20with%20frequencies%20in%20the%20desired%20range.%0A%20%20%20%20A%20common%20choice%20for%20filtering%20EMG%20data%20it's%20the%20Butterworth%20filter%20with%20zero%20lag.%20For%20more%20about%20data%20filtering%20%5Bclick%20here%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fdata_filtering.py).%3Cbr%3E%0A%20%20%20%20Let's%20use%20%60butter(N%3D2%2C%20...%2C%20output%3D'sos'%2C%20fs%3Dfreq)%60%20at%2010%E2%80%93400%20Hz%20and%20apply%20its%20second-order%20sections%20forward%20and%20backward.%20A%20band-pass%20transformation%20produces%20a%20fourth-order%20one-pass%20filter%3B%20forward-backward%20filtering%20doubles%20the%20order%20and%20squares%20its%20magnitude%20response.%20The%20specified%20cutoffs%20therefore%20have%20approximately%20%E2%88%926%20dB%20gain%20after%20both%20passes.%20Zero%20phase%20requires%20future%20samples%3A%20this%20is%20offline%20processing%2C%20and%20edge%20transients%20still%20need%20inspection.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20freq)%3A%0A%20%20%20%20from%20scipy.signal%20import%20butter%2C%20sosfiltfilt%0A%20%20%20%20_sos%20%3D%20butter(2%2C%20%5B10%2C%20400%5D%2C%20btype%3D'bandpass'%2C%20fs%3Dfreq%2C%20output%3D'sos')%0A%20%20%20%20dataf%20%3D%20sosfiltfilt(_sos%2C%20data)%0A%20%20%20%20return%20butter%2C%20dataf%2C%20sosfiltfilt%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%20Note%20that%20the%20high-pass%20filter%20removes%20any%20offset%20in%20the%20signal%20(a%20constant%20value%20or%20DC%20component%20can%20be%20thought%20as%20a%20signal%20with%20frequency%20equals%20to%20zero).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20dataf%2C%20plt%2C%20time)%3A%0A%20%20%20%20(_fig%2C%20(_ax1%2C%20_ax2))%20%3D%20plt.subplots(2%2C%201%2C%20sharex%3DTrue%2C%20sharey%3DTrue%2C%20figsize%3D(8%2C%205))%0A%20%20%20%20_ax1.plot(time%2C%20data%2C%20'r')%0A%20%20%20%20_ax1.set_title('Raw%20data'%2C%20fontsize%3D16)%0A%20%20%20%20_ax1.grid()%0A%20%20%20%20_ax2.plot(time%2C%20dataf%2C%20'b')%0A%20%20%20%20_ax2.set_xlabel('Time%20%5Bs%5D'%2C%20fontsize%3D12)%0A%20%20%20%20_ax2.set_xlim(time%5B0%5D%2C%20time%5B-1%5D)%0A%20%20%20%20_ax2.set_title('Band-pass%20filter%20at%20%5B10%20400%5D%20Hz'%2C%20fontsize%3D16)%0A%20%20%20%20_fig.text(0%2C%200.5%2C%20'EMG%20amplitude%20%5BV%5D'%2C%20ha%3D'center'%2C%20va%3D'center'%2C%20rotation%3D'vertical'%2C%20fontsize%3D12)%0A%20%20%20%20plt.locator_params(axis%3D'both'%2C%20nbins%3D4)%0A%20%20%20%20_ax2.grid()%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%201%20%E2%80%94%20Check%20a%20filter%20against%20sampling.**%20Can%20a%2010%E2%80%93400%20Hz%20band-pass%20be%20applied%20to%20data%20sampled%20at%20500%20Hz%3F%20What%20must%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%20Solution%20to%20Exercise%201%0A%0A%20%20%20%20No%3A%20the%20Nyquist%20frequency%20is%20250%20Hz%2C%20below%20the%20desired%20upper%20cutoff.%20Use%20data%20acquired%20at%20a%20sufficiently%20high%20sampling%20rate%20with%20appropriate%20anti-alias%20filtering%2C%20or%20choose%20a%20valid%20band%20within%20the%20acquired%20bandwidth.%20Upsampling%20cannot%20recover%20already%20aliased%20content.%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%20Linear%20envelope%0A%0A%20%20%20%20A%20common%20processing%20of%20the%20EMG%20signal%20is%20to%20calculate%20an%20activation%20level%20of%20the%20signal%2C%20a%20process%20called%20linear%20envelope.%20The%20linear%20envelope%20consists%20in%20two%20steps%3A%20the%20signal%20is%20full-wave%20rectified%20by%20computing%20the%20absolute%20value%20of%20the%20signal%3B%20then%20the%20rectified%20signal%20is%20low-pass%20filtered%20with%20a%20cutoff%20frequency%20typically%20at%20the%20range%20of%203%20to%208%20Hz%2C%20depending%20of%20the%20contraction%20muscle%20characteristics%20and%20the%20specific%20application%20for%20the%20linear%20envelope%20processing.%20A%20related%20smoothing%20estimate%20can%20be%20obtained%20with%20a%20moving%20average%20of%20the%20rectified%20signal%20with%20a%20moving%20window%20of%20100%20to%20200%20ms%20of%20duration.%20%20%20Alternatively%2C%20instead%20of%20the%20two%20previous%20steps%2C%20it%20can%20be%20applied%20a%20moving%20root%20mean%20square%20(RMS)%20with%20a%20moving%20window%20of%20100%20to%20200%20ms%20of%20duration%20to%20the%20original%20signal.%20These%20methods%20have%20different%20frequency%20responses%20and%20amplitude%20statistics%2C%20so%20their%20window%20lengths%20and%20values%20are%20not%20interchangeable.%0A%0A%20%20%20%20Let's%20calculate%20the%20linear%20envelope.%20First%2C%20the%20full-wave%20rectification%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dataf)%3A%0A%20%20%20%20%23%20full-wave%20rectification%0A%20%20%20%20datafr%20%3D%20abs(dataf)%0A%20%20%20%20return%20(datafr%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%20Second%2C%20a%20low-pass%20filter%20at%208%20Hz%20using%20a%20Butterworth%20filter%20of%20second%20order%20and%20zero%20lag%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(butter%2C%20datafr%2C%20freq%2C%20sosfiltfilt)%3A%0A%20%20%20%20_sos%20%3D%20butter(2%2C%208%2C%20btype%3D'lowpass'%2C%20fs%3Dfreq%2C%20output%3D'sos')%0A%20%20%20%20datafrle%20%3D%20sosfiltfilt(_sos%2C%20datafr)%0A%20%20%20%20return%20(datafrle%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%20Let's%20plot%20the%20results%20of%20each%20operation%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20dataf%2C%20datafr%2C%20datafrle%2C%20plt%2C%20time)%3A%0A%20%20%20%20(_fig%2C%20(_ax1%2C%20_ax2%2C%20ax3%2C%20ax4))%20%3D%20plt.subplots(4%2C%201%2C%20sharex%3DTrue%2C%20sharey%3DTrue%2C%20figsize%3D(8%2C%207))%0A%20%20%20%20_ax1.plot(time%2C%20data%2C%20'r')%0A%20%20%20%20_ax1.set_title('Raw%20data'%2C%20fontsize%3D12)%0A%20%20%20%20_ax1.set_xlim(time%5B0%5D%2C%20time%5B-1%5D)%0A%20%20%20%20_ax1.grid()%0A%20%20%20%20_ax2.plot(time%2C%20dataf%2C%20'b')%0A%20%20%20%20_ax2.set_title('Band-pass%20filter%20at%20%5B10%20400%5D%20Hz'%2C%20fontsize%3D12)%0A%20%20%20%20_ax2.grid()%0A%20%20%20%20ax3.plot(time%2C%20datafr%2C%20'b')%0A%20%20%20%20ax3.set_title('Band-pass%20filter%20at%20%5B10%20400%5D%20Hz%20and%20rectification'%2C%20fontsize%3D12)%0A%20%20%20%20ax3.grid()%0A%20%20%20%20ax4.plot(time%2C%20datafrle%2C%20'b'%2C%20linewidth%3D2)%0A%20%20%20%20ax4.set_xlabel('Time%20%5Bs%5D'%2C%20fontsize%3D12)%0A%20%20%20%20ax4.set_title('Band-pass%20filter%20at%20%5B10%20400%5D%20Hz%2C%20rectification%20and%20low-pass%20filter%20at%208%20Hz'%2C%20fontsize%3D12)%0A%20%20%20%20_fig.text(0%2C%200.5%2C%20'EMG%20amplitude%20%5BV%5D'%2C%20ha%3D'center'%2C%20va%3D'center'%2C%20rotation%3D'vertical'%2C%20fontsize%3D12)%0A%20%20%20%20plt.locator_params(axis%3D'both'%2C%20nbins%3D4)%0A%20%20%20%20ax4.grid()%0A%20%20%20%20plt.tight_layout(h_pad%3D0.1)%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20datafrle%2C%20plt%2C%20time)%3A%0A%20%20%20%20(_fig%2C%20_ax)%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(8%2C%203))%0A%20%20%20%20_ax.plot(time%2C%20data%2C%20'r')%0A%20%20%20%20_ax.plot(time%2C%20datafrle%2C%20'b'%2C%20linewidth%3D2)%0A%20%20%20%20_ax.set_xlabel('Time%20%5Bs%5D'%2C%20fontsize%3D12)%0A%20%20%20%20_ax.set_ylabel('EMG%20amplitude%20%5BV%5D'%2C%20fontsize%3D12)%0A%20%20%20%20_ax.set_title('Band-pass%20filter%20at%20%5B10%20400%5D%20Hz%2C%20rectification%20and%20low-pass%20filter%20at%208%20Hz'%2C%20fontsize%3D12)%0A%20%20%20%20_ax.set_xlim(time%5B0%5D%2C%20time%5B-1%5D)%0A%20%20%20%20plt.locator_params(axis%3D'both'%2C%20nbins%3D4)%0A%20%20%20%20plt.grid()%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%20Note%20that%20to%20calculate%20the%20linear%20envelope%20of%20the%20EMG%20signal%2C%20we%20(1)%20removed%20all%20frequencies%20below%2010%20Hz%20and%20above%20400%20Hz%3B%20(2)%20rectified%20the%20signal%3B%20and%20(3)%20removed%20all%20frequencies%20above%208%20Hz.%0A%0A%20%20%20%20How%20it%20is%20possible%20to%20have%20any%20signal%20with%20frequencies%20below%208%20Hz%20after%20the%20third%20step%20if%20we%20removed%20all%20frequencies%20below%2010%20Hz%20in%20the%20first%20step%3F%20The%20trick%20is%20that%20the%20rectification%20operation%20creates%20new%20frequencies%20in%20the%20signal.%0A%0A%20%20%20%20The%20band-pass%20filtering%20in%20the%20first%20step%20is%20not%20part%20of%20the%20linear%20envelope%20processing.%20One%20can%20calculate%20the%20linear%20envelope%20directly%20for%20the%20raw%20signal.%20However%20if%20we%20do%20that%20for%20an%20raw%20EMG%20signal%20with%20an%20offset%2C%20its%20linear%20envelope%20will%20always%20have%20values%20much%20greater%20than%20zero%20suggesting%20a%20muscle%20activation%20where%20there%20was%20not.%0A%0A%20%20%20%20Let's%20write%20a%20function%20that%20calculates%20the%20linear%20envelope%20of%20a%20signal%20and%20allows%20the%20user%20to%20input%20different%20parameters%20for%20the%20filtering%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(butter%2C%20np%2C%20sosfiltfilt)%3A%0A%20%20%20%20def%20linear_envelope(x%2C%20freq%3D1000%2C%20fc_bp%3D%5B10%2C%20400%5D%2C%20fc_lp%3D8)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Calculate%20the%20linear%20envelope%20of%20a%20signal.%0A%0A%20%20%20%20%20%20%20%20Parameters%0A%20%20%20%20%20%20%20%20----------%0A%20%20%20%20%20%20%20%20x%20%20%20%20%20%3A%201D%20array_like%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20raw%20signal%0A%20%20%20%20%20%20%20%20freq%20%20%3A%20number%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20sampling%20frequency%0A%20%20%20%20%20%20%20%20fc_bp%20%3A%20list%20%5Bfc_h%2C%20fc_l%5D%2C%20optional%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20cutoff%20frequencies%20for%20the%20band-pass%20filter%20(in%20Hz)%0A%20%20%20%20%20%20%20%20fc_lp%20%3A%20number%2C%20optional%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20cutoff%20frequency%20for%20the%20low-pass%20filter%20(in%20Hz)%0A%0A%20%20%20%20%20%20%20%20Returns%0A%20%20%20%20%20%20%20%20-------%0A%20%20%20%20%20%20%20%20x%20%20%20%20%20%3A%201D%20array_like%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20linear%20envelope%20of%20the%20signal%0A%0A%20%20%20%20%20%20%20%20Notes%0A%20%20%20%20%20%20%20%20-----%0A%20%20%20%20%20%20%20%20Butterworth%20filters%20use%20N%3D2%20prototypes%20and%20SOS%20forward-backward%20filtering.%0A%20%20%20%20%20%20%20%20Effective%20orders%20are%20eight%20for%20band-pass%20and%20four%20for%20low-pass%20filtering.%0A%20%20%20%20%20%20%20%20This%20zero-phase%20operation%20is%20noncausal%20and%20intended%20for%20offline%20data.%0A%0A%20%20%20%20%20%20%20%20See%20this%20notebook%20%5B1%5D_.%0A%0A%20%20%20%20%20%20%20%20References%0A%20%20%20%20%20%20%20%20----------%0A%20%20%20%20%20%20%20%20..%20%5B1%5D%20electromyography.py%0A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20if%20np.size(fc_bp)%20%3D%3D%202%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_sos%20%3D%20butter(2%2C%20fc_bp%2C%20btype%3D'bandpass'%2C%20fs%3Dfreq%2C%20output%3D'sos')%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20sosfiltfilt(_sos%2C%20x)%0A%20%20%20%20%20%20%20%20if%20np.size(fc_lp)%20%3D%3D%201%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20abs(x)%0A%20%20%20%20%20%20%20%20%20%20%20%20_sos%20%3D%20butter(2%2C%20fc_lp%2C%20btype%3D'lowpass'%2C%20fs%3Dfreq%2C%20output%3D'sos')%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20sosfiltfilt(_sos%2C%20x)%0A%20%20%20%20%20%20%20%20return%20x%0A%0A%20%20%20%20return%20(linear_envelope%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**Exercise%202%20%E2%80%94%20Compare%20mean%2C%20rectification%20and%20RMS.**%20For%20samples%20%24(-2%2C2%2C-2%2C2)%24%20mV%2C%20find%20the%20arithmetic%20mean%2C%20mean%20absolute%20value%20and%20RMS.%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%20Solution%20to%20Exercise%202%0A%0A%20%20%20%20The%20signed%20values%20cancel%2C%20giving%20mean%200%20mV.%20Every%20absolute%20value%20is%202%20mV%2C%20so%20mean%20absolute%20value%20is%202%20mV.%20The%20mean%20square%20is%204%20mV%C2%B2%20and%20RMS%20is%202%20mV.%20This%20resolves%20the%20opening%20challenge.%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%20Onset%20detection%0A%0A%20%20%20%20A%20common%20operation%20in%20the%20analysis%20of%20EMG%20data%20is%20to%20select%20the%20period%20of%20muscle%20activation%20based%20on%20the%20amplitude%20of%20the%20EMG%20signal.%20Of%20course%20this%20can%20be%20done%20manually%20by%20showing%20the%20EMG%20signal%20on%20a%20window%20and%20the%20expert%20user%20manually%20selects%20the%20events%20with%20mouse%20clicks%2C%20but%20a%20lot%20of%20research%20has%20been%20conducted%20to%20propose%20automatic%20methods%20to%20perform%20this%20task.%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%20Amplitude-threshold%20method%0A%0A%20%20%20%20The%20simplest%20method%20to%20automatically%20detect%20the%20onset%20in%20the%20EMG%20signal%20is%20based%20on%20amplitude%20threshold%2C%20where%20the%20EMG%20signal%20is%20considered%20to%20be%20'on'%20when%20it%20is%20above%20a%20certain%20threshold.%20This%20threshold%20is%20usually%20proportional%20to%20the%20amplitude%20of%20the%20baseline%20(the%20part%20of%20the%20EMG%20signal%20that%20we%20know%20there%20is%20no%20muscle%20activation)%2C%20for%20example%2C%20a%20threshold%20equals%20to%20three%20times%20the%20standard%20deviation%20of%20the%20baseline.%20To%20avoid%20the%20detection%20of%20spikes%20in%20the%20EMG%20signal%2C%20a%20second%20threshold%20can%20be%20used%20imposing%20that%20the%20signal%20should%20be%20above%20the%20threshold%20for%20at%20least%20a%20certain%20period%20(number%20of%20samples).%0A%0A%20%20%20%20The%20raw%20EMG%20signal%20can%20not%20be%20analyzed%20for%20the%20onset%20detection%20based%20on%20amplitude%20threshold%20because%20its%20amplitude%20varies%20from%20positive%20to%20negative%20very%20fast.%20The%20linear%20envelope%20of%20the%20EMG%20signal%20is%20usually%20analyzed%20or%20the%20amplitude%20detection%20is%20performed%20in%20terms%20of%20the%20standard%20deviation%20or%20RMS%20of%20the%20signal.%0A%0A%20%20%20%20The%20maintained%20%5Bdetecta%20package%5D(https%3A%2F%2Fgithub.com%2Fdemotu%2Fdetecta)%20provides%20%60detect_onset%60.%20Let's%20apply%20its%20threshold%20and%20duration%20rules%20to%20the%20envelope%3A%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%20from%20detecta%20import%20detect_onset%0A%0A%20%20%20%20return%20(detect_onset%2C)%0A%0A%0A%40app.cell%0Adef%20_(data%2C%20detect_onset%2C%20freq%2C%20linear_envelope%2C%20np)%3A%0A%20%20%20%20_threshold%20%3D%202%0A%20%20%20%20_window%20%3D%2050%0A%20%20%20%20data2%20%3D%20linear_envelope(data%2C%20freq%2C%20fc_bp%3D%5B20%2C%20400%5D%2C%20fc_lp%3D20)%0A%20%20%20%20inds%20%3D%20detect_onset(data2%2C%20threshold%3D_threshold%20*%20np.std(data%5B0%3A200%5D)%2C%20n_above%3D_window%2C%20n_below%3D10%2C%20show%3DTrue)%0A%20%20%20%20inds%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%20The%20Teager%E2%80%93Kaiser%20Energy%20operator%20to%20improve%20onset%20detection%0A%0A%20%20%20%20The%20Teager%E2%80%93Kaiser%20Energy%20(TKE)%20operator%20has%20been%20proposed%20to%20increase%20the%20accuracy%20of%20the%20onset%20detection%20by%20improving%20the%20SNR%20of%20the%20EMG%20signal%20(Li%20et%20al.%2C%202007).%3Cbr%3E%0A%20%20%20%20The%20TKE%20operator%20processes%20the%20data%20in%20the%20following%20way%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20y%5Bn%5D%3Dx%5E2%5Bn%5D-x%5Bn-1%5D%5Ccdot%20x%5Bn%2B1%5D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24x%5Bn%5D%24%20is%20the%20sample%20%24n%24%20of%20the%20signal%20%24x%24.%0A%0A%20%20%20%20Let's%20write%20a%20function%20that%20implements%20the%20TKE%20operator%3A%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%20def%20tkeo(x)%3A%0A%20%20%20%20%20%20%20%20r%22%22%22Calculates%20the%20Teager%E2%80%93Kaiser%20Energy%20operator.%0A%0A%20%20%20%20%20%20%20%20Parameters%0A%20%20%20%20%20%20%20%20----------%0A%20%20%20%20%20%20%20%20x%20%3A%201D%20array_like%0A%20%20%20%20%20%20%20%20%20%20%20%20raw%20signal%0A%0A%20%20%20%20%20%20%20%20Returns%0A%20%20%20%20%20%20%20%20-------%0A%20%20%20%20%20%20%20%20y%20%3A%201D%20array_like%0A%20%20%20%20%20%20%20%20%20%20%20%20signal%20processed%20by%20the%20Teager%E2%80%93Kaiser%20Energy%20operator%0A%0A%20%20%20%20%20%20%20%20Notes%0A%20%20%20%20%20%20%20%20-----%0A%0A%20%20%20%20%20%20%20%20See%20this%20notebook%20%5B1%5D_.%0A%0A%20%20%20%20%20%20%20%20References%0A%20%20%20%20%20%20%20%20----------%0A%20%20%20%20%20%20%20%20..%20%5B1%5D%20electromyography.py%0A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20x%20%3D%20np.asarray(x)%0A%20%20%20%20%20%20%20%20y%20%3D%20np.copy(x)%0A%20%20%20%20%20%20%20%20%23%20Teager%E2%80%93Kaiser%20Energy%20operator%0A%20%20%20%20%20%20%20%20y%5B1%3A-1%5D%20%3D%20x%5B1%3A-1%5D*x%5B1%3A-1%5D%20-%20x%5B%3A-2%5D*x%5B2%3A%5D%0A%20%20%20%20%20%20%20%20%23%20correct%20the%20data%20in%20the%20extremities%0A%20%20%20%20%20%20%20%20y%5B0%5D%2C%20y%5B-1%5D%20%3D%20y%5B1%5D%2C%20y%5B-2%5D%0A%0A%20%20%20%20%20%20%20%20return%20y%0A%0A%20%20%20%20return%20(tkeo%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%20And%20let's%20test%20it%20(with%20the%20same%20parameters%20as%20before)%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dataf%2C%20detect_onset%2C%20freq%2C%20linear_envelope%2C%20np%2C%20tkeo)%3A%0A%20%20%20%20_threshold%20%3D%202%0A%20%20%20%20_window%20%3D%2050%0A%20%20%20%20data3%20%3D%20tkeo(dataf)%0A%20%20%20%20data4%20%3D%20linear_envelope(data3%2C%20freq%2C%20fc_bp%3D%5B20%2C%20400%5D%2C%20fc_lp%3D20)%0A%20%20%20%20inds2%20%3D%20detect_onset(data4%2C%20threshold%3D_threshold%20*%20np.std(data3%5B0%3A200%5D)%2C%20n_above%3D_window%2C%20n_below%3D10%2C%20show%3DTrue)%0A%20%20%20%20inds2%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%20%E2%80%94%20Interpret%20a%20duration%20threshold.**%20At%201000%20Hz%2C%20an%20onset%20detector%20requires%2050%20samples%20above%20threshold%20and%20tolerates%2010%20below.%20What%20time%20scales%20do%20these%20counts%20represent%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%20Solution%20to%20Exercise%203%0A%0A%20%20%20%20They%20correspond%20to%20nominal%20sample-count%20durations%20of%2050%20ms%20and%2010%20ms.%20The%20exact%20first-to-last%20timestamp%20span%20of%2050%20consecutive%20samples%20is%2049%20ms.%20Threshold%20crossings%20are%20algorithmic%20events%2C%20not%20direct%20neural-onset%20measurements.%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%20Amplitude%20estimation%0A%0A%20%20%20%20Different%20measures%20can%20be%20used%20to%20describe%20the%20amplitude%20of%20the%20EMG%20signal%3A%0A%0A%20%20%20%20*%20Mean%2C%20median%0A%20%20%20%20*%20Standard%20deviation%20(SD)%0A%20%20%20%20*%20Root%20mean%20square%20(RMS)%0A%20%20%20%20*%20Maximum%2C%20minimum%0A%20%20%20%20*%20Etc.%0A%0A%20%20%20%20From%20the%20above%2C%20RMS%20is%20probably%20the%20less%20usual%20and%20it%20defined%20as%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20RMS(x)%20%3D%20%5Csqrt%7B%5Cfrac%7B1%7D%7BN%7D%5Csum_%7Bi%3D1%7D%5E%7BN%7D%5C%3B%20x_i%5E2%20%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20And%20in%20Python%2FNumpy%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dataf%2C%20np)%3A%0A%20%20%20%20rms%20%3D%20np.sqrt(np.mean(dataf**2))%20%23%20root%20mean%20square%2C%20the%20name%20makes%20sense.%0A%20%20%20%20rms%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%20the%20measures%20above%20have%20the%20same%20unit%20as%20the%20original%20EMG%20signal%20and%20they%20can%20be%20applied%20directly%20to%20the%20raw%20or%20to%20the%20linear%20evelope%20of%20the%20EMG%20signal.%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%20Integration%20(area)%0A%0A%20%20%20%20The%20area%20of%20the%20signal-versus-time%20curve%20(integral)%20is%20another%20measure%20used%20in%20EMG%20processing.%3Cbr%3E%0A%20%20%20%20The%20definition%20of%20the%20integral%20for%20a%20continuous%20function%20is%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Ctext%7BArea%7D%3D%20%5Cint_%7Bt_i%7D%5E%7Bt_f%7Dx(t)dt%0A%20%20%20%20%24%24%0A%0A%20%20%20%20For%20discrete%20data%2C%20as%20the%20EMG%20signal%2C%20the%20integral%20can%20be%20calculated%20by%20numerical%20integration%20using%20the%20rectangle%20rule%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Ctext%7BArea%7D%20%5C%3B%5Capprox%5C%3B%20%5CDelta%20t%20%5Csum_%7Bi%3D1%7D%5E%7BN%7D%5C%3Ax(t_i)%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Or%20more%20accurately%2C%20using%20the%20trapezoidal%20rule%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Ctext%7BArea%7D%20%5C%3B%5Capprox%5C%3B%20%5CDelta%20t%20%5Csum_%7Bi%3D1%7D%5E%7BN-1%7D%5Cfrac%7Bx(t_i)%2Bx(t_%7Bi%2B1%7D)%7D%7B2%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20And%20in%20Python%2FNumpy%20(SciPy%20calls%20these%20functions%20%60trapezoid%60%20and%20%60cumulative_trapezoid%60%3B%20the%20old%20names%20%60trapz%60%20and%20%60cumtrapz%60%20were%20removed)%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dataf%2C%20freq)%3A%0A%20%20%20%20from%20scipy.integrate%20import%20cumulative_trapezoid%2C%20trapezoid%0A%0A%20%20%20%20x%20%3D%20dataf%0A%20%20%20%20Dt%20%3D%201%20%2F%20freq%20%20%23%20time%20increment%20between%20samples%0A%20%20%20%20_Arect%20%3D%20Dt%20*%20sum(x)%0A%20%20%20%20_Atrap%20%3D%20Dt%20*%20trapezoid(x)%0A%20%20%20%20print('Rectangle%20rule%3A%20%25g%20Vs%2C%20trapezoidal%20rule%3A%20%25g%20Vs'%20%25%20(_Arect%2C%20_Atrap))%0A%20%20%20%20return%20Dt%2C%20cumulative_trapezoid%2C%20trapezoid%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%20The%20integrals%20above%20are%20single%20numbers%2C%20the%20total%20area%20under%20the%20EMG%20curve.%3Cbr%3E%0A%20%20%20%20If%20we%20want%20to%20calculate%20the%20integral%20as%20a%20function%20of%20time%2C%20we%20have%20to%20use%20the%20cumulative%20integration%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Dt%2C%20cumulative_trapezoid%2C%20np%2C%20x)%3A%0A%20%20%20%20_Arect%20%3D%20Dt%20*%20np.cumsum(x)%0A%20%20%20%20_Atrap%20%3D%20Dt%20*%20cumulative_trapezoid(x)%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%20As%20the%20raw%20EMG%20signal%20has%20zero%20mean%2C%20one%20is%20usually%20interested%20in%20the%20area%20under%20the%20curve%20of%20the%20linear%20envelope.%3Cbr%3E%0A%20%20%20%20Let's%20experiment%20with%20these%20integration%20measures%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(datafrle%2C%20freq%2C%20trapezoid)%3A%0A%20%20%20%20_Arect%20%3D%20sum(datafrle)%20%2F%20freq%0A%20%20%20%20print('Total%20area%20by%20the%20rectangle%20rule%3A%20%25f%20Vs'%20%25%20_Arect)%0A%20%20%20%20_Atrap%20%3D%20trapezoid(datafrle)%20%2F%20freq%0A%20%20%20%20print('Total%20area%20by%20the%20trapezoid%20rule%3A%20%25f%20Vs'%20%25%20_Atrap)%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%20Note%20that%20the%20'Vs'%20above%20is%20not%20the%20plural%20of%20V%3B%20it's%20Volts%20times%20seconds.%3Cbr%3E%0A%20%20%20%20And%20now%20the%20cumulative%20integration%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cumulative_trapezoid%2C%20datafrle%2C%20freq%2C%20np)%3A%0A%20%20%20%20Arect2%20%3D%20np.cumsum(datafrle)%2Ffreq%0A%20%20%20%20Atrap2%20%3D%20cumulative_trapezoid(datafrle%2C%20initial%3D0)%2Ffreq%0A%20%20%20%20return%20Arect2%2C%20Atrap2%0A%0A%0A%40app.cell%0Adef%20_(Arect2%2C%20Atrap2%2C%20datafrle%2C%20plt%2C%20time)%3A%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(8%2C%205))%0A%20%20%20%20_ax1.plot(time%2C%20datafrle%2C%20'r')%0A%20%20%20%20_ax1.set_title('EMG%20signal%20(linear%20envelope)')%0A%20%20%20%20_ax1.set_ylabel('EMG%20amplitude%20%5BV%5D')%0A%20%20%20%20_ax1.set_xlim(time%5B0%5D%2C%20time%5B-1%5D)%0A%20%20%20%20_ax2.plot(time%2C%20Arect2%2C%20'b'%2C%20label%3D'Rectangle')%0A%20%20%20%20_ax2.plot(time%2C%20Atrap2%2C%20'g'%2C%20label%3D'Trapezoid')%0A%20%20%20%20_ax2.set_xlabel('Time%20%5Bs%5D')%0A%20%20%20%20_ax2.set_title('Integral%20of%20the%20EMG%20signal')%0A%20%20%20%20_ax2.set_ylabel('EMG%20integral%20%5BVs%5D')%0A%20%20%20%20_ax2.legend(loc%3D'upper%20left'%2C%20frameon%3DFalse)%0A%20%20%20%20plt.locator_params(axis%3D'both'%2C%20nbins%3D4)%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%20Note%20that%20the%20value%20of%20the%20total%20area%20is%20equal%20to%20the%20last%20value%20of%20the%20cumulative%20integration%20from%20the%20graph%20above.%0A%0A%20%20%20%20The%20trapezoidal%20rule%20is%20more%20accurate%20than%20the%20rectangle%20rule%20to%20calculate%20the%20integration%20of%20a%20function.%3Cbr%3E%0A%20%20%20%20So%2C%20why%20did%20they%20result%20in%20practically%20the%20same%20values%3F%3Cbr%3E%0A%20%20%20%20Writing%20out%20the%20two%20sums%2C%20the%20trapezoidal%20rule%20counts%20every%20sample%20once%2C%20except%20the%20first%20and%20the%20last%2C%20which%20it%20counts%20by%20half.%20The%20two%20totals%20differ%20only%20by%20%24%5CDelta%20t%5C%2C(x_1%20%2B%20x_N)%2F2%24%2C%20and%20the%20linear%20envelope%20is%20close%20to%20zero%20at%20the%20beginning%20and%20at%20the%20end%20of%20this%20recording.%20Along%20the%20way%2C%20the%20rectangle%20rule%20overestimates%20the%20area%20where%20the%20curve%20rises%20and%20underestimates%20it%20where%20the%20curve%20falls%2C%20and%20these%20errors%20cancel%20out%20in%20the%20total.%0A%20%20%20%20See%20that%20the%20rectangle%20and%20trapezoidal%20rules%20are%20indeed%20different%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cumulative_trapezoid%2C%20np)%3A%0A%20%20%20%20x_1%20%3D%20%5B0%2C%202%2C%204%2C%206%2C%208%5D%0A%20%20%20%20print('Integral%20by%20the%20rectangle%20rule%3A'%2C%20np.cumsum(x_1)%20*%201.0)%0A%20%20%20%20print('Integral%20by%20the%20trapezoid%20rule%3A'%2C%20cumulative_trapezoid(x_1%2C%20initial%3D0))%0A%20%20%20%20return%20(x_1%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%20But%20if%20the%20first%20and%20last%20values%20of%20x%20cancel%20out%2C%20as%20here%20when%20the%20mean%20is%20removed%2C%20the%20total%20values%20of%20the%20integral%20using%20the%20cumsum%20and%20cumulative_trapezoid%20methods%20will%20be%20the%20same%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cumulative_trapezoid%2C%20np%2C%20x_1)%3A%0A%20%20%20%20x_2%20%3D%20x_1%20-%20np.mean(x_1)%0A%20%20%20%20print('Integral%20by%20the%20rectangle%20rule%3A'%2C%20np.cumsum(x_2)%20*%201.0)%0A%20%20%20%20print('Integral%20by%20the%20trapezoid%20rule%3A'%2C%20cumulative_trapezoid(x_2%2C%20initial%3D0))%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%20Integrator%20with%20time%20reset%0A%0A%20%20%20%20Another%20measure%20used%20in%20EMG%20processing%20is%20to%20calculate%20the%20integral%20but%20reset%20the%20integrator%20from%20time%20to%20time%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cumulative_trapezoid%2C%20datafrle%2C%20freq%2C%20np)%3A%0A%20%20%20%20nreset%20%3D%20400%20%23%20reset%20after%20this%20amount%20of%20samples%0A%20%20%20%20area%20%3D%20%5B%5D%0A%20%20%20%20for%20i%20in%20range(int(np.ceil(np.size(datafrle)%2Fnreset)))%3A%0A%20%20%20%20%20%20%20%20area%20%3D%20np.hstack((area%2C%20cumulative_trapezoid(datafrle%5Bi*nreset%3A(i%2B1)*nreset%5D%2C%20initial%3D0)%2Ffreq))%0A%20%20%20%20return%20area%2C%20nreset%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%20And%20let's%20plot%20the%20data%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(area%2C%20datafrle%2C%20freq%2C%20nreset%2C%20plt%2C%20time)%3A%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(8%2C%205))%0A%20%20%20%20_ax1.plot(time%2C%20datafrle%2C%20'r')%0A%20%20%20%20_ax1.set_title('EMG%20signal%20(linear%20envelope)')%0A%20%20%20%20_ax1.set_ylabel('EMG%20amplitude%20%5BV%5D')%0A%20%20%20%20_ax1.set_xlim(time%5B0%5D%2C%20time%5B-1%5D)%0A%20%20%20%20_ax2.plot(time%2C%20area%2C%20'y'%2C%20label%3D'Trapezoid')%0A%20%20%20%20_ax2.set_xlabel('Time%20%5Bs%5D')%0A%20%20%20%20_ax2.set_title('Integral%20of%20the%20EMG%20signal%20with%20time%20reset%20(t%20%3D%20%25g%20ms)'%20%25%20(1000%20*%20nreset%20%2F%20freq))%0A%20%20%20%20_ax2.set_ylabel('EMG%20integral%20%5BVs%5D')%0A%20%20%20%20plt.locator_params(axis%3D'both'%2C%20nbins%3D4)%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%204%20%E2%80%94%20Keep%20integration%20units.**%20A%20rectified%20EMG%20remains%20at%200.2%20mV%20for%200.5%20s.%20Give%20the%20integrated%20EMG%20and%20explain%20why%20integrating%20the%20signed%20oscillatory%20trace%20answers%20another%20question.%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%20Solution%20to%20Exercise%204%0A%0A%20%20%20%20The%20rectified%20area%20is%20%240.2(0.5)%3D0.1%24%20mV%20s.%20Signed%20positive%20and%20negative%20lobes%20can%20cancel%2C%20so%20a%20signed%20integral%20measures%20net%20area%20rather%20than%20total%20amplitude%20activity.%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%20Frequency%20estimation%0A%0A%20%20%20%20%23%23%23%20Power%20Spectral%20Density%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dataf%2C%20freq)%3A%0A%20%20%20%20from%20psd2%20import%20psd2%0A%20%20%20%20fpcntile%2C%20mpf%2C%20fmax%2C%20Ptotal%2C%20f%2C%20P%20%3D%20psd2(dataf%2C%20fs%3Dfreq)%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%20%E2%80%94%20Trade%20time%20for%20frequency%20resolution.**%20At%201000%20Hz%2C%20compare%20128-sample%20and%20512-sample%20spectral%20windows.%20Which%20better%20localizes%20a%20short%20burst%2C%20and%20which%20has%20finer%20FFT-bin%20spacing%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%20Solution%20to%20Exercise%205%0A%0A%20%20%20%20The%20windows%20last%200.128%20and%200.512%20s%20and%20have%20bin%20spacings%207.8125%20and%201.953125%20Hz.%20The%20shorter%20window%20localizes%20change%20better%3B%20the%20longer%20window%20resolves%20closer%20frequencies%2C%20subject%20to%20window%20shape%20and%20stationarity.%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%20Short-time%20spectral%20power%0A%0A%20%20%20%20A%20spectrogram%20describes%20how%20spectral%20content%20changes%20over%20time.%20Here%20%60specgram%60%20uses%20its%20default%20PSD%20mode%3A%20the%20returned%20array%20contains%20power%20spectral%20density%2C%20rather%20than%20complex%20short-time%20Fourier%20coefficients.%20Window%20length%20trades%20time%20resolution%20against%20frequency%20resolution.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dataf%2C%20freq%2C%20plt)%3A%0A%20%20%20%20(_fig%2C%20_ax1)%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(8%2C%204))%0A%20%20%20%20(P_1%2C%20freqs%2C%20t%2C%20im)%20%3D%20plt.specgram(dataf%2C%20NFFT%3D128%2C%20Fs%3Dfreq%2C%20noverlap%3D64%2C%20cmap%3Dplt.cm.jet)%0A%20%20%20%20_ax1.set_title('EMG%20spectrogram%20(PSD)'%2C%20fontsize%3D18)%0A%20%20%20%20_ax1.set_xlabel('Time%20%5Bs%5D')%0A%20%20%20%20_ax1.set_ylabel('Frequency%20%5BHz%5D')%0A%20%20%20%20_ax1.set_xlim(t%5B0%5D%2C%20t%5B-1%5D)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20P_1%2C%20freqs%2C%20t%0A%0A%0A%40app.cell%0Adef%20_(P_1%2C%20freqs%2C%20np%2C%20plt%2C%20t)%3A%0A%20%20%20%20(t2%2C%20freqs2)%20%3D%20np.meshgrid(t%2C%20freqs)%0A%20%20%20%20_fig%20%3D%20plt.figure(figsize%3D(10%2C%206))%0A%20%20%20%20_ax%20%3D%20_fig.add_subplot(projection%3D'3d')%0A%20%20%20%20surf%20%3D%20_ax.plot_surface(t2%2C%20freqs2%2C%20P_1%2C%20rstride%3D1%2C%20cstride%3D1%2C%20cmap%3Dplt.cm.jet%2C%20linewidth%3D0%2C%20antialiased%3DFalse)%0A%20%20%20%20_ax.set_xlim(t%5B0%5D%2C%20t%5B-1%5D)%0A%20%20%20%20_ax.set_ylim(0%2C%20500)%0A%20%20%20%20_ax.set_zlim(np.min(P_1)%2C%20np.max(P_1))%0A%20%20%20%20_ax.set_xlabel('Time%20%5Bs%5D'%2C%20fontsize%3D12)%0A%20%20%20%20_ax.set_ylabel('Frequency%20%5BHz%5D'%2C%20fontsize%3D12)%0A%20%20%20%20_ax.set_zlabel('PSD%20%5BV%C2%B2%2FHz%5D'%2C%20fontsize%3D12)%0A%20%20%20%20_ax.set_title('EMG%20spectrogram%20(PSD)'%2C%20fontsize%3D16%2C%20y%3D1.04)%0A%20%20%20%20_ax.view_init(30%2C%2060)%0A%20%20%20%20plt.show()%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%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%201.%20Why%20rectify%20before%20a%20low-pass%20envelope%3F%0A%20%20%20%202.%20Can%20zero-phase%20filtering%20be%20used%20causally%20in%20real%20time%3F%0A%20%20%20%203.%20Does%20increasing%20onset%20threshold%20always%20improve%20detection%3F%0A%20%20%20%204.%20Does%20a%20spectral%20shift%20uniquely%20identify%20muscle%20fatigue%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%0A%20%20%20%201.%20Rectification%20prevents%20positive%20and%20negative%20oscillations%20from%20canceling%20in%20the%20smoother.%202.%20No%3B%20forward%E2%80%93backward%20filtering%20uses%20future%20samples.%203.%20No%3B%20it%20trades%20false%20detections%20against%20missed%2Fdelayed%20events.%204.%20No%3B%20contraction%2C%20conduction%2C%20geometry%2C%20noise%20and%20processing%20also%20affect%20the%20spectrum.%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.%20For%20a%20zero-mean%20sine%20EMG%20with%20peak%20amplitude%201%20mV%2C%20give%20RMS%20and%20mean%20rectified%20amplitude.%0A%0A%20%20%20%202.%20An%20EMG%20is%20sampled%20at%202000%20Hz.%20Give%20the%20nominal%20number%20of%20samples%20in%20a%2025%20ms%20detection%20window.%0A%0A%20%20%20%203.%20A%200.4%20mV%20rectified%20signal%20lasts%200.25%20s.%20Give%20its%20area.%0A%0A%20%20%20%204.%20A%20256-sample%20FFT%20is%20computed%20at%201024%20Hz.%20Give%20bin%20spacing%20and%20Nyquist%20frequency.%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%20to%20Problems%0A%0A%20%20%20%201.%20%241%2F%5Csqrt2%5Capprox0.7071%24%20mV%20and%20%242%2F%5Cpi%5Capprox0.6366%24%20mV.%0A%0A%20%20%20%202.%2050%20samples.%0A%0A%20%20%20%203.%200.1%20mV%20s.%0A%0A%20%20%20%204.%204%20Hz%20and%20512%20Hz.%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%20Use%20%5BKineticsFundamentalConcepts%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fkinetics_fundamental_concepts.py)%20to%20connect%20this%20guide%20with%20the%20surrounding%20mechanics%20and%20modeling%20material.%20Return%20to%20the%20assumptions%20behind%20each%20equation%20when%20applying%20it%20to%20a%20new%20experiment.%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%5BMIT%20OpenCourseWare%5D(https%3A%2F%2Focw.mit.edu%2Fcourses%2Fres-6-007-signals-and-systems-spring-2011%2Fvideo_galleries%2Fvideo-lectures%2F)%20%E2%80%94%20Signals%20and%20Systems%3A%20the%20Fourier-series%2C%20filtering%20and%20sampling%20lectures%20explain%20the%20frequency%20analysis%20and%20processing%20choices%20used%20here.%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-%20De%20Luca%20CJ%20(1997)%20%5BThe%20use%20of%20electromyography%20in%20biomechanics%5D(https%3A%2F%2Fdoi.org%2F10.1123%2Fjab.13.2.135).%20Journal%20of%20Applied%20Biomechanics%2C%20%2013%2C%20135-163.%0A%0A%20%20%20%20-%20Konrad%20P%20(2006)%20%5BThe%20ABC%20of%20EMG%3A%20A%20Practical%20Introduction%20to%20Kinesiological%20Electromyography%5D(https%3A%2F%2Fwww.noraxon.com%2Fdownload%2Fabcs-of-emg%2F).%20Noraxon%20U.S.A.%2C%20Inc.%0A%0A%20%20%20%20-%20Li%20X%2C%20Zhou%20P%2C%20Aruin%20AS%20(2007)%20%5BTeager%E2%80%93Kaiser%20Energy%20Operation%20of%20surface%20EMG%20improves%20muscle%20activity%20onset%20detection%5D(https%3A%2F%2Fpubmed.ncbi.nlm.nih.gov%2F17473984).%20Annals%20of%20Biomedical%20Engineering%2C%2035%2C%201532%E2%80%931538.%0A%0A%20%20%20%20-%20Merletii%20R%20(1999)%20%5BStandards%20for%20Reporting%20EMG%20data%5D(http%3A%2F%2Fwww.isek-online.org%2Fstandards_emg.html).%20Journal%20of%20Electromyography%20and%20Kinesiology%2C%20February%201999%3B%209(1)%3AIII-IV.%0A%0A%20%20%20%20-%20Merletii%20R%2C%20Parker%20PJ%20(2004)%20%5BElectromyography%3A%20Physiology%2C%20Engineering%2C%20and%20Non-Invasive%20Applications%5D(https%3A%2F%2Fbooks.google.com.br%2Fbooks%3Fid%3DSQthgVMil3YC).%20Wiley-IEEE%20Press.%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%26printsec%3Dfrontcover%26source%3Dgbs_ge_summary_r%26cad%3D0%23v%3Donepage%26q%26f%3Dfalse).%204%20ed.%20Hoboken%2C%20EUA%3A%20Wiley.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A