```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 -> angle step -> total angle -> 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 angle step which is derived from an instantaneous frequency profile. We will generate a chirp in this notebook. An overview of this technique is given [here](https://www.youtube.com/watch?v=RQplkt0bw_c).
The angle 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 angle requires squaring time which will cause numeric errors as the time increases. Another approach is to implement the formula for angle 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 angle always stay between 0 and two pi by subtracting two phi whenever the angle exceeds the value. We will work with the cummlative sum of angle, but then compare it to the integral to determine how accurate the cummulative sum is.




```python

#Tc is the max depth we are interested in
Tc_sec=5
#Tc_sec=2

f_start_Hz=10e3
#talk about difference and similarity of sine wave example, answer why not 32 samples
f_stop_Hz=990e3

#We choose 8 samples per cycle at the maximum frequency to not require steep pulse shaping filter profiles on the output of the 
#digital to analog converter
#fs=16e6
fs=2e6
samplesPerCycle=fs/f_stop_Hz

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
#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/TDS_Part_1_chirp_compare_fs_{int(fs)}_samplesPerCycle_{int(samplesPerCycle)}_Tc_sec_{Tc_sec}_complex64.bin')
print(f'The sample rate is {int(fs)}')
print(f'The samples per cycle is {samplesPerCycle}')
```

    The sample rate is 2000000
    The samples per cycle is 2.0202020202020203



```python

```

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