```python
%matplotlib notebook
```


```python
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from matplotlib.ticker import ScalarFormatter
import math
```

This notebook assumes you have completed the notebook [Introduction of sine waves](TDS_Introduction-sine_waves.ipynb). This notebook follows the same pattern of time domain waveform generation: instantaneous frequency -> phase step -> total phase -> time domain waveform. 

Our goal is to track features of different acoustic impedance in material using a low power time domain waveform. Time delay spectrometry (TDS) is one implementation of this goal. To understand TDS we need to understand the waveform which is used by TDS called a chirp. A chirp is a sinusoid that is constantly varying in frequency. The chirp is generated by integrating a varying phase step which is derived from an instantaneous frequency profile. We will generate a chirp in this notebook.
The phase of the chirp can be found by integrating the instantaneous frequency:

\begin{equation}
	f(t)=\frac{f_{end}-f_{start}}{T_c}t + f_{start}
\end{equation}

\begin{equation}
\Delta\phi(t) = 2\pi f(t)\Delta t
\end{equation}

\begin{equation}
\phi (t)=\int_{}^{} \Delta\phi(t) = \int_{}^{} 2\pi f(t) dt = \int_{}^{}\frac{f_{end}-f_{start}}{T_c}tdt + \int_{}^{}f_{start}dt
\end{equation}

\begin{equation}
 \phi (t)= \frac{f_{end}-f_{start}}{T_c}\int_{}^{}tdt + f_{start}\int_{}^{}dt
\end{equation}

\begin{equation}
 \phi (t)= \frac{f_{end}-f_{start}}{T_c}\frac{t^2}{2} + f_{start}t
\end{equation}

This gives the time series value of

\begin{equation}
x(t) = e^{j\phi (t)} = e^{j(\frac{f_{end}-f_{start}}{T_c}\frac{t^2}{2} + f_{start}t)} 
\end{equation}

But the formula for phase requires squaring time which will cause numeric errors as the time increases. Another approach is to implement the formula for phase as a cummulative summation. 

\begin{equation}
\phi_{sum} (N)=\sum_{k=1}^{N} \Delta\phi(k) = \sum_{k=1}^{N} 2\pi f(k) t_s = \sum_{k=1}^{N}(\frac{f_{end}-f_{start}}{T_c}k + f_{start})t_s
\end{equation}


This allow for the phase always stay between 0 and two pi by subtracting two phi whenever the phase exceeds the value. We will work with the cummlative sum of phase, but then compare it to the integral to determine how accurate the cummulative sum is.




```python
speed_of_ultrasound_m_per_sec={ 'tissue':1540, #In most tissue
                                'water':1481, 
                                'air': 343
                              }
```


```python
#Try out differment chirp paramters to determine which would produce the best seperation of distances for the gel pads
c_m_per_sec = speed_of_ultrasound_m_per_sec['water']


f_start_Hz=200_000
#talk about difference and similarity of sine wave example, answer why not 32 samples
f_stop_Hz=1_800_000

for Tc_sec in [0.01, 0.1, 1.0, 10.0]:
    S_Hz_per_sec = (f_stop_Hz - f_start_Hz)/Tc_sec
    d_m = 0.015 #1 cm
    tau_sec = 2 *d_m/c_m_per_sec
    f_location_Hz = S_Hz_per_sec * tau_sec
    print(f"Tc = {Tc_sec}. For a target at {d_m*100} cm the expected freq is {f_location_Hz} Hz with a time delay of {tau_sec} sec")
    
```

    Tc = 0.01. For a target at 1.5 cm the expected freq is 3241.053342336259 Hz with a time delay of 2.025658338960162e-05 sec
    Tc = 0.1. For a target at 1.5 cm the expected freq is 324.1053342336259 Hz with a time delay of 2.025658338960162e-05 sec
    Tc = 1.0. For a target at 1.5 cm the expected freq is 32.41053342336259 Hz with a time delay of 2.025658338960162e-05 sec
    Tc = 10.0. For a target at 1.5 cm the expected freq is 3.241053342336259 Hz with a time delay of 2.025658338960162e-05 sec



```python
#Tc is the max depth we are interested in
Tc_sec=1.0

f_start_Hz=200_000
#talk about difference and similarity of sine wave example, answer why not 32 samples
f_stop_Hz=1_800_000

#sample rate
fs=2_000_000
ts=1/fs

total_samples= math.ceil(fs*Tc_sec)
n = np.arange(0,total_samples, step=1, dtype=np.float64)
t_sec=n*ts
t_usec = t_sec *1e6

#This is the frequency of the chirp over time. We assume linear change in frequency
chirp_freq_slope_HzPerSec=(f_stop_Hz-f_start_Hz)/Tc_sec

#Compute the instantaneous frequency which is a linear function
chirp_instantaneous_freq_Hz=chirp_freq_slope_HzPerSec*t_sec+f_start_Hz
chirp_instantaneous_angular_freq_radPerSec=2*np.pi*chirp_instantaneous_freq_Hz

#Since frequency is a change in phase the we can plot it as a phase step
chirp_phase_step_rad=chirp_instantaneous_angular_freq_radPerSec*ts

#The phase step can be summed (or integrated) to produce the total phase which is the phase value 
#for each point in time for the chirp function
chirp_phase_rad=np.cumsum(chirp_phase_step_rad)
#The time domain chirp function
chirp = np.exp(1j*chirp_phase_rad)
```


```python
#We can see, unlike the complex exponential, the chirp's instantaneous frequency is linearly increasing. 
#This corresponds with the linearly increasing phase step. 
fig, ax = plt.subplots(2, 1, sharex=True,figsize = [8, 8])
lns1=ax[0].plot(t_usec,chirp_instantaneous_freq_Hz,linewidth=4, label='instantanous frequency');
ax[0].set_title('Comparing the instantaneous frequency and phase step')
ax[0].set_ylabel('instantaneous frequency (Hz)')

axt = ax[0].twinx()
lns2=axt.plot(t_usec,chirp_phase_step_rad,linewidth=2,color='black', linestyle=':', label='phase step');
axt.set_ylabel('phase step (rad)')

#ref: https://stackoverflow.com/questions/5484922/secondary-axis-with-twinx-how-to-add-to-legend
lns = lns1+lns2
labs = [l.get_label() for l in lns]
ax[0].legend(lns, labs, loc=0)

#We see that summing or integrating the linearly increasing phase step gives a quadratic function of total phase.
ax[1].plot(t_usec,chirp_phase_rad,linewidth=4,label='chirp');
ax[1].plot([t_usec[0], t_usec[-1]],[chirp_phase_rad[0], chirp_phase_rad[-1]],linewidth=1, linestyle=':',label='linear (x=y)');
ax[1].set_title('Cumulative quandratic phase function of chirp')
ax[1].set_xlabel('time ($\mu$sec)')
ax[1].set_ylabel('total phase (rad)')
ax[1].legend();


```


    <IPython.core.display.Javascript object>







```python
print(f'chirp_phase_step_rad[0] = {chirp_phase_step_rad[0]} rad and the end is chirp_phase_step_rad[-1] = {chirp_phase_step_rad[-1]}')
print(f'amplitude = {chirp_phase_step_rad[-1] - chirp_phase_step_rad[0]} rad')
print(f'offset = {(chirp_phase_step_rad[-1] - chirp_phase_step_rad[0])/2+chirp_phase_step_rad[0]} rad')
```

    chirp_phase_step_rad[0] = 0.6283185307179585 rad and the end is chirp_phase_step_rad[-1] = 5.654864263187504
    amplitude = 5.026545732469545 rad
    offset = 3.141591396952731 rad



```python
#We now what to load this into GNURadio. So we will only write the raw file and then load that into GNURadio. This will cause the array size to be lost
chirp.astype('complex64').tofile(f'../data/make_chirp_for_GNURadio_conf_fs_{int(fs)}_Tc_sec_{Tc_sec}_complex64.bin')
print(f'The sample rate is {int(fs)}')
```

    The sample rate is 2000000


be0bb4d58272679d813b037a8d309762b5f61f3d
ipynb
Jupyter Notebook
course/tds-200/week_01/notebooks/make_chirp_for_GNURadio_conf.ipynb
potto216/tds-tutorials
2acf2002ac5514dc60781c3e2e6797a4595104e6
MIT
2020-07-12T19:17:59.000Z
2020-09-24T22:19:02.000Z
course/tds-200/week_01/notebooks/make_chirp_for_GNURadio_conf.ipynb
potto216/tds-tutorials
2acf2002ac5514dc60781c3e2e6797a4595104e6
MIT
2020-09-16T12:18:01.000Z
2020-12-17T23:04:37.000Z
course/tds-200/week_01/notebooks/make_chirp_for_GNURadio_conf.ipynb
potto216/tds-tutorials
2acf2002ac5514dc60781c3e2e6797a4595104e6
MIT
Qwen/Qwen-72B
1. YES
2. YES