%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%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(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%20Fourier%20transform%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%20A%20signal%20recorded%20in%20time%20can%20also%20be%20described%20by%20the%20frequencies%20it%20contains.%20The%20Fourier%20transform%20is%20the%20mathematical%20tool%20that%20moves%20between%20the%20two%20descriptions%2C%20and%20the%20fast%20Fourier%20transform%20(FFT)%20is%20the%20algorithm%20that%20makes%20it%20practical.%20Filters%2C%20sampling%20rates%2C%20cutoff%20frequencies%20and%20the%20fatigue%20of%20a%20muscle%20are%20all%20discussed%20in%20terms%20of%20frequency%2C%20so%20it%20pays%20to%20be%20able%20to%20compute%20and%20read%20a%20spectrum%20yourself.%0A%0A%20%20%20%20This%20notebook%20continues%20%5BFourier%20series%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Ffourier_series.py).%20It%20defines%20the%20Fourier%20transform%20and%20its%20discrete%20version%2C%20computes%20the%20FFT%20of%20simple%20signals%20and%20reads%20their%20amplitude%20and%20phase%2C%20shows%20the%20traps%20of%20a%20finite%20recording%20(frequency%20resolution%20and%20leakage)%2C%20and%20then%20estimates%20the%20power%20spectral%20density%20of%20noisy%20signals%20and%20of%20a%20real%20muscle.%20It%20ends%20with%20the%20spectrogram%2C%20which%20shows%20how%20the%20frequency%20content%20changes%20in%20time.%20You%20will%20also%20need%20%5BBasic%20properties%20of%20signals%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fsignal_basic_properties.py)%3A%20sampling%2C%20the%20Nyquist%20frequency%2C%20power%20and%20RMS.%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.%20The%20final%20Problems%20have%20collapsed%20final%20answers%20for%20independent%20practice.%20Where%20you%20find%20a%20**Challenge**%20or%20a%20set%20of%20**Guiding%20questions**%2C%20stop%20and%20answer%20on%20a%20scratchpad%20before%20moving%20on.%20Several%20of%20them%20ask%20you%20to%20predict%20a%20number%20*before*%20the%20code%20prints%20it%3B%20the%20prediction%20is%20the%20point%2C%20and%20being%20wrong%20is%20the%20most%20useful%20thing%20that%20can%20happen%20to%20you%20here.%0A%0A%20%20%20%20**Challenge%200.**%20Pick%20a%20signal%20from%20your%20own%20field%3A%20walking%2C%20a%20tremor%2C%20an%20EMG%2C%20the%20force%20under%20a%20runner's%20foot.%20Write%20down%20the%20range%20of%20frequencies%20you%20think%20it%20contains%2C%20from%20the%20lowest%20to%20the%20highest%20that%20matters.%20Keep%20it%3B%20by%20the%20end%20you%20will%20be%20able%20to%20check%20a%20guess%20like%20this%20with%20a%20few%20lines%20of%20code.%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%20has%20the%20FFT%20itself%2C%20SciPy%20the%20spectral%20estimators%2C%20pandas%20reads%20the%20data%20file%2C%20and%20Matplotlib%20draws%20the%20plots.%20The%20declared%20dependencies%20also%20include%20%60psd2%60%20for%20spectral%20summaries.%20Open%20this%20file%20in%20marimo%3B%20the%20EMG%20example%20needs%20access%20to%20the%20repository%20data.%20The%20synthetic%20examples%20need%20no%20data%20download.%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%20%20%20%20from%20scipy%20import%20integrate%2C%20signal%0A%0A%20%20%20%20matplotlib.rc(%22axes%22%2C%20labelsize%3D13%2C%20titlesize%3D14)%0A%20%20%20%20matplotlib.rc(%22xtick%22%2C%20labelsize%3D11)%0A%20%20%20%20matplotlib.rc(%22ytick%22%2C%20labelsize%3D11)%0A%20%20%20%20matplotlib.rc(%22legend%22%2C%20fontsize%3D11)%0A%20%20%20%20return%20np%2C%20pd%2C%20plt%2C%20signal%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%20What%20frequencies%20are%20in%20a%20muscle's%20activity%3F%0A%0A%20%20%20%20Here%20is%20the%20surface%20electromyography%20(EMG)%20used%20in%20%5BBasic%20properties%20of%20signals%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fsignal_basic_properties.py)%20and%20%5BData%20filtering%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fdata_filtering.py)%3A%203.4%20s%20of%20the%20electrical%20activity%20of%20a%20muscle%2C%20sampled%20at%201000%20Hz%2C%20with%20its%20mean%20removed.%20The%20bottom%20panel%20zooms%20in%20on%20100%20ms%20of%20the%20second%20burst%20of%20activity.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20pd%2C%20plt%2C%20repo_file)%3A%0A%20%20%20%20_data%20%3D%20pd.read_csv(%0A%20%20%20%20%20%20%20%20repo_file('data%2Femg.csv')%2C%20header%3DNone%0A%20%20%20%20).to_numpy()%0A%20%20%20%20emg_time%20%3D%20_data%5B%3A%2C%200%5D%0A%20%20%20%20emg%20%3D%20_data%5B%3A%2C%201%5D%20-%20np.mean(_data%5B%3A%2C%201%5D)%20%20%23%20remove%20the%20DC%20component%0A%20%20%20%20emg_fs%20%3D%201%20%2F%20np.mean(np.diff(emg_time))%0A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(10%2C%205.5))%0A%20%20%20%20_axs%5B0%5D.plot(emg_time%2C%20emg%2C%20linewidth%3D0.8)%0A%20%20%20%20_axs%5B0%5D.axvspan(1.5%2C%201.6%2C%20color%3D%22tab%3Aorange%22%2C%20alpha%3D0.3)%0A%20%20%20%20_axs%5B0%5D.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_axs%5B0%5D.set_ylabel(%22EMG%20%5BV%5D%22)%0A%20%20%20%20_zoom%20%3D%20(emg_time%20%3E%3D%201.5)%20%26%20(emg_time%20%3C%3D%201.6)%0A%20%20%20%20_axs%5B1%5D.plot(emg_time%5B_zoom%5D%2C%20emg%5B_zoom%5D%2C%20%22.-%22%2C%20color%3D%22tab%3Aorange%22%2C%20linewidth%3D1)%0A%20%20%20%20_axs%5B1%5D.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_axs%5B1%5D.set_ylabel(%22EMG%20%5BV%5D%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20emg%2C%20emg_fs%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**Guiding%20questions%200.**%0A%0A%20%20%20%201.%20Look%20at%20the%20zoom.%20Count%20roughly%20how%20many%20times%20the%20signal%20crosses%20zero%20in%20those%20100%20ms.%20What%20frequency%20does%20that%20suggest%3F%0A%20%20%20%202.%20Is%20there%20a%20single%20frequency%20in%20this%20signal%2C%20or%20many%3F%20Could%20you%20tell%20their%20relative%20sizes%20from%20the%20plot%3F%0A%20%20%20%203.%20When%20a%20muscle%20fatigues%2C%20its%20EMG%20is%20known%20to%20shift%20towards%20lower%20frequencies.%20How%20would%20you%20measure%20such%20a%20shift%20from%20a%20plot%20like%20this%20one%3F%0A%0A%20%20%20%20The%20time%20plot%20cannot%20answer%20the%20last%20two%20questions.%20A%20spectrum%20can%2C%20and%20by%20the%20end%20of%20this%20notebook%20you%20will%20compute%20this%20EMG's%20median%20frequency%2C%20the%20number%20used%20to%20follow%20muscle%20fatigue.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20The%20Fourier%20transform%0A%0A%20%20%20%20In%20continuation%20to%20%5BFourier%20series%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Ffourier_series.py)%2C%20the%20%5BFourier%20transform%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FFourier_transform)%20is%20a%20mathematical%20transformation%20of%20functions%20between%20the%20time%20(or%20spatial)%20domain%20and%20the%20frequency%20domain.%20Going%20from%20time%20to%20frequency%20is%20called%20Fourier%20analysis%3B%20the%20inverse%20is%20Fourier%20synthesis.%0A%0A%20%20%20%20The%20Fourier%20transform%20of%20a%20continuous%20function%20%24x(t)%24%20is%20by%20definition%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20X(f)%20%3D%20%5Cint_%7B-%5Cinfty%7D%5E%7B%5Cinfty%7D%20x(t)%5C%3A%5Cmathrm%7Be%7D%5E%7B-i2%5Cpi%20ft%7D%20%5C%3A%5Cmathrm%7Bd%7Dt%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%20the%20inverse%20Fourier%20transform%20is%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20x(t)%20%3D%20%5Cint_%7B-%5Cinfty%7D%5E%7B%5Cinfty%7D%20X(f)%5C%3A%5Cmathrm%7Be%7D%5E%7B%5C%3Ai2%5Cpi%20tf%7D%20%5C%3A%5Cmathrm%7Bd%7Df%0A%20%20%20%20%24%24%0A%0A%20%20%20%20A%20Fourier%20series%20describes%20a%20*periodic*%20function%20by%20a%20discrete%20set%20of%20harmonics%20of%20its%20fundamental%20frequency.%20The%20Fourier%20transform%20extends%20the%20idea%20to%20functions%20that%20need%20not%20be%20periodic%3A%20think%20of%20it%20as%20the%20period%20growing%20without%20limit%2C%20so%20that%20the%20harmonics%20become%20infinitely%20close%20and%20the%20sum%20becomes%20an%20integral%20over%20a%20continuum%20of%20frequencies.%20%24X(f)%24%20is%20complex%3A%20its%20magnitude%20says%20how%20much%20of%20the%20frequency%20%24f%24%20the%20signal%20contains%2C%20and%20its%20angle%20says%20the%20phase%20of%20that%20component.%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%20Discrete%20Fourier%20transform%0A%0A%20%20%20%20A%20recorded%20signal%20is%20neither%20continuous%20nor%20infinite%3A%20it%20is%20%24N%24%20samples%20%24x%5Bn%5D%24%2C%20taken%20every%20%24%5CDelta%20t%20%3D%201%2Ff_s%24%20seconds.%20Its%20Discrete%20Fourier%20Transform%20(DFT)%20is%20another%20sequence%20%24X%24%2C%20also%20with%20%24N%24%20elements%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20X%5Bk%5D%20%3D%20%5Csum_%7Bn%3D0%7D%5E%7BN-1%7D%20%20x%5Bn%5D%5C%2C%20%5Cmathrm%7Be%7D%5E%7B-i2%5Cpi%20kn%2FN%7D%20%5C%3B%2C%5Cquad%200%20%5Cleq%20k%20%5Cleq%20N-1%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20Inverse%20Discrete%20Fourier%20Transform%20(IDFT)%20inverts%20this%20operation%20and%20gives%20back%20the%20original%20data%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20x%5Bn%5D%20%3D%20%5Cfrac%7B1%7D%7BN%7D%20%5Csum_%7Bk%3D0%7D%5E%7BN-1%7D%20%20X%5Bk%5D%5C%2C%20%5Cmathrm%7Be%7D%5E%7Bi2%5Cpi%20kn%2FN%7D%20%5C%3B%2C%5Cquad%200%20%5Cleq%20n%20%5Cleq%20N-1%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20element%20%24X%5Bk%5D%24%20corresponds%20to%20the%20frequency%0A%0A%20%20%20%20%24%24%0A%20%20%20%20f_k%20%3D%20%5Cfrac%7Bk%7D%7BN%5CDelta%20t%7D%20%3D%20k%5C%2C%5Cfrac%7Bf_s%7D%7BN%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20so%20the%20frequencies%20in%20a%20DFT%20are%20spaced%20by%20%24%5CDelta%20f%20%3D%201%2F(N%5CDelta%20t)%24%2C%20the%20inverse%20of%20the%20duration%20of%20the%20recording.%20This%20is%20the%20**DFT%20bin%20spacing**%3A%20a%205%20s%20record%20has%20bins%200.2%20Hz%20apart%2C%20whatever%20its%20sampling%20rate.%20Resolving%20two%20nearby%20components%20also%20depends%20on%20the%20window%2C%20their%20amplitudes%20and%20noise%3B%20bin%20spacing%20alone%20does%20not%20guarantee%20that%20they%20will%20be%20distinguishable.%20The%20elements%20above%20%24k%20%3D%20N%2F2%24%2C%20beyond%20the%20Nyquist%20frequency%20%24f_s%2F2%24%2C%20are%20the%20negative%20frequencies%20%24f_k%20-%20f_s%24%20(see%20the%20section%20on%20aliasing%20in%20%5BBasic%20properties%20of%20signals%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fsignal_basic_properties.py)).%0A%0A%20%20%20%20Writing%20the%20complex%20exponential%20of%20the%20inverse%20transform%20as%20%24%5Ccos%20%2B%20i%5Csin%24%2C%20a%20real%20signal%20can%20be%20written%20as%20a%20sum%20of%20sinusoids%2C%20as%20in%20a%20Fourier%20series%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20x%5Bn%5D%20%3D%20%5Csum_%7Bk%3D0%7D%5E%7BN-1%7D%20a_k%5Ccos%5Cleft(%5Cfrac%7B2%5Cpi%20kn%7D%7BN%7D%5Cright)%2Bb_k%5Csin%5Cleft(%5Cfrac%7B2%5Cpi%20kn%7D%7BN%7D%5Cright)%20%5C%3B%2C%5Cquad%200%20%5Cleq%20n%20%5Cleq%20N-1%0A%20%20%20%20%24%24%0A%0A%20%20%20%20with%20the%20Fourier%20coefficients%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20a_k%20%3D%20%5Cfrac%7B%5Ctext%7BReal%7D(X%5Bk%5D)%7D%7BN%7D%20%5C%3B%2C%20%5Cqquad%20b_k%20%3D%20-%5Cfrac%7B%5Ctext%7BImag%7D(X%5Bk%5D)%7D%7BN%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%20%24a_0%20%3D%20X%5B0%5D%2FN%24%20is%20the%20mean%20of%20the%20data%2C%20its%20DC%20component.%0A%0A%20%20%20%20**Guiding%20questions%201.**%0A%0A%20%20%20%201.%20What%20is%20the%20frequency%20resolution%20of%20a%202%20s%20recording%20sampled%20at%201000%20Hz%3F%20And%20of%20a%2020%20s%20recording%20sampled%20at%20100%20Hz%3F%0A%20%20%20%202.%20For%20the%20EMG%20above%2C%20%24N%20%3D%203360%24%20and%20%24f_s%20%3D%201000%24%20Hz.%20What%20frequency%20does%20%24X%5B100%5D%24%20correspond%20to%3F%20And%20%24X%5B3000%5D%24%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Exercise%201.**%20With%20500%20samples%20at%20100%20Hz%2C%20locate%20the%20DFT%20bins%20for%20DC%2C%205%20Hz%2C%20and%20Nyquist.%20What%20record%20duration%20would%20halve%20the%20bin%20spacing%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%201%0A%20%20%20%20The%20duration%20is%20%24N%2Ff_s%3D5%24%20s%20and%20%24%5CDelta%20f%3Df_s%2FN%3D0.2%24%20Hz.%20The%20indices%20are%20%24k%3D0%24%2C%20%24k%3D25%24%20and%20%24k%3D250%24.%20Double%20the%20record%20duration%20to%2010%20s%20to%20obtain%200.1%20Hz%20spacing.%20Increasing%20only%20the%20sampling%20rate%20while%20retaining%205%20s%20does%20not%20do%20this.%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%20Fast%20Fourier%20Transform%20(FFT)%0A%0A%20%20%20%20Computed%20from%20its%20definition%2C%20the%20DFT%20takes%20about%20%24N%5E2%24%20operations.%20The%20%5BFFT%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FFast_Fourier_transform)%20is%20a%20family%20of%20fast%20algorithms%20that%20compute%20exactly%20the%20same%20DFT%20in%20about%20%24N%5Clog_2%20N%24%20operations%3A%20for%20a%20million%20samples%2C%20the%20difference%20between%20minutes%20and%20milliseconds.%20NumPy%20has%20it%20in%20%60numpy.fft%60%20(SciPy%20has%20an%20equivalent%20%60scipy.fft%60%3B%20the%20older%20%60scipy.fftpack%60%20is%20legacy).%0A%0A%20%20%20%20Our%20first%20signal%20is%20a%20sine%20wave%20with%20an%20amplitude%20of%202%2C%20a%20frequency%20of%205%20Hz%20and%20a%20phase%20of%20%2445%5Eo%24%2C%20plus%20a%20DC%20component%20of%201%2C%20sampled%20at%20100%20Hz%20for%205%20s%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt)%3A%0A%20%20%20%20fs%20%3D%20100%20%20%23%20sampling%20frequency%20%5BHz%5D%0A%20%20%20%20dt%20%3D%201%20%2F%20fs%0A%20%20%20%20t%20%3D%20np.arange(0%2C%20500)%20%2F%20fs%20%20%23%205%20s%0A%20%20%20%20x%20%3D%202%20*%20np.sin(2%20*%20np.pi%20*%205%20*%20t%20%2B%20np.pi%20%2F%204)%20%2B%201%0A%20%20%20%20N%20%3D%20x.size%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(9%2C%203))%0A%20%20%20%20_ax.plot(t%2C%20x%2C%20linewidth%3D2)%0A%20%20%20%20_ax.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22Amplitude%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20N%2C%20dt%2C%20fs%2C%20t%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%20Its%20FFT%20is%20simply%20%60np.fft.fft(x)%60%2C%20and%20%60np.fft.fftfreq%60%20gives%20the%20frequency%20of%20each%20element.%20The%20amplitude%20of%20each%20component%20is%20the%20magnitude%20of%20%24X%5Bk%5D%24%20divided%20by%20%24N%24%2C%20and%20its%20phase%20is%20the%20angle%20of%20%24X%5Bk%5D%24.%20The%20phase%20of%20a%20component%20with%20no%20amplitude%20is%20meaningless%20(it%20is%20the%20angle%20of%20a%20number%20that%20is%20zero%20up%20to%20rounding%20errors)%2C%20so%20we%20set%20it%20to%20zero%20wherever%20the%20amplitude%20is%20negligible.%0A%0A%20%20%20%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20at%20which%20frequencies%20there%20will%20be%20peaks%2C%20and%20their%20heights.%20Think%20of%20%24%5Csin%5Ctheta%20%3D%20(%5Cmathrm%7Be%7D%5E%7Bi%5Ctheta%7D%20-%20%5Cmathrm%7Be%7D%5E%7B-i%5Ctheta%7D)%2F2i%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(N%2C%20dt%2C%20np%2C%20plt%2C%20x)%3A%0A%20%20%20%20X%20%3D%20np.fft.fft(x)%20%20%23%20FFT%0A%20%20%20%20freqs%20%3D%20np.fft.fftfreq(N%2C%20dt)%20%20%23%20frequency%20of%20each%20element%20%5BHz%5D%0A%20%20%20%20amp%20%3D%20np.abs(X)%20%2F%20N%20%20%23%20amplitude%0A%20%20%20%20phase%20%3D%20np.where(amp%20%3E%201e-6%20*%20amp.max()%2C%20np.angle(X%2C%20deg%3DTrue)%2C%200)%20%20%23%20phase%20%5Bdegrees%5D%0A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(9%2C%205)%2C%20sharex%3DTrue)%0A%20%20%20%20_axs%5B0%5D.plot(np.fft.fftshift(freqs)%2C%20np.fft.fftshift(amp)%2C%20%22.-%22)%0A%20%20%20%20_axs%5B1%5D.plot(np.fft.fftshift(freqs)%2C%20np.fft.fftshift(phase)%2C%20%22.-%22)%0A%20%20%20%20_axs%5B0%5D.set_ylabel(%22Amplitude%22)%0A%20%20%20%20_axs%5B0%5D.set_ylim(-0.01%2C%201.1)%0A%20%20%20%20_axs%5B1%5D.set_ylabel(%22Phase%20%5B%24%5Eo%24%5D%22)%0A%20%20%20%20_axs%5B1%5D.set_xlabel(%22Frequency%20%5BHz%5D%22)%0A%20%20%20%20_axs%5B1%5D.set_xlim(-50%2C%2050)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%0A%20%20%20%20for%20_k%20in%20np.nonzero(amp%20%3E%200.01)%5B0%5D%3A%0A%20%20%20%20%20%20%20%20print(f%22f%20%3D%20%7Bfreqs%5B_k%5D%3A5.1f%7D%20Hz%3A%20amplitude%20%3D%20%7Bamp%5B_k%5D%3A.3f%7D%2C%20phase%20%3D%20%7Bphase%5B_k%5D%3A6.1f%7D%20degrees%22)%0A%20%20%20%20return%20(X%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%20Three%20peaks%3A%20the%20DC%20component%2C%20with%20its%20full%20value%20of%201%2C%20and%20two%20peaks%20of%201%2C%20half%20the%20amplitude%20of%20the%20sine%2C%20at%20%24%2B5%24%20and%20%24-5%24%20Hz.%20A%20real%20sinusoid%20is%20the%20sum%20of%20two%20complex%20exponentials%2C%20one%20at%20each%20frequency%2C%20and%20each%20carries%20half%20of%20it.%0A%0A%20%20%20%20The%20phase%20is%20%24-45%5Eo%24%20at%20%24%2B5%24%20Hz%2C%20not%20%24%2B45%5Eo%24.%20The%20FFT%20measures%20phase%20relative%20to%20a%20*cosine*%2C%20and%20%242%5Csin(%5Comega%20t%20%2B%2045%5Eo)%20%3D%202%5Ccos(%5Comega%20t%20-%2045%5Eo)%24.%20Keep%20that%20convention%20in%20mind%20whenever%20you%20read%20a%20phase%20off%20an%20FFT.%0A%0A%20%20%20%20For%20any%20real%20signal%2C%20%24X%5B-k%5D%24%20is%20the%20complex%20conjugate%20of%20%24X%5Bk%5D%24%3A%20the%20amplitudes%20at%20negative%20and%20positive%20frequencies%20are%20the%20same%2C%20and%20the%20phases%20have%20opposite%20signs%2C%20as%20in%20the%20plot.%20(If%20the%20signal%20is%20also%20even%2C%20%24X%24%20is%20real%3B%20if%20it%20is%20odd%2C%20%24X%24%20is%20imaginary.)%20The%20negative%20frequencies%20carry%20no%20new%20information%2C%20so%20we%20usually%20plot%20only%20the%20positive%20ones%20and%20double%20their%20amplitudes%2C%20except%20for%20DC%20and%2C%20for%20an%20even%20number%20of%20samples%2C%20the%20Nyquist%20component%2C%20each%20of%20which%20appears%20only%20once.%20This%20is%20the%20one-sided%20spectrum%2C%20which%20NumPy%20computes%20directly%20with%20%60np.fft.rfft%60%20and%20%60np.fft.rfftfreq%60%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(N%2C%20dt%2C%20np%2C%20plt%2C%20x)%3A%0A%20%20%20%20_X2%20%3D%20np.fft.rfft(x)%0A%20%20%20%20_freqs2%20%3D%20np.fft.rfftfreq(N%2C%20dt)%0A%20%20%20%20_amp2%20%3D%20np.abs(_X2)%20%2F%20N%0A%20%20%20%20_amp2%5B1%3A-1%20if%20N%20%25%202%20%3D%3D%200%20else%20None%5D%20*%3D%202%0A%20%20%20%20_phase2%20%3D%20np.angle(_X2%2C%20deg%3DTrue)%0A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(9%2C%205)%2C%20sharex%3DTrue)%0A%20%20%20%20_axs%5B0%5D.plot(_freqs2%2C%20_amp2%2C%20%22.-%22)%0A%20%20%20%20_axs%5B1%5D.plot(_freqs2%2C%20_phase2%2C%20%22.-%22)%0A%20%20%20%20_axs%5B0%5D.set_ylabel(%22Amplitude%22)%0A%20%20%20%20_axs%5B0%5D.set_ylim(-0.01%2C%202.1)%0A%20%20%20%20_axs%5B1%5D.set_ylabel(%22Phase%20%5B%24%5Eo%24%5D%22)%0A%20%20%20%20_axs%5B1%5D.set_xlabel(%22Frequency%20%5BHz%5D%22)%0A%20%20%20%20_axs%5B1%5D.set_xlim(-0.1%2C%2050)%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%20Now%20the%20amplitude%20at%205%20Hz%20is%20the%20amplitude%20of%20the%20sine%2C%202%2C%20and%20the%20DC%20component%20is%201.%0A%0A%20%20%20%20%23%23%23%20Fourier%20synthesis%0A%0A%20%20%20%20We%20can%20get%20the%20data%20back%20from%20the%20Fourier%20coefficients%2C%20by%20adding%20up%20all%20%24N%24%20sinusoids%20of%20the%20formula%20above.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20how%20close%20the%20synthesis%20will%20be%20to%20the%20original%20data.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(N%2C%20X%2C%20dt%2C%20np%2C%20plt%2C%20t%2C%20x)%3A%0A%20%20%20%20%23%20Fourier%20coefficients%0A%20%20%20%20_a%20%3D%20np.real(X)%20%2F%20N%0A%20%20%20%20_b%20%3D%20-np.imag(X)%20%2F%20N%0A%0A%20%20%20%20%23%20Fourier%20synthesis%3A%20a%20sum%20of%20N%20sinusoids%2C%20one%20per%20column%0A%20%20%20%20_w%20%3D%202%20*%20np.pi%20*%20np.arange(N)%20%2F%20(N%20*%20dt)%20%20%23%20angular%20frequencies%0A%20%20%20%20_y%20%3D%20_a%20*%20np.cos(np.outer(t%2C%20_w))%20%2B%20_b%20*%20np.sin(np.outer(t%2C%20_w))%0A%20%20%20%20_xfft%20%3D%20np.sum(_y%2C%20axis%3D1)%0A%0A%20%20%20%20%23%20or%2C%20simply%2C%20the%20inverse%20FFT%0A%20%20%20%20_xifft%20%3D%20np.real(np.fft.ifft(X))%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(9%2C%203))%0A%20%20%20%20_ax.plot(t%2C%20x%2C%20linewidth%3D3%2C%20label%3D%22Original%20data%22)%0A%20%20%20%20_ax.plot(t%2C%20_xfft%2C%20%22r--%22%2C%20linewidth%3D2%2C%20label%3D%22Fourier%20synthesis%22)%0A%20%20%20%20_ax.plot(t%2C%20_xifft%2C%20%22k%3A%22%2C%20linewidth%3D2%2C%20label%3D%22Inverse%20FFT%22)%0A%20%20%20%20_ax.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22Amplitude%22)%0A%20%20%20%20_ax.legend(framealpha%3D0.7%2C%20loc%3D%22upper%20right%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20print(f%22Largest%20difference%2C%20synthesis%3A%20%7Bnp.max(np.abs(_xfft%20-%20x))%3A.1e%7D%22)%0A%20%20%20%20print(f%22Largest%20difference%2C%20inverse%20FFT%3A%20%7Bnp.max(np.abs(_xifft%20-%20x))%3A.1e%7D%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Identical%2C%20to%20the%20rounding%20errors%20of%20the%20computer.%20The%20DFT%20loses%20nothing%3A%20the%20%24N%24%20complex%20numbers%20%24X%5Bk%5D%24%20are%20another%2C%20complete%20description%20of%20the%20%24N%24%20samples.%20(Between%20the%20samples%20it%20is%20another%20matter%3A%20the%20sum%20of%20sinusoids%20at%20frequencies%20above%20the%20Nyquist%20frequency%20does%20not%20interpolate%20the%20data%20sensibly.%20The%20DFT%20describes%20the%20samples%2C%20not%20the%20continuous%20signal%20they%20came%20from.)%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.**%20For%20%24x(t)%3D1%2B2%5Csin(2%5Cpi5t%2B%5Cpi%2F4)%24%2C%20explain%20the%20positive-frequency%20FFT%20amplitude%20and%20phase.%20What%20must%20a%20one-sided%20plot%20do%20to%20DC%20and%20Nyquist%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%202%0A%20%20%20%20Write%20the%20sine%20as%20%242%5Ccos(2%5Cpi5t-%5Cpi%2F4)%24.%20The%20conjugate%20frequency%20pair%20each%20has%20amplitude%201%20after%20division%20by%20%24N%24%3B%20combining%20the%20pair%20gives%20amplitude%202%20and%20cosine%20phase%20%24-45%5E%5Ccirc%24.%20DC%20remains%201.%20DC%20and%20the%20even-length%20Nyquist%20bin%20have%20no%20distinct%20partners%2C%20so%20doubling%20either%20would%20be%20an%20error.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Frequency%20resolution%20and%20leakage%0A%0A%20%20%20%20The%205%20Hz%20sine%20fitted%20exactly%2025%20periods%20in%20the%205%20s%20of%20data%2C%20and%205%20Hz%20is%20exactly%20one%20of%20the%20DFT%20frequencies%2C%20%24k%5C%2C%5CDelta%20f%24%20with%20%24%5CDelta%20f%20%3D%200.2%24%20Hz.%20Real%20signals%20are%20never%20that%20polite.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20what%20the%20one-sided%20spectrum%20of%20a%20sine%20of%20amplitude%202%20at%205.1%20Hz%2C%20halfway%20between%20two%20DFT%20frequencies%2C%20looks%20like.%20Will%20its%20peak%20still%20be%202%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(N%2C%20fs%2C%20np%2C%20plt%2C%20signal%2C%20t)%3A%0A%20%20%20%20_f%20%3D%20np.fft.rfftfreq(N%2C%201%20%2F%20fs)%0A%20%20%20%20_hann%20%3D%20signal.windows.hann(N)%0A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(11%2C%203.8)%2C%20sharey%3DTrue)%0A%20%20%20%20for%20_ax%2C%20_freq%20in%20zip(_axs%2C%20(5%2C%205.1))%3A%0A%20%20%20%20%20%20%20%20_x%20%3D%202%20*%20np.sin(2%20*%20np.pi%20*%20_freq%20*%20t)%0A%20%20%20%20%20%20%20%20_amp%20%3D%202%20*%20np.abs(np.fft.rfft(_x))%20%2F%20N%0A%20%20%20%20%20%20%20%20_amp_hann%20%3D%202%20*%20np.abs(np.fft.rfft(_x%20*%20_hann))%20%2F%20np.sum(_hann)%0A%20%20%20%20%20%20%20%20_amp%5B%5B0%2C%20-1%5D%5D%20%2F%3D%202%20%20%23%20N%20is%20even%3A%20do%20not%20double%20DC%20or%20Nyquist%0A%20%20%20%20%20%20%20%20_amp_hann%5B%5B0%2C%20-1%5D%5D%20%2F%3D%202%0A%20%20%20%20%20%20%20%20_ax.plot(_f%2C%20_amp%2C%20%22o-%22%2C%20label%3D%22no%20window%22)%0A%20%20%20%20%20%20%20%20_ax.plot(_f%2C%20_amp_hann%2C%20%22s--%22%2C%20label%3D%22Hann%20window%22)%0A%20%20%20%20%20%20%20%20_ax.set_xlim(3%2C%207)%0A%20%20%20%20%20%20%20%20_ax.set_title(f%22Sine%20of%20amplitude%202%20at%20%7B_freq%7D%20Hz%22)%0A%20%20%20%20%20%20%20%20_ax.set_xlabel(%22Frequency%20%5BHz%5D%22)%0A%20%20%20%20%20%20%20%20print(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22%7B_freq%7D%20Hz%3A%20peak%20%7B_amp.max()%3A.2f%7D%20(Hann%20%7B_amp_hann.max()%3A.2f%7D)%3B%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22amplitude%20at%2010%20Hz%20%7B_amp%5B_f%20%3D%3D%2010%5D%5B0%5D%3A.4f%7D%20(Hann%20%7B_amp_hann%5B_f%20%3D%3D%2010%5D%5B0%5D%3A.5f%7D)%22%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20_axs%5B0%5D.set_ylabel(%22Amplitude%22)%0A%20%20%20%20_axs%5B0%5D.legend(loc%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%20At%205.1%20Hz%20the%20peak%20falls%20to%20about%201.3%2C%20and%20the%20rest%20of%20the%20amplitude%20*leaks*%20into%20the%20neighbouring%20frequencies%2C%20spreading%20over%20the%20whole%20spectrum%3A%20even%2010%20Hz%20gets%20some.%20Nothing%20is%20wrong%20with%20the%20sine.%20The%20DFT%20implicitly%20treats%20the%205%20s%20of%20data%20as%20one%20period%20of%20a%20periodic%20signal%2C%20and%205.1%20Hz%20does%20not%20fit%20a%20whole%20number%20of%20periods%20in%205%20s%3A%20the%20implied%20periodic%20signal%20jumps%20where%20the%20end%20meets%20the%20beginning%2C%20and%20a%20jump%20contains%20all%20frequencies.%0A%0A%20%20%20%20The%20usual%20remedy%20is%20a%20**window**%3A%20multiplying%20the%20data%20by%20a%20function%20that%20goes%20smoothly%20to%20zero%20at%20both%20ends%2C%20such%20as%20the%20Hann%20window%2C%20removes%20the%20jump.%20The%20leakage%20far%20from%20the%20peak%20falls%20by%20orders%20of%20magnitude%3B%20the%20peak%20itself%20becomes%20wider%2C%20and%20with%20a%20correction%20for%20the%20window's%20area%20its%20height%20is%20closer%20to%202.%20Windows%20trade%20resolution%20for%20leakage%2C%20and%20the%20spectral%20estimators%20below%20use%20one%20by%20default.%0A%0A%20%20%20%20**Challenge%201.**%20Change%205.1%20Hz%20to%205.2%20Hz%20in%20the%20cell%20above.%20Then%20keep%205.1%20Hz%20but%20use%2010%20s%20of%20data%20instead%20of%205%20s.%20Explain%20both%20results%20with%20%24%5CDelta%20f%20%3D%201%2F(N%5CDelta%20t)%24.%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.**%20Answer%20Challenge%201%3A%20why%20does%205.2%20Hz%20fit%20the%205%20s%20record%2C%20and%20why%20does%205.1%20Hz%20fit%20a%2010%20s%20record%3F%20Does%20a%20Hann%20window%20create%20more%20frequency%20resolution%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%203%0A%20%20%20%20%245.2%5Ctimes5%3D26%24%20and%20%245.1%5Ctimes10%3D51%24%20are%20whole%20numbers%20of%20cycles%2C%20so%20both%20tones%20fall%20on%20their%20respective%20DFT%20grids.%20A%20Hann%20window%20reduces%20distant%20sidelobes%20by%20tapering%20the%20boundaries%2C%20but%20widens%20the%20main%20lobe.%20It%20trades%20leakage%20suppression%20for%20separation%20of%20close%20components%3B%20it%20does%20not%20add%20information.%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%20A%20signal%20buried%20in%20noise%0A%0A%20%20%20%20Another%20example%3A%20the%20sum%20of%20two%20sines%2C%20of%20amplitudes%202%20and%201%20at%205%20and%2020%20Hz%2C%20plus%20random%20noise%20with%20a%20standard%20deviation%20of%201%2C%20sampled%20at%20100%20Hz%20for%205%20s.%20**Before%20you%20run%20the%20next%20two%20cells**%2C%20predict%20whether%20you%20will%20be%20able%20to%20see%20the%20two%20sines%20in%20the%20time%20plot%2C%20and%20in%20the%20spectrum.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt)%3A%0A%20%20%20%20freq_n%20%3D%20100.0%20%20%23%20sampling%20frequency%20%5BHz%5D%0A%20%20%20%20t_n%20%3D%20np.arange(0%2C%205%2C%201%20%2F%20freq_n)%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(seed%3D42)%0A%20%20%20%20y_n%20%3D%20(%0A%20%20%20%20%20%20%20%202%20*%20np.sin(2%20*%20np.pi%20*%205%20*%20t_n)%0A%20%20%20%20%20%20%20%20%2B%20np.sin(2%20*%20np.pi%20*%2020%20*%20t_n)%0A%20%20%20%20%20%20%20%20%2B%20_rng.standard_normal(t_n.size)%0A%20%20%20%20)%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(9%2C%203))%0A%20%20%20%20_ax.set_title(%22Time%20domain%22)%0A%20%20%20%20_ax.plot(t_n%2C%20y_n%2C%20%22b%22%2C%20linewidth%3D1)%0A%20%20%20%20_ax.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22y%20%5BV%5D%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20freq_n%2C%20t_n%2C%20y_n%0A%0A%0A%40app.cell%0Adef%20_(freq_n%2C%20np%2C%20plt%2C%20y_n)%3A%0A%20%20%20%20%23%20one-sided%20amplitude%20spectrum%0A%20%20%20%20_amp%20%3D%20np.abs(np.fft.rfft(y_n))%20%2F%20y_n.size%0A%20%20%20%20_amp%5B1%3A-1%20if%20y_n.size%20%25%202%20%3D%3D%200%20else%20None%5D%20*%3D%202%0A%20%20%20%20_freqs%20%3D%20np.fft.rfftfreq(y_n.size%2C%201%20%2F%20freq_n)%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(9%2C%203))%0A%20%20%20%20_ax.set_title(%22Frequency%20domain%22)%0A%20%20%20%20_ax.plot(_freqs%2C%20_amp%2C%20%22r%22%2C%20linewidth%3D1.5)%0A%20%20%20%20_ax.set_xlabel(%22Frequency%20%5BHz%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22Amplitude%20%5BV%5D%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20_top%20%3D%20np.argsort(_amp)%5B%3A%3A-1%5D%5B%3A3%5D%0A%20%20%20%20print(%22Three%20largest%20peaks%3A%22%2C%20%22%2C%20%22.join(f%22%7B_freqs%5B_i%5D%3A.1f%7D%20Hz%20(%7B_amp%5B_i%5D%3A.2f%7D%20V)%22%20for%20_i%20in%20_top))%0A%20%20%20%20print(f%22Median%20amplitude%3A%20%7Bnp.median(_amp)%3A.3f%7D%20V%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20In%20the%20time%20plot%20the%20sines%20are%20hard%20to%20make%20out%3B%20in%20the%20spectrum%20they%20stand%20out%20clearly%2C%20at%205%20and%2020%20Hz%2C%20with%20amplitudes%20close%20to%202%20and%201.%20The%20noise%20has%20twice%20the%20power%20of%20the%2020%20Hz%20sine%20(a%20variance%20of%201%2C%20against%20%241%5E2%2F2%24)%2C%20but%20it%20is%20spread%20over%20all%20250%20frequencies%20of%20the%20spectrum%2C%20while%20each%20sine%20is%20concentrated%20in%20one.%20That%20is%20the%20power%20of%20the%20frequency%20domain%3A%20it%20separates%20what%20is%20concentrated%20in%20frequency%20from%20what%20is%20spread%20out.%0A%0A%20%20%20%20%23%23%23%20FFTW%2C%20the%20Fastest%20Fourier%20Transform%20in%20the%20West%0A%0A%20%20%20%20%5BFFTW%5D(https%3A%2F%2Fwww.fftw.org%2F)%20is%20a%20free%20collection%20of%20fast%20C%20routines%20for%20computing%20the%20DFT%2C%20among%20the%20fastest%20available.%20NumPy%20and%20SciPy%20use%20their%20own%20fast%20implementations%2C%20which%20are%20more%20than%20enough%20for%20the%20signals%20in%20biomechanics%3B%20if%20speed%20ever%20becomes%20a%20concern%2C%20the%20Python%20wrapper%20%5BpyFFTW%5D(https%3A%2F%2Fpypi.org%2Fproject%2FpyFFTW%2F)%20gives%20access%20to%20FFTW.%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.**%20A%201%20V%20amplitude%20sine%20is%20mixed%20with%20zero-mean%20noise%20of%20variance%20%240.2%24%20V%C2%B2.%20What%20mean-square%20power%20should%20a%20PSD%20recover%20on%20average%3F%20Is%20the%20units%20label%20V%C2%B2%20correct%20for%20PSD%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%204%0A%20%20%20%20The%20sine%20contributes%20%241%5E2%2F2%3D0.5%24%20V%C2%B2%20and%20independent%20zero-mean%20noise%20contributes%20%240.2%24%20V%C2%B2%2C%20so%20expected%20total%20power%20is%20%240.7%24%20V%C2%B2.%20A%20density%20has%20units%20V%C2%B2%2FHz%3B%20multiplying%20each%20bin%20by%20its%20width%20in%20Hz%20and%20summing%20gives%20V%C2%B2.%20Finite%20records%20fluctuate%20around%20the%20expectation.%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%20Power%20spectral%20density%0A%0A%20%20%20%20For%20noisy%20signals%2C%20the%20amplitude%20of%20each%20frequency%20is%20less%20useful%20than%20the%20**power%20spectral%20density**%20(PSD)%2C%20the%20power%20of%20the%20signal%20per%20unit%20of%20frequency%2C%20in%20units%20squared%20per%20hertz%20(for%20instance%2C%20V%C2%B2%2FHz).%20Its%20integral%20over%20all%20frequencies%20is%20the%20total%20power%20of%20the%20signal%2C%20that%20is%2C%20its%20mean%20squared%20value%20or%2C%20without%20the%20DC%20component%2C%20its%20variance%20(%5BParseval's%20theorem%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FParseval%2527s_theorem)).%0A%0A%20%20%20%20The%20simplest%20estimate%20of%20the%20PSD%20is%20the%20%5Bperiodogram%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FPeriodogram)%2C%20the%20squared%20magnitude%20of%20the%20FFT%2C%20scaled%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20P%5Bk%5D%20%3D%20%5Cfrac%7B%7CX%5Bk%5D%7C%5E2%7D%7BN%20f_s%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20doubled%20for%20the%20positive%20frequencies%20other%20than%20DC%20and%20Nyquist%2C%20as%20for%20the%20one-sided%20amplitude.%20The%20periodogram%20of%20a%20noisy%20signal%20is%20itself%20very%20noisy%3A%20each%20of%20its%20values%20is%20computed%20from%20a%20single%20%22sample%22%20of%20the%20noise.%20%5BWelch's%20method%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FWelch%2527s_method)%20reduces%20that%20variance%20by%20splitting%20the%20data%20into%20overlapping%20segments%2C%20computing%20a%20windowed%20periodogram%20of%20each%2C%20and%20averaging%20them%3B%20the%20price%20is%20a%20coarser%20frequency%20resolution%2C%20set%20by%20the%20length%20of%20the%20segments.%0A%0A%20%20%20%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20the%20total%20power%20of%20the%20noisy%20signal%20above%3A%20the%20power%20of%20a%20sine%20of%20amplitude%20%24A%24%20is%20%24A%5E2%2F2%24%2C%20and%20the%20power%20of%20the%20noise%20is%20its%20variance.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(freq_n%2C%20np%2C%20plt%2C%20signal%2C%20t_n%2C%20y_n)%3A%0A%20%20%20%20_N%20%3D%20y_n.size%0A%20%20%20%20_fp%2C%20_Pp%20%3D%20signal.periodogram(y_n%2C%20freq_n%2C%20window%3D%22boxcar%22%2C%20nfft%3D_N)%0A%20%20%20%20_fw%2C%20_Pw%20%3D%20signal.welch(y_n%2C%20freq_n%2C%20window%3D%22hann%22%2C%20nperseg%3D_N%20%2F%2F%204)%0A%0A%20%20%20%20%23%20quick%20and%20simple%20PSD%2C%20which%20is%20the%20periodogram%0A%20%20%20%20_P%20%3D%20np.abs(np.fft.rfft(y_n%20-%20np.mean(y_n)))%20**%202%20%2F%20_N%20%2F%20freq_n%0A%20%20%20%20_P%5B1%3A-1%5D%20%3D%202%20*%20_P%5B1%3A-1%5D%20%20%23%20one-sided%2C%20N%20even%3A%20DC%20and%20Nyquist%20appear%20once%0A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(3%2C%201%2C%20figsize%3D(10%2C%208))%0A%20%20%20%20_axs%5B0%5D.set_title(%22Time%20domain%22)%0A%20%20%20%20_axs%5B0%5D.plot(t_n%2C%20y_n%2C%20%22b%22%2C%20linewidth%3D1)%0A%20%20%20%20_axs%5B0%5D.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_axs%5B0%5D.set_ylabel(%22y%20%5BV%5D%22)%0A%20%20%20%20_axs%5B1%5D.set_title(%22Periodogram%22)%0A%20%20%20%20_axs%5B1%5D.plot(_fp%2C%20_Pp%2C%20%22r%22%2C%20linewidth%3D1.5)%0A%20%20%20%20_axs%5B1%5D.set_ylabel(%22PSD(y)%20%5BV%24%5E2%24%2FHz%5D%22)%0A%20%20%20%20_axs%5B2%5D.set_title(f%22Welch's%20method%2C%20segments%20of%20%7B_N%20%2F%2F%204%7D%20samples%22)%0A%20%20%20%20_axs%5B2%5D.plot(_fw%2C%20_Pw%2C%20%22r%22%2C%20linewidth%3D1.5)%0A%20%20%20%20_axs%5B2%5D.set_xlabel(%22Frequency%20%5BHz%5D%22)%0A%20%20%20%20_axs%5B2%5D.set_ylabel(%22PSD(y)%20%5BV%24%5E2%24%2FHz%5D%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%0A%20%20%20%20print(%22Quick%20PSD%20equals%20the%20periodogram%3A%22%2C%20np.allclose(_P%2C%20_Pp))%0A%20%20%20%20print(f%22Variance%20of%20the%20data%3A%20%20%20%20%20%20%7Bnp.var(y_n)%3A.3f%7D%20V2%22)%0A%20%20%20%20print(f%22Bin-sum%20of%20periodogram%3A%20%20%20%20%7Bnp.sum(_Pp)%20*%20(_fp%5B1%5D%20-%20_fp%5B0%5D)%3A.3f%7D%20V2%22)%0A%20%20%20%20print(f%22Bin-sum%20of%20Welch's%20PSD%3A%20%20%20%20%7Bnp.sum(_Pw)%20*%20(_fw%5B1%5D%20-%20_fw%5B0%5D)%3A.3f%7D%20V2%22)%0A%20%20%20%20print(f%22Bin%20spacing%3A%20periodogram%20%7B_fp%5B1%5D%3A.2f%7D%20Hz%2C%20Welch%20%7B_fw%5B1%5D%3A.2f%7D%20Hz%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%20expected%20power%20is%20%242%5E2%2F2%20%2B%201%5E2%2F2%20%2B%201%20%3D%203.5%24%20V%C2%B2.%20The%20data%20have%20a%20variance%20of%20about%203.4%20V%C2%B2%20(this%20particular%20noise%20sample%20has%20a%20little%20less%20power%20than%20its%20nominal%20value)%2C%20and%20the%20periodogram%20bin%20sum%20equals%20that%20variance%20(its%20default%20detrending%20removes%20DC).%20Welch%20estimates%20a%20weighted%20segment%20average%20and%20need%20not%20equal%20the%20full-record%20variance%20exactly%3A%20the%20PSD%20tells%20you%20*where%20in%20frequency*%20the%20power%20of%20a%20signal%20is.%20Welch's%20estimate%20is%20smoother%20and%20its%20peaks%20wider%2C%20because%20each%20segment%20is%20a%20quarter%20of%20the%20data%20and%20has%20bins%20four%20times%20farther%20apart%2C%20as%20well%20as%20a%20wider%20main%20lobe%20from%20the%20Hann%20window.%20For%20the%20discrete%20periodogram%2C%20summing%20bin%20power%20preserves%20Parseval%20exactly%3B%20trapezoidal%20integration%20would%20give%20half%20weight%20to%20the%20endpoint%20bins.%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%20Frequency%20characteristics%20of%20a%20PSD%0A%0A%20%20%20%20A%20PSD%20is%20often%20summarized%20by%20a%20few%20frequencies%3A%0A%0A%20%20%20%20-%20the%20**peak%20frequency**%2C%20%24F_%7Bmax%7D%24%2C%20where%20the%20PSD%20is%20largest%3B%0A%20%20%20%20-%20the%20**mean%20frequency**%2C%20its%20centroid%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20F_%7Bmean%7D%20%3D%20%5Cfrac%7B%20%5Csum_%7Bi%3D1%7D%5E%7BN%7D%20F_i%5C%2CP_i%20%7D%7B%20%5Csum_%7Bi%3D1%7D%5E%7BN%7D%20P_i%20%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20-%20and%20the%20frequency%20**percentiles**%3A%20%24F_%7B50%5C%25%7D%24%2C%20the%20**median%20frequency**%2C%20below%20which%20lies%20half%20of%20the%20total%20power%2C%20and%20%24F_%7B95%5C%25%7D%24%2C%20below%20which%20lies%2095%25%20of%20it%2C%20an%20estimate%20of%20the%20bandwidth%20of%20the%20signal.%0A%0A%20%20%20%20The%20function%20%60psd2%60%2C%20from%20the%20Python%20package%20%5Bpsd2%5D(https%3A%2F%2Fpypi.org%2Fproject%2Fpsd2%2F)%20(included%20in%20this%20notebook%E2%80%99s%20dependencies)%2C%20estimates%20the%20PSD%20with%20Welch's%20method%20(it%20is%20a%20wrapper%20of%20%60scipy.signal.welch%60)%20and%20computes%20these%20characteristics%3B%20it%20returns%20the%20frequency%20percentiles%20(for%20example%2C%20%60fpcntile%5B50%5D%60%20is%20the%20median%20frequency)%2C%20the%20mean%20frequency%2C%20the%20peak%20frequency%2C%20the%20total%20power%2C%20and%20the%20PSD%20itself%2C%20and%20plots%20the%20results.%20Its%20signature%20is%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20fpcntile%2C%20mpf%2C%20fmax%2C%20Ptotal%2C%20f%2C%20P%20%3D%20psd2(x%2C%20fs%3D1.0%2C%20window%3D'hann'%2C%20nperseg%3DNone%2C%20noverlap%3DNone%2C%20nfft%3DNone%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%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20detrend%3D'constant'%2C%20show%3DTrue%2C%20ax%3DNone%2C%20scales%3D'linear'%2C%20xlim%3DNone%2C%20units%3D'V')%0A%20%20%20%20%60%60%60%0A%0A%20%20%20%20Run%20%60help(psd2)%60%20to%20read%20its%20documentation%2C%20or%20read%20its%20code%20on%20%5BGitHub%5D(https%3A%2F%2Fgithub.com%2Fdemotu%2Fpsd2).%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%20psd2%20import%20psd2%0A%0A%20%20%20%20return%20(psd2%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%20try%20it%20on%20the%20example%20from%20the%20documentation%20of%20%60scipy.signal.welch%60%3A%20a%20sine%20wave%20of%202%20V%20RMS%20at%201234%20Hz%2C%20corrupted%20by%20white%20noise%20of%200.001%20V%C2%B2%2FHz%2C%20sampled%20at%2010%20kHz%20for%2010%20s.%20**Before%20you%20run%20the%20next%20cell**%2C%20predict%20the%20peak%20frequency%20and%20the%20total%20power.%20Will%20the%20median%20and%20the%20mean%20frequency%20also%20be%20near%201234%20Hz%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20psd2)%3A%0A%20%20%20%20_fs%20%3D%2010e3%0A%20%20%20%20_N%20%3D%20100_000%0A%20%20%20%20_amp%20%3D%202%20*%20np.sqrt(2)%0A%20%20%20%20_freq%20%3D%201234.0%0A%20%20%20%20_noise_power%20%3D%200.001%20*%20_fs%20%2F%202%0A%20%20%20%20_time%20%3D%20np.arange(_N)%20%2F%20_fs%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(seed%3D42)%0A%20%20%20%20_x%20%3D%20_amp%20*%20np.sin(2%20*%20np.pi%20*%20_freq%20*%20_time)%0A%20%20%20%20_x%20%2B%3D%20_rng.normal(scale%3Dnp.sqrt(_noise_power)%2C%20size%3D_time.shape)%0A%0A%20%20%20%20_fpcntile%2C%20_mpf%2C%20_fmax%2C%20_Ptotal%2C%20_%2C%20_%20%3D%20psd2(_x%2C%20fs%3D_fs)%0A%20%20%20%20print(f%22Peak%20frequency%3A%20%20%20%7B_fmax%3A7.1f%7D%20Hz%22)%0A%20%20%20%20print(f%22Median%20frequency%3A%20%7B_fpcntile%5B50%5D%3A7.1f%7D%20Hz%22)%0A%20%20%20%20print(f%22Mean%20frequency%3A%20%20%20%7B_mpf%3A7.1f%7D%20Hz%22)%0A%20%20%20%20print(f%22Total%20power%3A%20%20%20%20%20%20%7B_Ptotal%3A7.2f%7D%20V2%20(variance%20of%20the%20data%3A%20%7Bnp.var(_x)%3A.2f%7D%20V2)%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%20peak%20and%20the%20median%20are%20at%201234%20Hz%2C%20and%20the%20total%20power%20is%20about%209%20V%C2%B2%3A%204%20V%C2%B2%20from%20the%20sine%20(2%20V%20RMS%2C%20squared)%20plus%205%20V%C2%B2%20from%20the%20noise%20(0.001%20V%C2%B2%2FHz%20over%205000%20Hz).%20The%20mean%20frequency%2C%20though%2C%20is%20near%201940%20Hz%2C%20where%20there%20is%20nothing%20but%20noise.%20The%20noise%20floor%20looks%20negligible%20next%20to%20the%20peak%2C%20but%20it%20extends%20over%205000%20Hz%2C%20and%20the%20mean%20frequency%20weighs%20every%20hertz%20of%20it.%20The%20median%20is%20much%20less%20sensitive.%20Choose%20the%20summary%20frequency%20with%20the%20noise%20in%20mind.%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%20Back%20to%20the%20muscle%0A%0A%20%20%20%20Now%20we%20can%20answer%20the%20questions%20from%20the%20beginning.%20Let's%20estimate%20the%20PSD%20of%20the%20EMG%2C%20with%20Welch's%20method%20and%20segments%20of%20256%20samples%20(a%20resolution%20of%20about%204%20Hz).%20**Before%20you%20run%20the%20next%20cell**%2C%20go%20back%20to%20your%20answer%20to%20Guiding%20question%200.1%20and%20predict%20the%20median%20frequency%20of%20this%20EMG.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(emg%2C%20emg_fs%2C%20psd2)%3A%0A%20%20%20%20_fpcntile%2C%20_mpf%2C%20_fmax%2C%20_Ptotal%2C%20_%2C%20_%20%3D%20psd2(emg%2C%20fs%3Demg_fs%2C%20nperseg%3D256)%0A%20%20%20%20print(f%22Peak%20frequency%3A%20%20%20%7B_fmax%3A6.1f%7D%20Hz%22)%0A%20%20%20%20print(f%22Median%20frequency%3A%20%7B_fpcntile%5B50%5D%3A6.1f%7D%20Hz%22)%0A%20%20%20%20print(f%22Mean%20frequency%3A%20%20%20%7B_mpf%3A6.1f%7D%20Hz%22)%0A%20%20%20%20print(f%2295%25%20of%20the%20power%20below%20%7B_fpcntile%5B95%5D%3A.1f%7D%20Hz%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%20power%20of%20this%20EMG%20is%20spread%20over%20a%20broad%20band%2C%20from%20a%20few%20tens%20of%20hertz%20to%20about%20400%20Hz%2C%20with%20a%20median%20frequency%20of%20about%20164%20Hz%20and%20a%20mean%20frequency%20a%20little%20higher.%20No%20single%20frequency%20dominates%3A%20an%20EMG%20is%20the%20sum%20of%20the%20electrical%20activity%20of%20many%20motor%20units%2C%20each%20firing%20irregularly%2C%20and%20its%20spectrum%20is%20broad.%0A%0A%20%20%20%20That%20breadth%20is%20what%20makes%20the%20median%20frequency%20useful.%20During%20a%20sustained%20contraction%2C%20as%20the%20muscle%20fatigues%2C%20the%20conduction%20velocity%20of%20its%20fibres%20decreases%20and%20the%20whole%20spectrum%20shifts%20to%20lower%20frequencies%3B%20the%20median%20frequency%20of%20successive%20windows%20of%20the%20EMG%20falls%2C%20and%20this%20is%20one%20commonly%20used%20indicator%20of%20muscle%20fatigue%20(De%20Luca%2C%201997).%20Electrode%20geometry%2C%20recruitment%2C%20motion%20and%20analysis%20choices%20can%20also%20change%20the%20spectrum%2C%20so%20a%20median-frequency%20change%20alone%20does%20not%20establish%20fatigue.%20Here%2095%25%20of%20the%20recorded%20power%20is%20below%20400%20Hz.%20That%20describes%20the%20stored%20signal%3B%20it%20cannot%20prove%20that%20an%20analog%20anti-alias%20filter%20removed%20frequencies%20above%20500%20Hz%20before%20acquisition.%0A%0A%20%20%20%20**Guiding%20questions%202.**%0A%0A%20%20%20%201.%20How%20close%20was%20your%20estimate%20from%20counting%20zero%20crossings%20to%20the%20median%20frequency%3F%20Why%20might%20the%20two%20differ%3F%0A%20%20%20%202.%20Compute%20the%20PSD%20again%20with%20%60nperseg%3D64%60%20and%20with%20%60nperseg%3D1024%60.%20How%20do%20the%20curve%20and%20the%20median%20frequency%20change%3F%20Which%20would%20you%20report%3F%0A%20%20%20%203.%20This%20recording%20has%20three%20bursts%20of%20activity.%20How%20would%20you%20check%20whether%20the%20median%20frequency%20changes%20from%20one%20burst%20to%20the%20next%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%205.**%20You%20have%202%20s%20of%20EMG%20sampled%20at%201000%20Hz.%20Compare%20Welch%20windows%20of%20250%20and%201000%20samples%2C%20and%20explain%20why%20a%20zero-crossing%20estimate%20need%20not%20equal%20the%20median%20frequency.%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%205%0A%20%20%20%20The%20segment%20lengths%20are%200.25%20and%201%20s%2C%20giving%20bin%20spacings%204%20and%201%20Hz.%20Short%20windows%20permit%20more%20averages%20and%20track%20shorter%20changes%2C%20but%20broaden%20spectral%20features.%20Median%20frequency%20divides%20spectral%20power%20into%20equal%20halves%3B%20zero%20crossings%20depend%20on%20the%20waveform%20and%20noise%20and%20are%20not%20that%20power-weighted%20statistic.%20Use%20equal%20windowing%20and%20preprocessing%20when%20comparing%20bursts.%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%20Short-time%20Fourier%20transform%0A%0A%20%20%20%20The%20spectrum%20of%20a%20whole%20recording%20averages%20over%20its%20whole%20duration%2C%20and%20the%20last%20guiding%20question%20needs%20something%20else%3A%20how%20the%20frequency%20content%20*changes%20in%20time*.%20The%20short-time%20Fourier%20transform%20answers%20it%20by%20computing%20the%20FFT%20of%20short%2C%20overlapping%20windows%20of%20the%20data%2C%20one%20after%20the%20other.%20The%20result%2C%20plotted%20as%20a%20colour%20map%20of%20power%20against%20time%20and%20frequency%2C%20is%20a%20**spectrogram**.%0A%0A%20%20%20%20A%20simple%20test%20is%20a%20%5Bchirp%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FChirp)%2C%20a%20sinusoid%20whose%20frequency%20rises%20linearly%2C%20here%20from%20100%20Hz%20to%20300%20Hz%20in%20about%204%20s%2C%20sampled%20at%201000%20Hz.%20**Before%20you%20run%20the%20next%20cell**%2C%20sketch%20what%20its%20spectrogram%20should%20look%20like.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20signal)%3A%0A%20%20%20%20_fs%20%3D%201000%0A%20%20%20%20_t%20%3D%20np.arange(2**12)%20%2F%20_fs%0A%20%20%20%20_c%20%3D%20signal.chirp(_t%2C%20f0%3D100%2C%20f1%3D300%2C%20t1%3D_t%5B-1%5D%2C%20method%3D%22linear%22)%0A%0A%20%20%20%20_%2C%20_ax%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(10%2C%204.5))%0A%20%20%20%20_ax.specgram(_c%2C%20NFFT%3D256%2C%20Fs%3D_fs%2C%20noverlap%3D128%2C%20cmap%3D%22gist_heat%22)%0A%20%20%20%20_ax.set_title(%22Spectrogram%20of%20a%20chirp%22)%0A%20%20%20%20_ax.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_ax.set_ylabel(%22Frequency%20%5BHz%5D%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20A%20straight%20line%20rising%20from%20100%20to%20300%20Hz%2C%20as%20it%20should%20be.%20Each%20column%20is%20the%20spectrum%20of%20256%20samples%2C%20a%20quarter%20of%20a%20second%2C%20so%20the%20spectrogram%20can%20follow%20changes%20that%20are%20slow%20compared%20with%20a%20quarter%20of%20a%20second%2C%20and%20place%20DFT%20bins%20about%204%20Hz%20apart%20(the%20window%20main%20lobe%20is%20wider).%20A%20shorter%20window%20would%20follow%20faster%20changes%20and%20blur%20the%20frequencies%20more%3B%20a%20longer%20one%2C%20the%20reverse.%20This%20is%20the%20time-frequency%20trade-off%20of%20every%20spectrogram%2C%20the%20same%20%24%5CDelta%20f%20%3D%201%2F(N%5CDelta%20t)%24%20applied%20to%20each%20window.%0A%0A%20%20%20%20Now%20the%20EMG%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(emg%2C%20emg_fs%2C%20np%2C%20plt)%3A%0A%20%20%20%20_fig%2C%20_axs%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(10%2C%206)%2C%20sharex%3DTrue)%0A%20%20%20%20_axs%5B0%5D.plot(np.arange(emg.size)%20%2F%20emg_fs%2C%20emg%2C%20linewidth%3D0.8)%0A%20%20%20%20_axs%5B0%5D.set_ylabel(%22EMG%20%5BV%5D%22)%0A%20%20%20%20_axs%5B1%5D.specgram(emg%2C%20NFFT%3D128%2C%20Fs%3Demg_fs%2C%20noverlap%3D64%2C%20cmap%3D%22gist_heat%22)%0A%20%20%20%20_axs%5B1%5D.set_xlim(0%2C%20emg.size%20%2F%20emg_fs)%0A%20%20%20%20_axs%5B1%5D.set_ylim(0%2C%20500)%0A%20%20%20%20_axs%5B1%5D.set_xlabel(%22Time%20%5Bs%5D%22)%0A%20%20%20%20_axs%5B1%5D.set_ylabel(%22Frequency%20%5BHz%5D%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Each%20burst%20of%20activity%20is%20a%20vertical%20band%20of%20power%20spread%20over%20a%20broad%20range%20of%20frequencies%2C%20brightest%20roughly%20between%2050%20and%20250%20Hz%3B%20between%20the%20bursts%2C%20there%20is%20almost%20nothing.%20To%20follow%20the%20median%20frequency%20during%20a%20fatiguing%20contraction%2C%20you%20would%20compute%20it%20for%20each%20column%20of%20a%20spectrogram%20like%20this%20one%2C%20or%20for%20each%20window%20of%20a%20few%20hundred%20milliseconds%20of%20the%20EMG.%0A%0A%20%20%20%20**Challenge%202.**%20Use%20%60psd2%60%20on%20each%20of%20the%20three%20bursts%20of%20the%20EMG%20separately%20(they%20are%20roughly%20centred%20at%200.5%2C%201.5%20and%202.6%20s%3B%20take%200.5%20s%20around%20each).%20Does%20the%20median%20frequency%20change%20from%20burst%20to%20burst%3F%20Is%20the%20difference%20larger%20than%20what%20you%20get%20by%20moving%20each%20window%20by%200.1%20s%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%206.**%20For%20a%20256-sample%20window%20at%201000%20Hz%2C%20compute%20its%20duration%20and%20bin%20spacing.%20How%20does%20halving%20the%20window%20affect%20a%20changing-frequency%20signal%3F%0A%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Solution%20to%20Exercise%206%0A%20%20%20%20The%20duration%20is%20%24256%2F1000%3D0.256%24%20s%20and%20the%20bin%20spacing%20is%20%241000%2F256%3D3.90625%24%20Hz.%20A%20128-sample%20window%20spans%200.128%20s%20with%207.8125%20Hz%20bins.%20It%20localizes%20changes%20more%20closely%20in%20time%2C%20while%20spreading%20their%20frequency%20representation%20more%20broadly.%0A%20%20%20%20%2F%2F%2F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Checkpoint%20questions%0A%0A%20%20%20%20Pause%20here%20before%20the%20problems.%0A%0A%20%20%20%201.%20Go%20back%20to%20your%20guess%20in%20Challenge%200.%20How%20long%20would%20you%20need%20to%20record%20that%20signal%20to%20obtain%20DFT%20bins%200.5%20Hz%20apart%2C%20and%20how%20fast%20would%20you%20need%20to%20sample%20it%3F%0A%20%20%20%202.%20The%20FFT%20of%20a%203%20Hz%20sine%20sampled%20for%201%20s%20at%20100%20Hz%20has%20a%20clean%20single%20peak%3B%20sampled%20for%201.1%20s%2C%20it%20does%20not.%20Why%3F%0A%20%20%20%203.%20An%20FFT%20reports%20a%20phase%20of%20%24-90%5Eo%24%20for%20a%20component.%20Is%20that%20component%20a%20sine%20or%20a%20cosine%3F%0A%20%20%20%204.%20Why%20is%20the%20amplitude%20of%20the%20one-sided%20spectrum%20doubled%2C%20except%20at%20zero%20frequency%20and%20(for%20even-length%20data)%20Nyquist%3F%0A%20%20%20%205.%20Why%20is%20the%20periodogram%20of%20a%20noisy%20signal%20so%20noisy%2C%20and%20what%20does%20Welch's%20method%20give%20up%20to%20make%20it%20smoother%3F%0A%20%20%20%206.%20Two%20EMG%20recordings%20have%20the%20same%20median%20frequency%20but%20different%20mean%20frequencies.%20What%20might%20be%20different%20about%20them%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%2F%2F%2F%20details%20%7C%20Checkpoint%20answers%0A%20%20%20%201.%20Two%20seconds%20gives%200.5%20Hz%20bin%20spacing.%20Choose%20%24f_s%3E2f_%7B%5Cmax%7D%24%20with%20an%20analog%20anti-alias%20filter%20and%20transition%20margin%3B%20resolving%20close%20noisy%20tones%20may%20need%20a%20longer%20record.%0A%20%20%20%202.%20Three%20cycles%20fit%201%20s%3B%203.3%20cycles%20do%20not%20fit%201.1%20s%2C%20producing%20leakage%20on%20that%20DFT%20grid.%0A%20%20%20%203.%20A%20positive%20sine%2C%20since%20%24%5Ccos(%5Comega%20t-%5Cpi%2F2)%3D%5Csin(%5Comega%20t)%24.%0A%20%20%20%204.%20Positive%20and%20negative%20members%20of%20a%20conjugate%20pair%20carry%20equal%20amplitudes.%20DC%20and%20the%20even-length%20Nyquist%20bin%20are%20unpaired.%0A%20%20%20%205.%20A%20periodogram%20has%20high%20pointwise%20variability.%20Welch%20averages%20segments%2C%20sacrificing%20segment%20duration%20and%20frequency%20detail.%0A%20%20%20%206.%20Their%20power%20distributions%20can%20have%20different%20tails%20or%20asymmetry%20while%20retaining%20the%20same%20halfway%20point.%20This%20alone%20does%20not%20identify%20a%20physiological%20cause.%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.%20Compute%20the%20magnitude%20of%20the%20Fourier-transform%20approximation%20%60dt%20*%20np.fft.fft(pulse)%60%20for%20a%20rectangular%20pulse%20of%20unit%20height%20and%200.5%20s%20duration%2C%20in%20the%20middle%20of%2010%20s%20of%20zeros%20sampled%20at%20100%20Hz.%20Compare%20it%20with%20the%20magnitude%20of%20the%20analytical%20Fourier%20transform%20of%20a%20rectangular%20pulse%2C%20%24%7CT%5C%2C%5Cmathrm%7Bsinc%7D(fT)%7C%24%2C%20with%20%24T%20%3D%200.5%24%20s.%0A%0A%20%20%20%202.%20Check%20Parseval's%20theorem%20numerically%20for%20the%20noisy%20signal%20above%3A%20%24%5Csum_n%20%7Cx%5Bn%5D%7C%5E2%20%3D%20%5Cfrac%7B1%7D%7BN%7D%5Csum_k%20%7CX%5Bk%5D%7C%5E2%24.%0A%0A%20%20%20%203.%20Sample%20a%2070%20Hz%20sine%20at%20100%20Hz%20and%20compute%20its%20spectrum.%20At%20what%20frequency%20does%20the%20peak%20appear%3F%20Relate%20it%20to%20the%20section%20on%20aliasing%20in%20%5BBasic%20properties%20of%20signals%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fsignal_basic_properties.py).%0A%0A%20%20%20%204.%20Load%20Pezzack's%20benchmark%20angle%20from%20%5BData%20filtering%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fdata_filtering.py)%20and%20plot%20its%20PSD%20on%20a%20logarithmic%20scale%20(%60scales%3D'semilogy'%60).%20Above%20what%20frequency%20does%20it%20flatten%20into%20a%20noise%20floor%3F%20Compare%20with%20the%20optimal%20cutoff%20frequency%20found%20in%20%5BResidual%20analysis%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fresidual_analysis.py).%0A%0A%20%20%20%205.%20Zero%20padding%3A%20compute%20the%20spectrum%20of%20the%205.1%20Hz%20sine%20with%20%60np.fft.rfft(x%2C%20n%3D4*N)%60.%20Does%20the%20peak%20get%20closer%20to%202%3F%20Does%20the%20frequency%20resolution%20really%20improve%3F%0A%0A%20%20%20%206.%20Compute%20the%20spectrogram%20of%20the%20noisy%20two-sine%20signal%20with%20windows%20of%2032%2C%20128%20and%20%60min(512%2C%20len(y_n))%60%20samples%20(the%20record%20contains%20only%20500%20samples)%2C%20and%20describe%20the%20time-frequency%20trade-off%20you%20see.%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.%20The%20zero-frequency%20magnitude%20is%200.5%20s%3B%20the%20first%20nonzero%20nulls%20are%20at%20%24%5Cpm2%24%20Hz%2C%20then%20multiples%20of%202%20Hz.%20%60dt%20*%20fft(pulse)%60%20approximates%20the%20continuous%20transform%20%24T%5C%2C%5Cmathrm%7Bsinc%7D(fT)%24%20in%20magnitude%2C%20with%20discrete-time%20deviations%20away%20from%20low%20frequencies.%20A%20one-sided%20amplitude%20spectrum%20scaled%20by%20%242%2FN%24%20has%20different%20normalization%20and%20must%20not%20be%20compared%20directly%20with%20this%20transform.%0A%0A%20%20%20%202.%20Both%20sides%20are%20approximately%20%241698.51629126209%24%20for%20the%20supplied%20seed%20and%20record.%20Agreement%20is%20to%20floating-point%20precision.%0A%0A%20%20%20%203.%20The%20peak%20appears%20at%2030%20Hz.%20At%20these%20sample%20times%2C%20%24%5Csin(2%5Cpi70t)%3D-%5Csin(2%5Cpi30t)%24%3B%20amplitude%20alone%20cannot%20distinguish%20them.%0A%0A%20%20%20%204.%20The%20low-power%20tail%20is%20broadly%20visible%20around%2010%E2%80%9320%20Hz%2C%20with%20its%20appearance%20depending%20on%20the%20window%20and%20segment%20length.%20There%20is%20no%20unique%20sharp%20boundary.%20The%20residual-analysis%20choice%20is%20about%205.59%20Hz%20under%20its%20stated%20filter%20convention%2C%20below%20this%20tail.%20A%20visual%20PSD%20boundary%20and%20a%20residual%20criterion%20need%20not%20agree.%0A%0A%20%20%20%205.%20With%20normalization%20by%20the%20original%20%24N%3D500%24%2C%20the%20padded%20grid%20contains%205.1%20Hz%20and%20the%20peak%20is%202.000.%20The%20grid%20spacing%20becomes%200.05%20Hz%3B%20the%20information%20and%20main-lobe%20width%20of%20the%20original%205%20s%20record%20do%20not%20improve.%20Dividing%20by%20%244N%24%20would%20incorrectly%20reduce%20amplitudes%20fourfold.%0A%0A%20%20%20%206.%20The%2032-%2C%20128-%20and%20500-sample%20windows%20span%200.32%2C%201.28%20and%205%20s%2C%20with%20bin%20spacings%203.125%2C%200.78125%20and%200.2%20Hz.%20Short%20windows%20show%20fluctuations%20in%20time%20more%20clearly%3B%20long%20windows%20separate%20frequency%20components%20more%20clearly.%20A%20500-sample%20window%20gives%20only%20one%20full-record%20segment%20here.%20A%20literal%20request%20for%20512%20samples%20exceeds%20the%20available%20record.%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%5BFourier%20series%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Ffourier_series.py)%20%E2%80%94%20the%20periodic%20case%2C%20harmonics%20and%20their%20coefficients.%0A%20%20%20%20-%20%5BBasic%20properties%20of%20signals%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Fsignal_basic_properties.py)%20%E2%80%94%20sampling%2C%20aliasing%20and%20power.%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%20filters%20described%20by%20their%20frequency%20response.%0A%20%20%20%20-%20%5BElectromyography%5D(https%3A%2F%2Fmolab.marimo.io%2Fgithub%2FBMClab%2FBMC%2Fblob%2Fmaster%2Fnotebooks%2Felectromyography.py)%20%E2%80%94%20processing%20EMG%20in%20time%20and%20in%20frequency.%0A%20%20%20%20-%20Smith%20SW%20(1997)%20%5BThe%20Scientist%20and%20Engineer's%20Guide%20to%20Digital%20Signal%20Processing%5D(https%3A%2F%2Fwww.dspguide.com%2F)%2C%20chapters%208%20to%2012%20on%20the%20DFT%20and%20the%20FFT.%0A%20%20%20%20%23%23%23%20Video%20lectures%20on%20the%20Internet%0A%0A%20%20%20%20-%203Blue1Brown%3A%20%5BBut%20what%20is%20the%20Fourier%20Transform%3F%20A%20visual%20introduction%5D(https%3A%2F%2Fwww.youtube.com%2Fwatch%3Fv%3DspUNpyF58BY).%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-%20Bartlett%20MS%20(1950)%20Periodogram%20analysis%20and%20continuous%20spectra.%20Biometrika%2C%2037%2C%201-16.%0A%0A%20%20%20%20-%20De%20Luca%20CJ%20(1997)%20The%20use%20of%20surface%20electromyography%20in%20biomechanics.%20Journal%20of%20Applied%20Biomechanics%2C%2013%2C%20135-163.%0A%0A%20%20%20%20-%20Lyons%20RG%20(2010)%20%5BUnderstanding%20Digital%20Signal%20Processing%5D(https%3A%2F%2Fbooks.google.com.br%2Fbooks%3Fid%3DUBU7Y2tpwWUC%26hl).%203rd%20edition.%20Prentice%20Hall.%0A%0A%20%20%20%20-%20Smith%20SW%20(1997)%20%5BThe%20Scientist%20and%20Engineer's%20Guide%20to%20Digital%20Signal%20Processing%5D(https%3A%2F%2Fwww.dspguide.com%2F).%20California%20Technical%20Pub.%0A%0A%20%20%20%20-%20Welch%20P%20(1967)%20The%20use%20of%20the%20fast%20Fourier%20transform%20for%20the%20estimation%20of%20power%20spectra%3A%20a%20method%20based%20on%20time%20averaging%20over%20short%2C%20modified%20periodograms.%20IEEE%20Transactions%20on%20Audio%20and%20Electroacoustics%2C%2015%2C%2070-73.%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
4a85b7ba50d8a6f5886445b61761627d