This notebook holds the design parameters and generates an audio chirp


```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
#max free 8 points per sample

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

speed_of_sound_in_air_m_per_sec=343

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

print(f"The wavelength ranges from {(speed_of_sound_in_air_m_per_sec/f_start_Hz):.3f}m to {speed_of_sound_in_air_m_per_sec/f_stop_Hz:.3f} m")

#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=44.1e3

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

#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)
```

    The wavelength ranges from 0.114m to 0.017 m



```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_sec,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_sec,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_sec,chirp_phase_rad,linewidth=4,label='chirp');
ax[1].plot([t_sec[0], t_sec[-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 (sec)')
ax[1].set_ylabel('total phase (rad)')
ax[1].legend();


```


    <IPython.core.display.Javascript object>







```python
#Write a wav file that is 32-bit floating-point [-1.0,+1.0] np.float32
from scipy.io.wavfile import write
write(f'../data/audio_chirp_{t_sec}sec_{f_start_Hz}Hz_{f_stop_Hz}Hz.wav', int(fs), np.real(chirp).astype(np.float32))
```

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