filters chapter

import%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%20os%0A%20%20%20%20import%20time%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20from%20scipy%20import%20signal%0A%20%20%20%20from%20scipy.signal%20import%20firwin2%2C%20convolve%2C%20fftconvolve%2C%20lfilter%0A%0A%20%20%20%20return%20convolve%2C%20fftconvolve%2C%20firwin2%2C%20lfilter%2C%20np%2C%20plt%2C%20signal%2C%20time%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20signal)%3A%0A%20%20%20%20num_taps%20%3D%2051%0A%20%20%20%20cut_off%20%3D%203000%20%20%23%20Hz%0A%20%20%20%20sample_rate%20%3D%2032000%20%20%23%20Hz%0A%0A%20%20%20%20%23%20create%20our%20low%20pass%20filter%0A%20%20%20%20h%20%3D%20signal.firwin(num_taps%2C%20cut_off%2C%20fs%3Dsample_rate)%0A%0A%20%20%20%20%23%20plot%20the%20impulse%20response%0A%20%20%20%20plt.plot(h%2C%20'.-')%0A%20%20%20%20plt.show()%0A%0A%20%20%20%20%23%20plot%20the%20frequency%20response%0A%20%20%20%20H%20%3D%20np.abs(np.fft.fft(h%2C%201024))%20%20%23%20take%20the%201024-point%20FFT%20and%20magnitude%0A%20%20%20%20H%20%3D%20np.fft.fftshift(H)%20%20%23%20make%200%20Hz%20in%20the%20center%0A%20%20%20%20w%20%3D%20np.linspace(-sample_rate%20%2F%202%2C%20sample_rate%20%2F%202%2C%20len(H))%20%20%23%20x%20axis%0A%20%20%20%20plt.plot(w%2C%20H%2C%20'.-')%0A%20%20%20%20plt.show()%0A%0A%20%20%20%20%23%20(h%20was%20found%20using%20the%20first%20code%20snippet)%0A%0A%20%20%20%20%23%20Shift%20the%20filter%20in%20frequency%20by%20multiplying%20by%20exp(j*2*pi*f0*t)%0A%20%20%20%20f0%20%3D%2010e3%20%20%23%20amount%20we%20will%20shift%0A%20%20%20%20Ts%20%3D%201.0%20%2F%20sample_rate%20%20%23%20sample%20period%0A%20%20%20%20t%20%3D%20np.arange(0.0%2C%20Ts%20*%20len(h)%2C%20Ts)%20%20%23%20time%20vector.%20args%20are%20(start%2C%20stop%2C%20step)%0A%20%20%20%20exponential%20%3D%20np.exp(2j%20*%20np.pi%20*%20f0%20*%20t)%20%20%23%20this%20is%20essentially%20a%20complex%20sine%20wave%0A%0A%20%20%20%20h_band_pass%20%3D%20h%20*%20exponential%20%20%23%20do%20the%20shift%0A%0A%20%20%20%20%23%20plot%20impulse%20response%0A%20%20%20%20plt.subplot(121)%0A%20%20%20%20plt.plot(np.real(h_band_pass)%2C%20'.-')%0A%20%20%20%20plt.plot(np.imag(h_band_pass)%2C%20'.-')%0A%20%20%20%20plt.legend(%5B'real'%2C%20'imag'%5D%2C%20loc%3D1)%0A%0A%20%20%20%20%23%20plot%20the%20frequency%20response%0A%20%20%20%20H%20%3D%20np.abs(np.fft.fft(h_band_pass%2C%201024))%20%20%23%20take%20the%201024-point%20FFT%20and%20magnitude%0A%20%20%20%20H%20%3D%20np.fft.fftshift(H)%20%20%23%20make%200%20Hz%20in%20the%20center%0A%20%20%20%20w%20%3D%20np.linspace(-sample_rate%20%2F%202%2C%20sample_rate%20%2F%202%2C%20len(H))%20%20%23%20x%20axis%0A%20%20%20%20plt.subplot(122)%0A%20%20%20%20plt.plot(w%2C%20H%2C%20'.-')%0A%20%20%20%20plt.xlabel('Frequency%20%5BHz%5D')%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(convolve%2C%20fftconvolve%2C%20firwin2%2C%20lfilter%2C%20np)%3A%0A%20%20%20%20%23%20comparing%20different%20convolution%20methods%0A%20%20%20%20%23%20Create%20a%20test%20signal%2C%20we'll%20use%20Gaussian%20noise%0A%20%20%20%20_sample_rate%20%3D%201e6%20%23%20Hz%0A%20%20%20%20_N%20%3D%201000%20%23%20samples%20to%20simulate%0A%20%20%20%20_x%20%3D%20np.random.randn(_N)%20%2B%201j%20*%20np.random.randn(_N)%0A%0A%20%20%20%20%23%20Create%20an%20FIR%20filter%2C%20same%20one%20as%202nd%20example%20above%0A%20%20%20%20freqs%20%3D%20%5B0%2C%20100e3%2C%20110e3%2C%20190e3%2C%20200e3%2C%20300e3%2C%20310e3%2C%20500e3%5D%0A%20%20%20%20gains%20%3D%20%5B1%2C%201%2C%20%20%20%20%200%2C%20%20%20%20%200%2C%20%20%20%20%200.5%2C%20%20%200.5%2C%20%20%200%2C%20%20%20%20%200%5D%0A%20%20%20%20h2%20%3D%20firwin2(101%2C%20freqs%2C%20gains%2C%20fs%3D_sample_rate)%0A%0A%20%20%20%20%23%20Apply%20filter%20using%20the%20four%20different%20methods%0A%20%20%20%20x_numpy%20%3D%20np.convolve(h2%2C%20_x)%0A%20%20%20%20x_scipy%20%3D%20convolve(h2%2C%20_x)%20%23%20scipys%20convolve%0A%20%20%20%20x_fft_convolve%20%3D%20fftconvolve(h2%2C%20_x)%0A%20%20%20%20x_lfilter%20%3D%20lfilter(h2%2C%201%2C%20_x)%20%23%202nd%20arg%20is%20always%201%20for%20FIR%20filters%0A%0A%20%20%20%20%23%20Prove%20they%20are%20all%20giving%20the%20same%20output%0A%20%20%20%20print(x_numpy%5B0%3A2%5D)%0A%20%20%20%20print(x_scipy%5B0%3A2%5D)%0A%20%20%20%20print(x_fft_convolve%5B0%3A2%5D)%0A%20%20%20%20print(x_lfilter%5B0%3A2%5D)%0A%20%20%20%20return%20(h2%2C)%0A%0A%0A%40app.cell%0Adef%20_(convolve%2C%20fftconvolve%2C%20lfilter%2C%20np%2C%20plt%2C%20time)%3A%0A%20%20%20%20%23%20Make%20plot%20comparing%20methods%20using%20different%20taps%2C%20this%20takes%20a%20while%20to%20run%20(even%20longer%20in%20the%20web-based%20marimo%20notebook)%20so%20it's%20off%20by%20default%0A%20%20%20%20if%20False%3A%0A%20%20%20%20%20%20%20%20%23_taps_sizes%20%3D%20np.logspace(1%2C%203%2C%20num%3D100)%20%23%2010**x%20is%20the%20start%2Fstop%0A%20%20%20%20%20%20%20%20_taps_sizes%20%3D%20np.linspace(10%2C%2010000%2C%20num%3D100)%0A%20%20%20%20%20%20%20%20_N%20%3D%20100000%20%23%20signal%20length%0A%20%20%20%20%20%20%20%20%23_N%20%3D%201000%20%23%20signal%20length%0A%20%20%20%20%20%20%20%20_x%20%3D%20np.random.randn(_N)%20%2B%201j%20*%20np.random.randn(_N)%20%23%20complex%20signal%0A%20%20%20%20%20%20%20%20_results%20%3D%20np.zeros((4%2C%20len(_taps_sizes)))%0A%20%20%20%20%20%20%20%20for%20_method_i%20in%20range(4)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20_i%20in%20range(len(_taps_sizes))%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_h%20%3D%20np.random.randn(int(_taps_sizes%5B_i%5D))%20%23%20real%20taps%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_num_trials%20%3D%20int(10000%20%2F%20_taps_sizes%5B_i%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23_num_trials%20%3D%20int(1000000%20%2F%20_taps_sizes%5B_i%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_start_t%20%3D%20time.time()%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20_method_i%20%3D%3D%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20_ii%20in%20range(_num_trials)%3A%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%3D%20np.convolve(_h%2C%20_x)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20elif%20_method_i%20%3D%3D%201%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20_ii%20in%20range(_num_trials)%3A%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%3D%20convolve(_h%2C%20_x)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20elif%20_method_i%20%3D%3D%202%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20_ii%20in%20range(_num_trials)%3A%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%3D%20fftconvolve(_h%2C%20_x)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20elif%20_method_i%20%3D%3D%203%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20_ii%20in%20range(_num_trials)%3A%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%3D%20lfilter(_h%2C%201%2C%20_x)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_results%5B_method_i%2C%20_i%5D%20%3D%20(time.time()%20-%20_start_t)%2F_num_trials*1e3%20%23%20in%20ms%0A%0A%20%20%20%20%20%20%20%20%23%20Lots%20of%20folks%20use%20semilogx%2C%20but%20its%20nice%20to%20see%20the%20linear%20aspect%0A%20%20%20%20%20%20%20%20plt.plot(_taps_sizes%2C%20_results%5B0%2C%3A%5D%2C%20'.-')%0A%20%20%20%20%20%20%20%20plt.plot(_taps_sizes%2C%20_results%5B1%2C%3A%5D%2C%20'.-')%0A%20%20%20%20%20%20%20%20plt.plot(_taps_sizes%2C%20_results%5B2%2C%3A%5D%2C%20'.-')%0A%20%20%20%20%20%20%20%20plt.plot(_taps_sizes%2C%20_results%5B3%2C%3A%5D%2C%20'.-')%0A%20%20%20%20%20%20%20%20plt.title('Input%20Signal%20Length%3A%20'%20%2B%20str(_N)%20%2B%20'%20samples')%0A%20%20%20%20%20%20%20%20plt.xlabel('Number%20of%20taps')%0A%20%20%20%20%20%20%20%20plt.ylabel('Time%20per%20call%20(ms)')%0A%20%20%20%20%20%20%20%20plt.legend(%5B'np.convolve'%2C%20'scipy.signal.convolve'%2C%20'scipy.signal.fftconvolve'%2C%20'scipy.signal.lfilter'%5D)%0A%20%20%20%20%20%20%20%20%23if%20os.path.isdir('..%2F_images')%3A%0A%20%20%20%20%20%20%20%20%23%20%20%20%20plt.savefig('..%2F_images%2Fconvolve_comparison_'%2Bstr(_N)%2B'.svg'%2C%20bbox_inches%3D'tight')%0A%20%20%20%20%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(fftconvolve%2C%20h2%2C%20np%2C%20plt)%3A%0A%20%20%20%20_sample_rate%20%3D%201e6%0A%0A%20%20%20%20%23%20(h2%20was%20created%20above)%0A%0A%20%20%20%20%23%20Simulate%20signal%20comprising%20of%20Gaussian%20noise%0A%20%20%20%20_N%20%3D%20100000%20%23%20signal%20length%0A%20%20%20%20_x%20%3D%20np.random.randn(_N)%20%2B%201j%20*%20np.random.randn(_N)%20%23%20complex%20signal%0A%0A%20%20%20%20%23%20Save%20PSD%20of%20the%20input%20signal%0A%20%20%20%20PSD_input%20%3D%2010*np.log10(np.fft.fftshift(np.abs(np.fft.fft(_x))**2)%2Flen(_x))%0A%0A%20%20%20%20%23%20Apply%20filter%0A%20%20%20%20_x%20%3D%20fftconvolve(_x%2C%20h2%2C%20'same')%0A%0A%20%20%20%20%23%20Look%20at%20PSD%20of%20the%20output%20signal%0A%20%20%20%20PSD_output%20%3D%2010*np.log10(np.fft.fftshift(np.abs(np.fft.fft(_x))**2)%2Flen(_x))%0A%20%20%20%20f%20%3D%20np.linspace(-_sample_rate%2F2%2F1e6%2C%20_sample_rate%2F2%2F1e6%2C%20len(PSD_output))%0A%20%20%20%20plt.plot(f%2C%20PSD_input%2C%20alpha%3D0.8)%0A%20%20%20%20plt.plot(f%2C%20PSD_output%2C%20alpha%3D0.8)%0A%20%20%20%20plt.xlabel('Frequency%20%5BMHz%5D')%0A%20%20%20%20plt.ylabel('PSD%20%5BdB%5D')%0A%20%20%20%20plt.axis(%5B_sample_rate%2F-2%2F1e6%2C%20_sample_rate%2F2%2F1e6%2C%20-40%2C%2020%5D)%0A%20%20%20%20plt.legend(%5B'Input'%2C%20'Output'%5D%2C%20loc%3D1)%0A%20%20%20%20plt.grid()%0A%20%20%20%20%23if%20os.path.isdir('..%2F_images')%3A%0A%20%20%20%20%23%20%20%20%20plt.savefig('..%2F_images%2Ffftconvolve.svg'%2C%20bbox_inches%3D'tight')%0A%20%20%20%20plt.show()%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