Signal Processing Basics

Open in molab View on GitHub Download .py

mo.md(%22%5Cn%20%20%20%20%23%20Signal%20Processing%20Basics%5Cn%20%20%20%20*Written%20by%20Luke%20Chang*%5Cn%5Cn%20%20%20%20In%20this%20lab%2C%20we%20will%20cover%20the%20basics%20of%20convolution%2C%20sine%20waves%2C%20and%20fourier%20transforms.%20This%20lab%20is%20largely%20based%20on%20exercises%20from%20Mike%20X%20Cohen's%20excellent%20book%2C%20%5BAnalyzing%20Neural%20Data%20Analysis%3A%20Theory%20and%20Practice%5D(https%3A%2F%2Fwww.amazon.com%2FAnalyzing-Neural-Time-Data-Practice%2Fdp%2F0262019876).%20If%20you%20are%20interested%20in%20learning%20in%20more%20detail%20about%20the%20basics%20of%20EEG%20and%20time-series%20analyses%20I%20highly%20recommend%20his%20accessible%20introduction.%20I%20also%20encourage%20you%20to%20watch%20his%20accompanying%20freely%20available%20%5B*lecturelets*%5D(https%3A%2F%2Fwww.youtube.com%2Fchannel%2FUCUR_LsXk7IYyueSnXcNextQ)%20to%20learn%20more%20about%20each%20topic%20introduced%20in%20this%20notebook.%5Cn%5Cn%20%20%20%20%23%23%20Time%20Domain%5Cn%5Cn%20%20%20%20First%20we%20will%20work%20on%20signals%20in%20the%20time%20domain.%20This%20requires%20measuring%20a%20signal%20at%20a%20constant%20interval%20over%20time.%20The%20frequency%20with%20which%20we%20measure%20a%20signal%20is%20referred%20to%20as%20the%20sampling%20frequency.%20The%20units%20of%20this%20are%20typically%20described%20in%20%24Hz%24%20-%20or%20the%20number%20of%20cycles%20per%20second.%20It%20is%20critical%20that%20the%20sampling%20frequency%20is%20consistent%20over%20the%20entire%20measurement%20of%20the%20time%20series.%5Cn%5Cn%20%20%20%20%23%23%23%20Dot%20Product%5Cn%20%20%20%20To%20understand%20convolution%2C%20we%20first%20need%20to%20familiarize%20ourselves%20with%20the%20dot%20product.%20The%20dot%20product%20is%20simply%20the%20sum%20of%20the%20elements%20of%20a%20vector%20weighted%20by%20the%20elements%20of%20another%20vector.%20This%20method%20is%20commonly%20used%20in%20signal%20processing%2C%20and%20also%20in%20statistics%20as%20a%20measure%20of%20similarity%20between%20two%20vectors.%20Finally%2C%20there%20is%20also%20a%20geometric%20interpretation%20which%20is%20a%20mapping%20between%20vectors%20(i.e.%2C%20the%20product%20of%20the%20magnitudes%20of%20the%20two%20vectors%20scaled%20by%20the%20cosine%20of%20the%20angle%20between%20them).%20For%20a%20more%20in%20depth%20overview%20of%20the%20dot%20product%20and%20its%20relation%20to%20convolution%2C%20you%20can%20watch%20this%20optional%20%5Bvideo%5D(https%3A%2F%2Fyoutu.be%2Frea6M1oagmA).%5Cn%5Cn%20%20%20%20%24%24dotproduct_%7Bab%7D%3D%5C%5Csum%5C%5Climits_%7Bi%3D1%7D%5En%20a_i%20b_i%24%24%5Cn%5Cn%20%20%20%20Let's%20create%20some%20vectors%20of%20random%20numbers%20and%20see%20how%20the%20dot%20product%20works.%20First%2C%20the%20two%20vectors%20need%20to%20be%20of%20the%20same%20length.%5Cn%20%20%20%20%22)

Dot Product: 644

Dot Product: 666.66

You can see that after convolution, each spike has now become the shape of the kernel. Signal: 100, Kernel: 10, Convolved: 109 samples

_t_true%20%3D%20arange(-1%2C%201%2C%201%20%2F%201000)%0A_t%20%3D%20arange(-1%2C%201%20%2B%201%20%2F%20sf_slider.value%2C%201%20%2F%20sf_slider.value)%0A_waves_true%20%3D%20np.array(%5Bmulti_amps%5B_j%5D%20*%20sin(2%20*%20pi%20*%20_f%20*%20_t_true%20%2B%20multi_phases%5B_j%5D)%20for%20_j%2C%20_f%20in%20enumerate(multi_freqs)%5D)%0A_waves%20%3D%20np.array(%5Bmulti_amps%5B_j%5D%20*%20sin(2%20*%20pi%20*%20_f%20*%20_t%20%2B%20multi_phases%5B_j%5D)%20for%20_j%2C%20_f%20in%20enumerate(multi_freqs)%5D)%0Acombined_signal%20%3D%20np.sum(_waves%2C%20axis%3D0)%20%2B%20noise_slider.value%20*%20np.random.randn(len(_t))%0A_combined_true%20%3D%20np.sum(_waves_true%2C%20axis%3D0)%0A_nyquist%20%3D%20sf_slider.value%20%2F%202%0A%0Adef%20_apparent(f%2C%20sf)%3A%0A%20%20%20%20folded%20%3D%20f%20%25%20sf%0A%20%20%20%20return%20folded%20if%20folded%20%3C%3D%20sf%20%2F%202%20else%20sf%20-%20folded%0A_fig%2C%20_axes%20%3D%20plt.subplots(nrows%3D6%2C%20figsize%3D(12%2C%2011)%2C%20sharex%3DTrue)%0Afor%20_j%2C%20_f%20in%20enumerate(multi_freqs)%3A%0A%20%20%20%20_ax%20%3D%20_axes%5B_j%5D%0A%20%20%20%20_ax.plot(_t_true%2C%20_waves_true%5B_j%5D%2C%20color%3D'gray'%2C%20alpha%3D0.45%2C%20linewidth%3D1%2C%20label%3D'true%20signal')%0A%20%20%20%20_is_aliased%20%3D%20_f%20%3E%20_nyquist%0A%20%20%20%20_color%20%3D%20'tab%3Ared'%20if%20_is_aliased%20else%20'tab%3Ablue'%0A%20%20%20%20_ax.plot(_t%2C%20_waves%5B_j%5D%2C%20'o-'%2C%20color%3D_color%2C%20linewidth%3D1.3%2C%20markersize%3D4%2C%20label%3D'sampled')%0A%20%20%20%20_ylabel%20%3D%20f'%7B_f%7D%20Hz'%0A%20%20%20%20if%20_is_aliased%3A%0A%20%20%20%20%20%20%20%20_ylabel%20%2B%3D%20f'%5Cn%E2%86%92%20looks%20like%5Cn%7B_apparent(_f%2C%20sf_slider.value)%3A.1f%7D%20Hz'%0A%20%20%20%20_ax.set_ylabel(_ylabel%2C%20fontsize%3D10)%0A%20%20%20%20_ax.grid(alpha%3D0.3)%0A_axes%5B0%5D.set_title(f'Sampling%3A%20%7Bsf_slider.value%7D%20Hz%20%20%20%7C%20%20%20Nyquist%3A%20%7B_nyquist%7D%20Hz'%2C%20fontsize%3D16)%0A_axes%5B0%5D.legend(loc%3D'upper%20right'%2C%20fontsize%3D8%2C%20ncol%3D2)%0A_axes%5B5%5D.plot(_t_true%2C%20_combined_true%2C%20color%3D'gray'%2C%20alpha%3D0.45%2C%20linewidth%3D1%2C%20label%3D'true%20sum')%0A_axes%5B5%5D.plot(_t%2C%20combined_signal%2C%20'o-'%2C%20color%3D'purple'%2C%20linewidth%3D1.3%2C%20markersize%3D4%2C%20label%3Df'sampled%20%2B%20noise%20(%7Bnoise_slider.value%7D)')%0A_axes%5B5%5D.set_ylabel('Sum'%2C%20fontsize%3D10)%0A_axes%5B5%5D.set_xlabel('Time%20(s)'%2C%20fontsize%3D14)%0A_axes%5B5%5D.legend(loc%3D'upper%20right'%2C%20fontsize%3D8%2C%20ncol%3D2)%0A_axes%5B5%5D.grid(alpha%3D0.3)%0A_axes%5B0%5D.set_xlim(-0.25%2C%200.25)%0Aplt.tight_layout()%0Aplt.close()%0A_aliased%20%3D%20%5B(f%2C%20_apparent(f%2C%20sf_slider.value))%20for%20f%20in%20multi_freqs%20if%20f%20%3E%20_nyquist%5D%0Aif%20_aliased%3A%0A%20%20%20%20_msg%20%3D%20'**Aliasing!**%20'%20%2B%20'%2C%20'.join((f'%7Bf%7D%20Hz%20%E2%86%92%20appears%20as%20%7Ba%3A.1f%7D%20Hz'%20for%20f%2C%20a%20in%20_aliased))%20%2B%20f'%20%20(Nyquist%20%3D%20%7B_nyquist%7D%20Hz)'%0A%20%20%20%20_warn%20%3D%20mo.callout(mo.md(_msg)%2C%20kind%3D'warn')%0Aelse%3A%0A%20%20%20%20_warn%20%3D%20mo.callout(mo.md(f'All%20frequencies%20below%20Nyquist%20(%7B_nyquist%7D%20Hz)%20%E2%80%94%20no%20aliasing.')%2C%20kind%3D'success')%0Amo.vstack(%5Bmo.hstack(%5Bsf_slider%2C%20noise_slider%5D)%2C%20_warn%2C%20_fig%5D)

mo.md(%22%5Cn%20%20%20%20%23%23%20Frequency%20Domain%5Cn%5Cn%20%20%20%20In%20the%20previous%20example%2C%20we%20generated%20a%20complex%20signal%20composed%20of%20multiple%20sine%20waves%20oscillating%20at%20different%20frequencies.%20Typically%20in%20data%20analysis%2C%20we%20only%20observe%20the%20signal%20and%20are%20trying%20to%20uncover%20the%20generative%20processes%20that%20gave%20rise%20to%20the%20signal.%20In%20this%20section%2C%20we%20will%20introduce%20the%20frequency%20domain%20and%20how%20we%20can%20identify%20if%20there%20are%20any%20frequencies%20oscillating%20at%20a%20consistent%20frequency%20in%20our%20signal%20using%20the%20fourier%20transform.%20The%20fourier%20transform%20convolves%20different%20frequencies%20of%20sine%20waves%20with%20our%20data%20to%20identify%20oscillatory%20components.%5Cn%5Cn%20%20%20%20One%20important%20assumption%3A%20**stationarity**%20%E2%80%94%20the%20generative%20processes%20don't%20vary%20over%20time.%5Cn%5Cn%20%20%20%20See%20this%20%5Bvideo%5D(https%3A%2F%2Fyoutu.be%2Frea6M1oagmA)%20or%20a%20more%20in%20depth%20discussion%20on%20stationarity.%20In%20practice%2C%20this%20assumption%20is%20rarely%20true.%20Often%20it%20can%20be%20useful%20to%20use%20other%20techniques%20such%20as%20wavelets%20to%20look%20at%20time%20x%20frequency%20representations.%20We%20will%20not%20be%20covering%20wavelets%20here%2C%20but%20see%20this%20series%20of%20%5Bvideos%5D(https%3A%2F%2Fwww.youtube.com%2Fwatch%3Fv%3D7ahrcB5HL0k)%20for%20more%20information.%5Cn%5Cn%20%20%20%20%23%23%23%20Discrete%20Time%20Fourier%20Transform%5Cn%20%20%20%20We%20will%20gain%20an%20intution%20of%20how%20the%20fourier%20transform%20works%20by%20building%20our%20own%20discrete%20time%20fourier%20transform.%5Cn%5Cn%20%20%20%20Let%E2%80%99s%20watch%20this%20short%20video%20about%20the%20fourier%20transform%20by%20Mike%20X%20Cohen.%20Don%E2%80%99t%20worry%20too%20much%20about%20the%20details%20of%20the%20discussion%20on%20the%20matlab%20code%20as%20we%20will%20be%20exploring%20these%20concepts%20in%20python%20below.%5Cn%20%20%20%20%22)

mo.md('%5Cn%20%20%20%20The%20discrete%20Fourier%20transform%20of%20variable%20%24x%24%20at%20frequency%5Cn%20%20%20%20%24f%24%20can%20be%20defined%20as%3A%5Cn%5Cn%20%20%20%20%24%24X_f%20%3D%20%5C%5Csum%5C%5Climits_%7Bk%3D0%7D%5E%7Bn-1%7D%20x_k%20%5C%5Ccdot%20e%5E%7B%5C%5Cfrac%7B-i2%5C%5Cpi%20fk%7D%7Bn%7D%7D%24%24%5Cn%5Cn%20%20%20%20where%20%24n%24%20refers%20to%20the%20number%20of%20data%20points%20in%20vector%20%24x%24%2C%20and%20the%20capital%20letter%20%24X_f%24%20is%20the%20fourier%20coefficient%20of%20time%20series%20variable%20%24x%24%20at%20frequency%20%24f%24.%5Cn%5Cn%20%20%20%20Essentially%2C%20we%20create%20a%20bank%20of%20complex%20sine%20waves%20at%20different%20frequencies%20that%20are%20linearly%20spaced.%20The%20zero%20frequency%20component%20reflects%20the%20mean%20offset%20over%20the%20entire%20signal%20and%20will%20simply%20be%20zero%20in%20our%20example.%5Cn%5Cn%20%20%20%20%23%23%23%23%20Complex%20Sine%20Waves%5Cn%20%20%20%20You%20may%20have%20noticed%20that%20we%20are%20computing%20complex%20sine%20waves%20using%20the%20%60np.exp%60%20function%20instead%20of%20the%20%60np.sin%60%20function.%5Cn%5Cn%20%20%20%20%24%24e%5E%7Bi(2%5C%5Cpi%20ft%20%2B%20%5C%5Ctheta)%7D%24%24%5Cn%5Cn%20%20%20%20We%20will%20not%20spend%20too%20much%20time%20on%20the%20details%2C%20but%20basically%20complex%20sine%20waves%20have%20three%20components%3A%20time%2C%20a%20real%20part%20of%20the%20sine%20wave%2C%20and%20the%20imaginary%20part%20of%20the%20sine%20wave%2C%20which%20are%20basically%20phase%20shifted%20by%20%24%5C%5Cfrac%7B%5C%5Cpi%7D%7B2%7D%24.%20%241j%24%20is%20how%20we%20can%20specify%20a%20complex%20number%20in%20python.%20We%20can%20extract%20the%20real%20components%20using%20%60np.real%60%20or%20the%20imaginary%20using%20%60np.imag%60.%5Cn%5Cn%20%20%20%20We%20can%20visualize%20complex%20sine%20waves%20in%20three%20dimensions.%20For%20more%20information%2C%20watch%20this%20%5Bvideo%5D(https%3A%2F%2Fwww.youtube.com%2Fwatch%3Fv%3DiZCDOuzfsY0).%20If%20you%20need%20a%20refresher%20on%20complex%20numbers%2C%20you%20may%20want%20to%20watch%20this%20%5Bvideo%5D(https%3A%2F%2Fwww.youtube.com%2Fwatch%3Fv%3DfNfXKiIIufY).%5Cn%5Cn%20%20%20%20In%20this%20plot%2C%20we%20show%20this%20complex%20signal%20in%203%20dimensions%20and%20also%20project%20on%20two%20dimensional%20planes%20to%20show%20that%20the%20real%20and%20imaginary%20create%20a%20unit%20circle%2C%20and%20are%20phase%20offset%20by%20%24%5C%5Cfrac%7B%5C%5Cpi%7D%7B2%7D%24%20with%20respect%20to%20time.%5Cn%20%20%20%20')

from%20matplotlib.collections%20import%20LineCollection%0Aimport%20matplotlib.patches%20as%20mpatches%0A_time%20%3D%20arange(-1%2C%201%20%2B%201%20%2F%20osc_sf%2C%201%20%2F%20osc_sf)%0A_z%20%3D%20exp(1j%20*%20(2%20*%20pi%20*%20cx_freq_slider.value%20*%20_time%20%2B%20cx_theta_slider.value))%0A_phi%20%3D%20cx_theta_slider.value%0A_t0%20%3D%20int(np.argmin(np.abs(_time)))%0A_x0%2C%20_y0%20%3D%20(real(_z)%5B_t0%5D%2C%20imag(_z)%5B_t0%5D)%0A_fig%20%3D%20plt.figure(figsize%3D(15%2C%2010))%0A_ax1%20%3D%20_fig.add_subplot(2%2C%202%2C%201%2C%20projection%3D'3d')%0A_ax1.plot(np.arange(len(_time))%20%2F%20osc_sf%2C%20real(_z)%2C%20imag(_z))%0A_ax1.set_xlabel('Time%20(sec)'%2C%20fontsize%3D14)%0A_ax1.set_ylabel('Real(z)'%2C%20fontsize%3D14)%0A_ax1.set_zlabel('Imaginary(z)'%2C%20fontsize%3D14)%0A_ax1.set_title('Complex%20Sine%20Wave'%2C%20fontsize%3D16)%0A_ax1.view_init(15%2C%20250)%0A_ax2%20%3D%20_fig.add_subplot(2%2C%202%2C%202)%0A_pts%20%3D%20np.column_stack(%5Breal(_z)%2C%20imag(_z)%5D)%0A_segs%20%3D%20np.concatenate(%5B_pts%5B%3A-1%2C%20None%2C%20%3A%5D%2C%20_pts%5B1%3A%2C%20None%2C%20%3A%5D%5D%2C%20axis%3D1)%0A_lc%20%3D%20LineCollection(_segs%2C%20cmap%3D'viridis'%2C%20array%3D_time%5B%3A-1%5D%2C%20linewidth%3D1.8)%0A_ax2.add_collection(_lc)%0A_u%20%3D%20np.linspace(0%2C%202%20*%20pi%2C%20200)%0A_ax2.plot(np.cos(_u)%2C%20np.sin(_u)%2C%20color%3D'gray'%2C%20alpha%3D0.25%2C%20linestyle%3D'--'%2C%20linewidth%3D1)%0A_ax2.axhline(0%2C%20color%3D'gray'%2C%20alpha%3D0.3%2C%20linewidth%3D0.6)%0A_ax2.axvline(0%2C%20color%3D'gray'%2C%20alpha%3D0.3%2C%20linewidth%3D0.6)%0A_ax2.plot(%5B0%2C%20_x0%5D%2C%20%5B0%2C%20_y0%5D%2C%20color%3D'red'%2C%20linewidth%3D2.2%2C%20label%3D'phasor%20at%20t%3D0')%0A_ax2.scatter(%5B_x0%5D%2C%20%5B_y0%5D%2C%20color%3D'red'%2C%20s%3D50%2C%20zorder%3D5)%0A_phi_deg%20%3D%20np.degrees(_phi)%0A_theta1%2C%20_theta2%20%3D%20(0%2C%20_phi_deg)%20if%20_phi_deg%20%3E%3D%200%20else%20(_phi_deg%2C%200)%0A_arc%20%3D%20mpatches.Arc((0%2C%200)%2C%200.6%2C%200.6%2C%20theta1%3D_theta1%2C%20theta2%3D_theta2%2C%20color%3D'red'%2C%20linewidth%3D1.8)%0A_ax2.add_patch(_arc)%0A_ax2.annotate(f'%24%5C%5Ctheta%3D%7B_phi%3A.2f%7D%24'%2C%20xy%3D(0.42%20*%20np.cos(_phi%20%2F%202)%2C%200.42%20*%20np.sin(_phi%20%2F%202))%2C%20color%3D'red'%2C%20fontsize%3D12%2C%20ha%3D'center'%2C%20va%3D'center')%0A_ax2.set_xlim(-1.35%2C%201.35)%0A_ax2.set_ylim(-1.35%2C%201.35)%0A_ax2.set_aspect('equal')%0A_ax2.set_xlabel('Real(z)'%2C%20fontsize%3D14)%0A_ax2.set_ylabel('Imaginary(z)'%2C%20fontsize%3D14)%0A_ax2.set_title(f'Phase%20plane%20%20(%24%5C%5Ctheta%3D%7B_phi%3A.2f%7D%24%20rad)'%2C%20fontsize%3D16)%0A_ax2.legend(loc%3D'upper%20right'%2C%20fontsize%3D9)%0Aplt.colorbar(_lc%2C%20ax%3D_ax2%2C%20fraction%3D0.046%2C%20pad%3D0.04%2C%20label%3D'Time%20(s)')%0A_ax3%20%3D%20_fig.add_subplot(2%2C%202%2C%203)%0A_ax3.plot(np.arange(len(_time))%20%2F%20osc_sf%2C%20real(_z))%0A_ax3.set_xlabel('Time%20(sec)'%2C%20fontsize%3D14)%0A_ax3.set_ylabel('Real(z)'%2C%20fontsize%3D14)%0A_ax3.set_title('Real%20vs%20Time'%2C%20fontsize%3D16)%0A_ax4%20%3D%20_fig.add_subplot(2%2C%202%2C%204)%0A_ax4.plot(np.arange(len(_time))%20%2F%20osc_sf%2C%20imag(_z))%0A_ax4.set_xlabel('Time%20(sec)'%2C%20fontsize%3D14)%0A_ax4.set_ylabel('Imaginary(z)'%2C%20fontsize%3D14)%0A_ax4.set_title('Imaginary%20vs%20Time'%2C%20fontsize%3D16)%0Aplt.tight_layout()%0Aplt.close()%0Amo.vstack(%5Bmo.hstack(%5Bcx_freq_slider%2C%20cx_theta_slider%5D)%2C%20_fig%5D)

Recall: freq = [3, 10, 5, 15, 35], amplitude = [5, 15, 10, 5, 7]

Let's learn a few more important details about the DFT:

Convolution Theorem

Convolution in the time domain is the same as multiplication in the frequency domain. This means that time domain convolution computations can be performed much more efficiently in the frequency domain via simple multiplication. (The opposite is also true that multiplication in the time domain is the same as convolution in the frequency domain. Watch this Video for an overview of the convolution theorem and convolution in the frequency domain.

Let's prove it: