%23%20gallery_category%3A%20FDN%20Design%20%26%20Analysis%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.23.13%22%0Aapp%20%3D%20marimo.App()%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%20Filter%20feedback%20delay%20network%20(FFDN)%20with%20paraunitary%20feedback%20matrix%0A%0A%20%20%20%20An%20FDN%20with%20a%20*paraunitary*%20(FIR%2C%20lossless)%20scattering%20matrix%20in%20the%20loop.%0A%20%20%20%20The%20example%20computes%20the%20impulse%20response%20by%20time-domain%20recursion%20and%20by%0A%20%20%20%20modal%20decomposition%2C%20and%20verifies%20that%20the%20system%20is%20lossless%20(all%20poles%0A%20%20%20%20on%20the%20unit%20circle).%0A%0A%20%20%20%20Reference%3A%20*Schlecht%2C%20S.%2C%20Habets%2C%20E.%20(2020).%20Scattering%20in%20Feedback%20Delay%0A%20%20%20%20Networks.%20IEEE%2FACM%20Transactions%20on%20Audio%2C%20Speech%2C%20and%20Language%20Processing.*%0A%20%20%20%20%5Bdoi%3A10.1109%2Ftaslp.2020.3001395%5D(https%3A%2F%2Fdx.doi.org%2F10.1109%2Ftaslp.2020.3001395)%0A%0A%20%20%20%20Original%20MATLAB%3A%20%60example_paraunitaryFDN.m%60%2C%20Sebastian%20J.%20Schlecht%2C%0A%20%20%20%2023%20April%202018.%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%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20plotly.graph_objects%20as%20go%0A%20%20%20%20import%20plotly.io%20as%20pio%0A%20%20%20%20import%20torch%0A%0A%20%20%20%20import%20pyFDN%0A%0A%20%20%20%20pio.renderers.default%20%3D%20%22sphinx_gallery%22%0A%20%20%20%20return%20go%2C%20np%2C%20plt%2C%20pyFDN%2C%20torch%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%20Define%20FFDN%20with%20paraunitary%20feedback%20matrix%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20pyFDN)%3A%0A%20%20%20%20np.random.seed(2)%0A%20%20%20%20fs%20%3D%2048000%0A%20%20%20%20ir_len%20%3D%20int(fs%20%2F%2010)%0A%0A%20%20%20%20num_delays%20%3D%204%0A%20%20%20%20num_stages%20%3D%203%0A%20%20%20%20delays%20%3D%20np.random.randint(150%2C%20501%2C%20num_delays)%0A%20%20%20%20input_gain%20%3D%20np.eye(num_delays%2C%201)%0A%20%20%20%20output_gain%20%3D%20np.eye(1%2C%20num_delays)%0A%20%20%20%20direct%20%3D%20np.zeros((1%2C%201))%0A%0A%20%20%20%20feedback_matrix%20%3D%20pyFDN.filter_matrix_gallery(%0A%20%20%20%20%20%20%20%20num_delays%2C%20%22RandomDense%22%2C%20num_stages%3Dnum_stages%2C%20stage_matrix_type%3D%22random%22%0A%20%20%20%20)%0A%20%20%20%20print(f%22Delays%3A%20%7Bdelays%7D%20(sum%20%3D%20%7Bdelays.sum()%7D)%22)%0A%20%20%20%20print(f%22Feedback%20matrix%3A%20%7Bfeedback_matrix.shape%5B2%5D%7D%20taps%22)%0A%20%20%20%20return%20delays%2C%20direct%2C%20feedback_matrix%2C%20fs%2C%20input_gain%2C%20ir_len%2C%20output_gain%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%20Feedback%20matrix%20impulse%20responses%0A%0A%20%20%20%20Each%20subplot%20shows%20the%20FIR%20of%20one%20matrix%20entry.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(feedback_matrix%2C%20np%2C%20plt%2C%20pyFDN)%3A%0A%20%20%20%20pyFDN.plot_impulse_response_matrix(%0A%20%20%20%20%20%20%20%20np.arange(feedback_matrix.shape%5B2%5D)%2C%0A%20%20%20%20%20%20%20%20pyFDN.mulaw_encode(feedback_matrix.transpose(2%2C%200%2C%201))%2C%0A%20%20%20%20%20%20%20%20xlabel%3D%22Time%20(samples)%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3D%22Amplitude%20(mu)%22%2C%0A%20%20%20%20%20%20%20%20title%3D%22Paraunitary%20feedback%20matrix%22%2C%0A%20%20%20%20)%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%20Verify%20paraunitarity%0A%0A%20%20%20%20A%20paraunitary%20matrix%20satisfies%20%24A%5ET(z%5E%7B-1%7D)%5C%2CA(z)%20%3D%20I%24%3B%20in%20the%20time%20domain%0A%20%20%20%20the%20matrix%20impulse%20response%20is%20lossless.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(feedback_matrix%2C%20pyFDN)%3A%0A%20%20%20%20is_pu%2C%20_%2C%20max_off%20%3D%20pyFDN.is_paraunitary(feedback_matrix.transpose(2%2C%200%2C%201))%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22Feedback%20matrix%20is%20paraunitary%3A%20%7Bbool(is_pu)%7D%20(max%20off-diagonal%20%7Bmax_off%3A.2e%7D)%22%0A%20%20%20%20)%0A%20%20%20%20assert%20is_pu%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%20Impulse%20response%20and%20modal%20decomposition%0A%0A%20%20%20%20The%20FIR%20feedback%20matrix%20runs%20directly%20in%20the%20time-domain%20recursion.%0A%20%20%20%20For%20the%20modal%20decomposition%20the%20FIR%20matrix%20is%20placed%20as%20a%20FLAMO%20Filter%0A%20%20%20%20module%20in%20the%20loop%20(%60dss_to_flamo%60)%20and%20%60flamo_to_pr%60%20refines%20the%20poles%0A%20%20%20%20with%20Ehrlich%E2%80%93Aberth%20iteration.%20The%20FIR%20feedback%20adds%20poles%20beyond%20the%0A%20%20%20%20delay%20count%2C%20so%20the%20number%20of%20root%20seeds%20is%20set%20to%20the%20degree%20of%20the%0A%20%20%20%20generalized%20characteristic%20polynomial.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20delays%2C%0A%20%20%20%20direct%2C%0A%20%20%20%20feedback_matrix%2C%0A%20%20%20%20fs%2C%0A%20%20%20%20input_gain%2C%0A%20%20%20%20ir_len%2C%0A%20%20%20%20np%2C%0A%20%20%20%20output_gain%2C%0A%20%20%20%20pyFDN%2C%0A%20%20%20%20torch%2C%0A)%3A%0A%20%20%20%20ir_time%20%3D%20pyFDN.dss_to_impz(%0A%20%20%20%20%20%20%20%20ir_len%2C%20delays%2C%20feedback_matrix%2C%20input_gain%2C%20output_gain%2C%20direct%0A%20%20%20%20)%5B%3A%2C%200%2C%200%5D%0A%0A%20%20%20%20%23%20pole%20count%20%3D%20degree%20of%20the%20generalized%20characteristic%20polynomial%20in%20w%20%3D%20z%5E%7B-1%7D%0A%20%20%20%20gcp%20%3D%20pyFDN.general_char_poly(delays%2C%20feedback_matrix)%0A%20%20%20%20num_poles%20%3D%20int(np.nonzero(np.abs(gcp)%20%3E%201e-12%20*%20np.max(np.abs(gcp)))%5B0%5D%5B-1%5D)%0A%0A%20%20%20%20model%20%3D%20pyFDN.dss_to_flamo(%0A%20%20%20%20%20%20%20%20A%3Dfeedback_matrix%2C%0A%20%20%20%20%20%20%20%20B%3Dinput_gain%2C%0A%20%20%20%20%20%20%20%20C%3Doutput_gain%2C%0A%20%20%20%20%20%20%20%20D%3Ddirect%2C%0A%20%20%20%20%20%20%20%20m%3Ddelays%2C%0A%20%20%20%20%20%20%20%20Fs%3Dfs%2C%0A%20%20%20%20%20%20%20%20shell%3DFalse%2C%0A%20%20%20%20%20%20%20%20dtype%3Dtorch.float64%2C%0A%20%20%20%20)%0A%20%20%20%20residues%2C%20poles%2C%20direct_term%2C%20is_pair%2C%20_meta%20%3D%20pyFDN.flamo_to_pr(%0A%20%20%20%20%20%20%20%20model%2C%0A%20%20%20%20%20%20%20%20quality_threshold%3D1e-10%2C%0A%20%20%20%20%20%20%20%20refinement_tol%3D1e-10%2C%0A%20%20%20%20%20%20%20%20maximum_iterations%3D80%2C%0A%20%20%20%20%20%20%20%20deflation_type%3D%22fullDeflation%22%2C%0A%20%20%20%20%20%20%20%20verbose%3DFalse%2C%0A%20%20%20%20%20%20%20%20num_poles%3Dnum_poles%2C%0A%20%20%20%20)%0A%20%20%20%20ir_modal%20%3D%20pyFDN.pr_to_impz(%0A%20%20%20%20%20%20%20%20residues%2C%20poles%2C%20direct_term%2C%20is_pair%2C%20ir_len%2C%20mode%3D%22lowMemory%22%0A%20%20%20%20)%5B%3A%2C%200%2C%200%5D%0A%0A%20%20%20%20difference%20%3D%20ir_time%20-%20ir_modal%0A%20%20%20%20max_deviation%20%3D%20np.max(np.abs(difference))%0A%20%20%20%20print(f%22Number%20of%20poles%3A%20%7Bpoles.size%7D%22)%0A%20%20%20%20print(f%22Max%20%7CIR_time%20-%20IR_modal%7C%20%3D%20%7Bmax_deviation%3A.3e%7D%22)%0A%20%20%20%20assert%20pyFDN.is_almost_zero(difference%2C%20tol%3D1e-3)%0A%20%20%20%20return%20difference%2C%20ir_modal%2C%20ir_time%2C%20poles%2C%20residues%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%20Time-domain%20vs%20modal%20reconstruction%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(difference%2C%20ir_modal%2C%20ir_time%2C%20pyFDN)%3A%0A%20%20%20%20pyFDN.plot_impulse_response(%0A%20%20%20%20%20%20%20%20difference%2C%0A%20%20%20%20%20%20%20%20ir_time%2C%0A%20%20%20%20%20%20%20%20ir_modal%2C%0A%20%20%20%20%20%20%20%20labels%3D%5B%22Difference%22%2C%20%22Time%20domain%22%2C%20%22Poles%2Fresidues%22%5D%2C%0A%20%20%20%20%20%20%20%20title%3D%22Impulse%20response%3A%20time-domain%20recursion%20vs%20modal%20synthesis%22%2C%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Poles%20and%20residues%0A%0A%20%20%20%20The%20FFDN%20is%20lossless%3A%20all%20pole%20magnitudes%20are%200%20dB.%20The%20residue%20magnitudes%0A%20%20%20%20spread%20over%20a%20wide%20range.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(go%2C%20np%2C%20poles%2C%20pyFDN%2C%20residues)%3A%0A%20%20%20%20fig_pr%20%3D%20go.Figure()%0A%20%20%20%20fig_pr.add_trace(%0A%20%20%20%20%20%20%20%20go.Scatter(%0A%20%20%20%20%20%20%20%20%20%20%20%20x%3Dnp.angle(poles)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20y%3DpyFDN.lin_to_db(np.abs(poles))%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mode%3D%22markers%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20marker%3D%7B%22size%22%3A%204%7D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20name%3D%22Poles%22%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20)%0A%20%20%20%20fig_pr.add_trace(%0A%20%20%20%20%20%20%20%20go.Scatter(%0A%20%20%20%20%20%20%20%20%20%20%20%20x%3Dnp.angle(poles)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20y%3DpyFDN.lin_to_db(np.abs(residues%5B%3A%2C%200%2C%200%5D))%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mode%3D%22markers%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20marker%3D%7B%22size%22%3A%204%2C%20%22symbol%22%3A%20%22x%22%7D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20name%3D%22Residues%22%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20)%0A%20%20%20%20fig_pr.update_layout(%0A%20%20%20%20%20%20%20%20title%3D%22Pole%20and%20residue%20magnitudes%20over%20pole%20angle%22%2C%0A%20%20%20%20%20%20%20%20xaxis%3D%7B%22title%22%3A%20%22Pole%20angle%20(rad)%22%7D%2C%0A%20%20%20%20%20%20%20%20yaxis%3D%7B%22title%22%3A%20%22Magnitude%20(dB)%22%7D%2C%0A%20%20%20%20%20%20%20%20template%3D%22plotly_white%22%2C%0A%20%20%20%20%20%20%20%20height%3D420%2C%0A%20%20%20%20)%0A%20%20%20%20fig_pr.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%20Losslessness%20check%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20poles%2C%20pyFDN)%3A%0A%20%20%20%20max_pole_deviation%20%3D%20np.max(np.abs(np.abs(poles)%20-%201.0))%0A%20%20%20%20print(f%22Max%20%7C%20%7Cpole%7C%20-%201%20%7C%20%3D%20%7Bmax_pole_deviation%3A.3e%7D%22)%0A%20%20%20%20assert%20pyFDN.is_almost_zero(np.abs(poles)%20-%201.0)%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
04d1921110cbcad3e5cca2e403a9cdc0