beamforming mvdr interactive

%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%22matplotlib%3D%3D3.10.8%22%2C%20%22numpy%3D%3D2.4.3%22%5D%0A%23%20%0A%23%20%5Btool.marimo.export%5D%0A%23%20lock_kind%20%3D%20%22resolved%22%0A%23%20%2F%2F%2F%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.25.1%22%0Aapp%20%3D%20marimo.App()%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20plt.style.use(%22default%22)%0A%20%20%20%20BLUE%2C%20ORANGE%2C%20RED%2C%20GREY%2C%20GREEN%20%3D%20%22%231f77b4%22%2C%20%22%23ff7f0e%22%2C%20%22%23d62728%22%2C%20%220.55%22%2C%20%22%232ca02c%22%0A%20%20%20%20return%20BLUE%2C%20GREY%2C%20ORANGE%2C%20RED%2C%20mo%2C%20np%2C%20plt%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Listening%20in%20One%20Direction%3A%20Delay-and-Sum%20and%20MVDR%0A%0A%20%20%20%20A%20uniform%20line%20of%20%24M%24%20sensors%2C%20half%20a%20wavelength%20apart%2C%20hears%20a%20narrowband%20**target**%20from%20the%20look%0A%20%20%20%20direction%20%24%5Ctheta_0%24%20and%20a%20strong%20**interferer**%20from%20somewhere%20else%2C%20in%20white%20sensor%20noise.%20A%20beamformer%0A%20%20%20%20combines%20the%20sensors%2C%20%24%5Cwidehat%20Y%3D%5Cmathbf%20w%5E%5Cmathsf%7BH%7D%5Cmathbf%20X%24%2C%20and%20its%20response%20to%20a%20plane%20wave%20from%0A%20%20%20%20direction%20%24%5Ctheta%24%20is%20%24%5Clvert%5Cmathbf%20w%5E%5Cmathsf%7BH%7D%5Cmathbf%20a(%5Ctheta)%5Crvert%5E2%24.%0A%0A%20%20%20%20-%20**Delay-and-sum**%2C%20%24%5Cmathbf%20w%3D%5Cmathbf%20a(%5Ctheta_0)%2FM%24%3A%20the%20matched%20filter%20for%20the%20look%20direction%20in%0A%20%20%20%20%20%20white%20noise.%20It%20knows%20nothing%20about%20the%20interferer.%0A%20%20%20%20-%20**MVDR**%2C%20%24%5Cmathbf%20w%3D%5Cmathbf%20R_%7B%5Cmathbf%7BXX%7D%7D%5E%7B-1%7D%5Cmathbf%20a(%5Ctheta_0)%5C%2C%2F%5C%2C%5Cmathbf%20a%5E%5Cmathsf%7BH%7D(%5Ctheta_0)%5Cmathbf%20R_%7B%5Cmathbf%7BXX%7D%7D%5E%7B-1%7D%5Cmathbf%20a(%5Ctheta_0)%24%3A%0A%20%20%20%20%20%20unit%20gain%20toward%20%24%5Ctheta_0%24%2C%20and%20as%20little%20output%20power%20as%20possible%20otherwise.%20It%20learns%20where%20the%0A%20%20%20%20%20%20interferer%20is%20from%20the%20data%20covariance%20alone.%0A%0A%20%20%20%20Tick%20**estimate%20%24%5Cmathbf%20R_%7B%5Cmathbf%7BXX%7D%7D%24%20from%20snapshots**%20to%20replace%20the%20true%20covariance%20by%20a%20sample%0A%20%20%20%20average%2C%20and%20watch%20what%20a%20short%20record%20does%20to%20MVDR.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20M%20%3D%20mo.ui.slider(2%2C%2032%2C%20step%3D1%2C%20value%3D8%2C%20label%3D%22sensors%20M%22)%0A%20%20%20%20look%20%3D%20mo.ui.slider(-60%2C%2060%2C%20step%3D1%2C%20value%3D0%2C%20label%3D%22look%20direction%20%CE%B8%E2%82%80%20(%C2%B0)%22)%0A%20%20%20%20intf%20%3D%20mo.ui.slider(-89%2C%2089%2C%20step%3D1%2C%20value%3D-21%2C%20label%3D%22interferer%20direction%20(%C2%B0)%22)%0A%20%20%20%20inr%20%3D%20mo.ui.slider(-10%2C%2040%2C%20step%3D1%2C%20value%3D20%2C%20label%3D%22interferer-to-noise%20ratio%20(dB)%22)%0A%20%20%20%20snr%20%3D%20mo.ui.slider(-20%2C%2020%2C%20step%3D1%2C%20value%3D0%2C%20label%3D%22target-to-noise%20ratio%20(dB)%22)%0A%20%20%20%20sample%20%3D%20mo.ui.checkbox(label%3D%22estimate%20R_XX%20from%20snapshots%22%2C%20value%3DFalse)%0A%20%20%20%20L%20%3D%20mo.ui.slider(4%2C%20400%2C%20step%3D2%2C%20value%3D40%2C%20label%3D%22snapshots%20L%22)%0A%20%20%20%20mo.vstack(%5Bmo.hstack(%5BM%2C%20look%2C%20intf%5D)%2C%20mo.hstack(%5Binr%2C%20snr%5D)%2C%20mo.hstack(%5Bsample%2C%20L%5D)%5D)%0A%20%20%20%20return%20L%2C%20M%2C%20inr%2C%20intf%2C%20look%2C%20sample%2C%20snr%0A%0A%0A%40app.cell%0Adef%20_(L%2C%20M%2C%20inr%2C%20intf%2C%20look%2C%20np%2C%20sample%2C%20snr)%3A%0A%20%20%20%20m%20%3D%20np.arange(M.value)%0A%0A%20%20%20%20def%20steer(theta_deg)%3A%0A%20%20%20%20%20%20%20%20return%20np.exp(-1j%20*%20np.pi%20*%20m%20*%20np.sin(np.deg2rad(theta_deg)))%0A%0A%20%20%20%20a0%2C%20ai%20%3D%20steer(look.value)%2C%20steer(intf.value)%0A%20%20%20%20p_s%2C%20p_i%20%3D%2010%20**%20(snr.value%20%2F%2010)%2C%2010%20**%20(inr.value%20%2F%2010)%0A%20%20%20%20R_true%20%3D%20p_s%20*%20np.outer(a0%2C%20a0.conj())%20%2B%20p_i%20*%20np.outer(ai%2C%20ai.conj())%20%2B%20np.eye(M.value)%0A%20%20%20%20if%20sample.value%3A%0A%20%20%20%20%20%20%20%20_rng%20%3D%20np.random.default_rng(7)%0A%20%20%20%20%20%20%20%20_cn%20%3D%20lambda%20*s%3A%20(_rng.standard_normal(s)%20%2B%201j%20*%20_rng.standard_normal(s))%20%2F%20np.sqrt(2)%0A%20%20%20%20%20%20%20%20_X%20%3D%20(np.sqrt(p_s)%20*%20np.outer(a0%2C%20_cn(L.value))%20%2B%20np.sqrt(p_i)%20*%20np.outer(ai%2C%20_cn(L.value))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%2B%20_cn(M.value%2C%20L.value))%0A%20%20%20%20%20%20%20%20R%20%3D%20_X%20%40%20_X.conj().T%20%2F%20L.value%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20R%20%3D%20R_true%0A%20%20%20%20w_ds%20%3D%20a0%20%2F%20M.value%0A%20%20%20%20_Ra%20%3D%20np.linalg.solve(R%2C%20a0)%0A%20%20%20%20w_mv%20%3D%20_Ra%20%2F%20(a0.conj()%20%40%20_Ra)%0A%0A%20%20%20%20def%20sinr(w)%3A%0A%20%20%20%20%20%20%20%20sig%20%3D%20p_s%20*%20abs(w.conj()%20%40%20a0)%20**%202%0A%20%20%20%20%20%20%20%20rest%20%3D%20(w.conj()%20%40%20(R_true%20-%20p_s%20*%20np.outer(a0%2C%20a0.conj()))%20%40%20w).real%0A%20%20%20%20%20%20%20%20return%2010%20*%20np.log10(sig%20%2F%20rest)%0A%0A%20%20%20%20th%20%3D%20np.linspace(-90%2C%2090%2C%201801)%0A%20%20%20%20A%20%3D%20np.exp(-1j%20*%20np.pi%20*%20np.outer(m%2C%20np.sin(np.deg2rad(th))))%0A%20%20%20%20return%20A%2C%20R%2C%20sinr%2C%20th%2C%20w_ds%2C%20w_mv%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20A%2C%0A%20%20%20%20BLUE%2C%0A%20%20%20%20GREY%2C%0A%20%20%20%20ORANGE%2C%0A%20%20%20%20R%2C%0A%20%20%20%20RED%2C%0A%20%20%20%20intf%2C%0A%20%20%20%20look%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20sinr%2C%0A%20%20%20%20th%2C%0A%20%20%20%20w_ds%2C%0A%20%20%20%20w_mv%2C%0A)%3A%0A%20%20%20%20def%20patt(w)%3A%0A%20%20%20%20%20%20%20%20p%20%3D%20np.abs(w.conj()%20%40%20A)%20**%202%0A%20%20%20%20%20%20%20%20return%2010%20*%20np.log10(np.maximum(p%2C%201e-8))%0A%0A%20%20%20%20_cap%20%3D%20np.real(1%20%2F%20np.einsum(%22mt%2Cmn%2Cnt-%3Et%22%2C%20A.conj()%2C%20np.linalg.inv(R)%2C%20A))%0A%20%20%20%20_bart%20%3D%20np.real(np.einsum(%22mt%2Cmn%2Cnt-%3Et%22%2C%20A.conj()%2C%20R%2C%20A))%20%2F%20A.shape%5B0%5D%20**%202%0A%0A%20%20%20%20fig%2C%20(a1%2C%20a2)%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(12%2C%207.2)%2C%20sharex%3DTrue)%0A%20%20%20%20a1.plot(th%2C%20patt(w_ds)%2C%20color%3DGREY%2C%20lw%3D2%2C%20label%3Df%22delay-and-sum%20%20%20(output%20SINR%20%7Bsinr(w_ds)%3A.1f%7D%20dB)%22)%0A%20%20%20%20a1.plot(th%2C%20patt(w_mv)%2C%20color%3DBLUE%2C%20lw%3D2.2%2C%20label%3Df%22MVDR%20%20%20(output%20SINR%20%7Bsinr(w_mv)%3A.1f%7D%20dB)%22)%0A%20%20%20%20a1.axvline(look.value%2C%20color%3DORANGE%2C%20ls%3D%22--%22%2C%20lw%3D1.6)%0A%20%20%20%20a1.axvline(intf.value%2C%20color%3DRED%2C%20ls%3D%22--%22%2C%20lw%3D1.6)%0A%20%20%20%20a1.set_ylim(-60%2C%2010)%3B%20a1.set_ylabel(%22response%20(dB)%22)%0A%20%20%20%20a1.legend(fontsize%3D11%2C%20frameon%3DFalse%2C%20loc%3D%22lower%20right%22)%0A%20%20%20%20a1.set_title(%22beampatterns%20(orange%3A%20look%20direction%2C%20red%3A%20interferer)%22%2C%20fontsize%3D12)%0A%0A%20%20%20%20a2.plot(th%2C%2010%20*%20np.log10(_bart)%2C%20color%3DGREY%2C%20lw%3D2%2C%20label%3D%22delay-and-sum%20scan%20(Bartlett)%22)%0A%20%20%20%20a2.plot(th%2C%2010%20*%20np.log10(_cap)%2C%20color%3DBLUE%2C%20lw%3D2.2%2C%20label%3D%22MVDR%20scan%20(Capon)%22)%0A%20%20%20%20a2.axvline(look.value%2C%20color%3DORANGE%2C%20ls%3D%22--%22%2C%20lw%3D1.6)%0A%20%20%20%20a2.axvline(intf.value%2C%20color%3DRED%2C%20ls%3D%22--%22%2C%20lw%3D1.6)%0A%20%20%20%20a2.set_xlabel(r%22direction%20%24%5Ctheta%24%20(degrees)%22)%3B%20a2.set_ylabel(%22power%20(dB)%22)%0A%20%20%20%20a2.legend(fontsize%3D11%2C%20frameon%3DFalse%2C%20loc%3D%22upper%20right%22)%0A%20%20%20%20a2.set_title(%22scanning%20the%20look%20direction%3A%20where%20is%20the%20power%3F%22%2C%20fontsize%3D12)%0A%20%20%20%20fig.tight_layout()%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22Top%3A%20each%20beamformer's%20response%20toward%20every%20direction.%20Bottom%3A%20steer%20each%20beamformer%20over%20all%20directions%20and%20record%20its%20output%20power%3B%20the%20peaks%20locate%20the%20sources.%22)%2C%0A%20%20%20%20%20%20%20%20fig%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20What%20to%20try%0A%0A%20%20%20%20-%20**Move%20the%20interferer%20onto%20a%20sidelobe%20of%20delay-and-sum.**%20MVDR%20puts%20a%20null%20exactly%20on%20it%2C%20and%20its%20output%20SINR%20stays%20high%20while%20delay-and-sum's%20collapses.%0A%20%20%20%20-%20**Raise%20the%20interferer%20power.**%20The%20MVDR%20null%20gets%20deeper%3A%20the%20stronger%20the%20interferer%2C%20the%20more%20it%20dominates%20%24%5Cmathbf%20R_%7B%5Cmathbf%7BXX%7D%7D%24%2C%20and%20the%20more%20%24%5Cmathbf%20R_%7B%5Cmathbf%7BXX%7D%7D%5E%7B-1%7D%24%20suppresses%20that%20direction.%0A%20%20%20%20-%20**Move%20the%20interferer%20close%20to%20the%20look%20direction.**%20Once%20it%20is%20inside%20the%20main%20lobe%2C%20nulling%20it%20also%20costs%20target%20gain%3B%20resolution%20is%20set%20by%20the%20aperture%20%24M%24.%0A%20%20%20%20-%20**Estimate%20%24%5Cmathbf%20R_%7B%5Cmathbf%7BXX%7D%7D%24%20from%20few%20snapshots%2C%20with%20a%20strong%20target.**%20MVDR%20starts%20to%20cancel%20the%20target%20itself%3A%20the%20sample%20covariance%20is%20a%20noisy%20estimate%2C%20and%20a%20slightly%20wrong%20steering%20or%20covariance%20lets%20the%20target%20look%20like%20interference.%20More%20snapshots%2C%20or%20diagonal%20loading%2C%20cure%20it.%0A%20%20%20%20-%20**Compare%20the%20two%20scans.**%20Capon's%20peaks%20are%20much%20sharper%20than%20Bartlett's%20for%20the%20same%20array%3A%20the%20same%20%24%5Cmathbf%20R%5E%7B-1%7D%24%20that%20nulls%20also%20resolves.%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