# Radar FMCW and CSM algorithms

This notebook shows a naive (didactically intuitive) implementation of following algorithms:

- FMCW (Frequency Modulated Continuous Wave)
- CSM (Chirp Sequence Modulation)

Both algorithms are used for range and velocity measurements in automated/assisted driving domain.

## Introduction

### Doppler

An electromagnetic wave 

\begin{align}
w(t) &= a_t cos(2\pi f_0t)
\end{align}

reflected on a moving object with velocity $v$ leads to a new wave:

\begin{align}
u_r(t) &= a_r \text{cos}(2\pi (f_0 - f_{d}) t + \phi_r + \phi_t) \\

\end{align}

In order to understand what signal the algorithms operate on we have to take a closer look at the inner workings of a continuous wave radar (simplified view):



[Source](https://en.wikipedia.org/wiki/Continuous-wave_radar#/media/File:Bsp2_CW-Radar.EN.png)

### Antennas

The CW-RADAR has one transmitting antenna and at least one receiving antenna. The former emits a continuous radio wave, that's different to pulse-based RADAR systems. The receiving antenna transforms the reflected waves from the object of interests and converts the signal into eletrical domain (voltage).

Automotive long range radars use a $77.7\text{GHz}$ base frequency.

### RF-Generator

The RF-generator is a module which dynamically generates a sine voltage with a particular frequency. Since both methods presented here modulate (alter) the frequency of the tramitted radio wave, the generator is an essential part in the system.

### Mixer

The mixer is the most important component in the system described here. It allows us to measure the frequency change introduced by the moving object we are going to detect. We apply to the antenna the voltage $u_t(t)$ and receive the voltage $u_r(t)$.

\begin{align}
u_t(t) &= A_t \text{cos}(2\pi f_0 t+\phi_t) \\
u_r(t) &= A_r \text{cos}(2\pi (f_0 - f_{d}) t + \phi_r)
\end{align}

With Doppler frequency

\begin{align}
f_{d} = \frac{2v}{\lambda} = \frac{2 \dot r}{\lambda} = \frac{2 f_0 \dot r}{c}
\end{align}

Given base wavelength $\lambda$ and target velocity $v=\dot r$.

Theoretically, if we could sample $f_0 - f_{d} \approx 77.7GHz$ and we were done. But high frequency AD converters (according to Nyquist $\geq 2f_0 \approx 145GHz$ is required) are techinically expensive and the change introduced by $f_d$ is very small to detect changes. That's why we need the mixer to get rid of the high frequency $f_0$ and work only on the $f_d$. 

The mixer multiplies both signals, corresponding to trigonometry relations:

\begin{align}
\text{cos}(x)\text{cos}(y) = \frac{1}{2} [ \text{cos}(x-y) + \text{cos}(x+y)]
\end{align}

Applied to $u_t(t)$ and $u_r(t)$:

\begin{align}
u_t(t)u_r(t) &= \frac{A_t A_r}{2} [ \text{cos}(2\pi f_0 t+\phi_t-2\pi (f_0 - f_{d}) t - \phi_r) + \text{cos}(2\pi f_0 t+\phi_t+ 2\pi (f_0 - f_{d}) t + \phi_r)] \\
&= \frac{A_t A_r}{2} [ \text{cos}(2\pi f_{d} t + \phi_t - \phi_r) + \text{cos}(2\pi f_0 t+\phi_t + 2\pi (f_0- f_{d}) t + \phi_r)] \\
&= \frac{A_t A_r}{2} [ \text{cos}(2\pi f_{d} t + \phi_t - \phi_r) + \text{cos}(2\pi \underbrace{2 f_0}_{\text{2x 77GHz}} t - 2\pi f_{d}t + \phi_t + \phi_r)] \\
\end{align}

The product of both signals equals to a sum of two waves, the first, only dependent on the doppler frequency (which is velocity dependent) and the second a very high frequency signal (144 GHz) which we filter out with the low pass filter:

\begin{align}
u_{It,r}(t) = \frac{A_t A_r}{2} \text{cos}(2\pi f_{d} t + \phi_t - \phi_r)
\end{align}

After the AD converter we could sample $u_{It,r}(t)$ and measure the (unsigned) velocity. The index $I$ means in-phase. We later will introduce an other index in order to extract the sign of the velocity (approaching vs. departing).


```python
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.fftpack import fft
```

Constant for velocity of light within the air:


```python
c = 299792458
```

## FCMW

Instead of sending with one frequency, we will alter $f_0$ over time:

\begin{align}
u_t(t) &= A_t \text{cos}(2\pi f_t t+\phi_t)
\end{align}

with 

\begin{align}
f_t = f_t(t) = f_0 + m_w(t-t_0)
\end{align}


```python
f_0 = 77.7*1e9 # 77.7GHz

phi_t = 0
A_t = 1
f_ramp = 425*1e6 #Hz
T_ramp = 0.010
m_w = f_ramp/T_ramp
```


```python
def f_t(t):
    return f_0 + m_w*(t % T_ramp)
```


```python
def u_t(t):
    f_t = f_t(t)
    return A_t*np.cos(2*np.pi*f_b*t + phi_t)
```

The received signal depends on both the range $r$ and the velocity $v$:

\begin{align}
u_r(t) &= A_r \text{cos}(2\pi f_r t+\phi_r)
\end{align}

with 

\begin{align}
f_r = f_r(t) = f_0 + m_w(t- \frac{2r}{c}-t_0) - \frac{2vf_0}{c}
\end{align}


```python
v = 50/3.6 # 50km/h
r = 100 # 100m distance
```


```python
def f_r(t):
    return f_0 + m_w*(t % T_ramp - 2*r/c) - 2*v*f_0/c
```


```python
t = np.arange(0, 2*T_ramp, 1e-8)
```


```python
plt.figure(figsize=(15,5))
plt.plot(t, f_t(t), label="$f_t$")
plt.plot(t, f_r(t), label="$f_r$")
plt.legend()
plt.xlabel("t [s]")
plt.ylabel("frequency [Hz]")
plt.title("Transmitted and received frequencies")
plt.grid();
```

Because we don't see any difference let's zoom in 


```python
t_op = T_ramp
plt.figure(figsize=(10,10))
plt.plot(t, f_t(t), label="$f_t$")
plt.plot(t, f_r(t), label="$f_r$")
plt.legend()
plt.xlim([t_op-T_ramp/10000, t_op+T_ramp/10000])
plt.ylim([f_0-100e3, f_0+100e3])
plt.xlabel("t [s]")
plt.ylabel("frequency [Hz]")
plt.title("Transmitted and received frequencies")
plt.grid();
```

If we could measure $\Delta f$ directly we could derive $r$ and $v$:

\begin{align}
\Delta f = f_t-f_r = \frac{2m_w}{c}r + \frac{2 f_0}{c}\dot r
\end{align}

There is just one problem. We have one equation but two unknowns. The solution for that problem is to have a second ramp with a different $m_w$:

\begin{align}
f_t(t) = \{\begin{array}{lr}
        f_0 + m_{w,1}(t-t_0), & \text{for } 0 \leq t < T_{r,1}\\
        f_0 + m_w T_{r,1} + m_{w,2}(t-t_0), & \text{for } T_{r,1} \leq t < T_{r,1}+T_{r,2}\
        \end{array}
\end{align}

Define and implement both waveforms:


```python
T_ramp1 = T_ramp*0.7
T_ramp2 = T_ramp*0.3
m_w1 = f_ramp/T_ramp1
m_w2 = -f_ramp/T_ramp2
```


```python
def f_t(t):
    r1 = f_0 + m_w1*(t % (T_ramp1+T_ramp2))
    r1[(t%(T_ramp1+T_ramp2)) > T_ramp1] = 0
    
    r2 = f_0 + m_w1*T_ramp1 + m_w2*( (t-T_ramp1) % (T_ramp1+T_ramp2))
    r2[(t%(T_ramp1+T_ramp2)) <= T_ramp1] = 0
    
    return r1+r2
```


```python
def f_r(t):
    r1 = f_0 + m_w1*(t % (T_ramp1+T_ramp2) - 2*r/c) - 2*v*f_0/c
    r1[(t%(T_ramp1+T_ramp2)) > T_ramp1] = 0
    
    r2 = f_0 + m_w1*T_ramp1 + m_w2*( (t-T_ramp1) % (T_ramp1+T_ramp2) - 2*r/c) - 2*v*f_0/c
    r2[(t%(T_ramp1+T_ramp2)) <= T_ramp1] = 0
    return r1+r2
```


```python
t = np.arange(0, 2*T_ramp, 1e-7)
plt.figure(figsize=(15,5))
plt.plot(t, f_t(t), label="$f_t$")
plt.plot(t, f_r(t), label="$f_r$")
plt.legend()
plt.xlabel("t [s]")
plt.ylabel("frequency [Hz]")
plt.title("FMCW ramps with different $m_{w,1} \\neq m_{w,2}$")
plt.grid();
```

Zoom in again:


```python
t_op = T_ramp
plt.figure(figsize=(10,5))
plt.plot(t, f_t(t), label="$f_t$")
plt.plot(t, f_r(t), label="$f_r$")
plt.legend()
plt.xlim([t_op-T_ramp/5000, t_op+T_ramp/5000])
plt.ylim([f_0-100e3, f_0+600e3])
plt.xlabel("t [s]")
plt.ylabel("frequency [Hz]")
plt.title("FMCW ramps with different $m_{w,1} \\neq m_{w,2}$")
plt.grid();
```

Now we have two equations:

\begin{align}
\Delta f_1 = \frac{2m_{w,1}}{c}r + \frac{2 f_0}{c}\dot r \\
\Delta f_2 = \frac{2m_{w,2}}{c}r + \frac{2 f_0}{c}\dot r
\end{align}

Solving the equation system:


```python
t_f1 = np.asarray([T_ramp*0.1]) # beginning of rising ramp
t_f2 = np.asarray([T_ramp*0.8])  # beginning of falling ramp
delta_f_1 = f_t(t_f1) - f_r(t_f1)
delta_f_2 = f_t(t_f2) - f_r(t_f2)
```


```python
A = np.asarray([[2*m_w1/c, 2*f_0/c], [2*m_w2/c, 2*f_0/c]])
Y = np.asarray([delta_f_1, delta_f_2])
```


```python
x = np.linalg.solve(A, Y).flatten()
print("range r:", x[0])
print("velocity v:", x[1]*3.6)
```

    range r: 99.99999999693145
    velocity v: 49.99999988508472


As you could see we could successfully estimate the range and velocity.

**Note**: we get $\Delta f_i$ for free as a result of the FFT on the mixed and filtered signal 


\begin{align}
u_{It,r}(t) = \frac{A_t A_r}{2} \text{cos}(2\pi (f_t - f_r) t + \phi_t - \phi_r) = \frac{A_t A_r}{2} \text{cos}(2\pi \Delta f t + \phi_t - \phi_r)
\end{align}

In the FMCW chapter we did not address various problems like:

- multi target estimation (what happens if the FFT returns 2, 3, 4, ... major peaks?)
- angle estimation

We will solve the first problem in our second method.

## Chirp Sequence Modulation

This method is an extension of FMCW. Instead of sending one singe rising and falling frequency chirp we will repeat that process $n_r$ times in a sequence with chirp frequency $f_{chirp}$:


```python
# chirp sequence frequency
f_chirp = 50*1e3 #Hz

# ramp frequency
f_r = 200*1e6 #Hz
T_r = 1/f_chirp # duration of one cycle
m_w = f_r/T_r

n_r = 150 # number of chirps
T_M = T_r*n_r
```

Because we perform analog sampling, we have to configure our AD-converter:


```python
# sample settings
f_s = 50e6 #50 MHz
n_s = int(T_r*f_s)
```

Base frequency setup


```python
f_0 = 77.7*1e9

# some helpful
w_0 = 2*np.pi*f_0
lambda_0 = c/f_0
```


```python
def f_transmitted(t):
    return f_0 + m_w*(t%T_r)
```


```python
def chirp(t):
    return np.cos(2*np.pi*(f_transmitted(t))*t)
```

Lets visualize the chirp sequence:


```python
t = np.linspace(0, 3*T_r, 1e6)
```


```python
plt.figure(figsize=(15,5))
plt.plot(t, f_transmitted(t))
plt.xlabel("t [s]")
plt.ylabel("frequency [Hz]")
plt.title("Chirp sequence Modulation, transmitted signal $f_t(t)$");
```

### Reflector setup 

This is our target of interest


```python
r_0 = 50 # initial distance
v_veh = 36/3.6 # velocity
```


```python
def get_range(t):
    return r_0+v_veh*t
```

### Returned waveform

According to Eq. 17.46 in [1](https://www.springer.com/de/book/9783658057336) the returned waveform (after mixing and LF-filter) is:

\begin{align}
u_{It,r}(t) = \frac{A_t A_r}{2} \text{cos}\left(2\pi(\frac{2m_w}{c}r + \frac{2 f_0}{c}\dot r)t + \frac{4 \pi f_0 r}{c} + 2\pi(\frac{2 r}{c})^2m_w\right)
\end{align}


```python
4*np.pi*100*f_0/c, 2*np.pi*(2*100/c)**2*m_w
```




    (325694.31641129137, 27.963945936914552)



As you can see $2\pi(\frac{2 r}{c})^2m_w$ has not so much influence on the phase, so we ignore it:

\begin{align}
u_{It,r}(t) = \frac{A_t A_r}{2} \text{cos}\left(2\pi(\frac{2m_w}{c}r + \frac{2 f_0}{c}\dot r)t + \frac{4 \pi f_0 r}{c}\right)
\end{align}


```python
def itr(t):
    r = get_range(t)
    w_itr = 2*f_0*v_veh/c + 2*m_w*r/c
    # we do t%T_r because the eq. above only valid within the ramp
    v = np.cos(2*np.pi*w_itr*(t%T_r) +2*r*2*np.pi*f_0/c)
    return v
```

We build up a table of $n_r × n_s$ where $n_r$ is the number of chirps (ramps) and $n_s$ is the number of samples within a chirp.


```python
print(n_r, n_s)
```

    150 1000



```python
t_sample = np.linspace(0, T_M, n_r*n_s)
```


```python
v_sample = itr(t_sample)
```


```python
plt.figure(figsize=(15,5))
plt.plot(t_sample, v_sample, "+")
plt.xlim(0, 0.1*T_r)
plt.xlabel("t [s]")
plt.title("Samples visualized in time range [0, 10% $T_r$]");
```

We allocate memory for the sampling data


```python
table = np.zeros((n_r, n_s))
```

Fill in the values ...


```python
for chirp_nr in range(n_r):
    table[chirp_nr, :] = v_sample[(chirp_nr*n_s):(n_s*(chirp_nr+1))]
```

DF for pretty printing


```python
table_df = pd.DataFrame(data=table, 
                        columns=["sample_%03d"%i for i in range(n_s)], 
                        index=["chirp_%03d"%i for i in range(n_r)])
```


```python
table_df.head(10)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>sample_000</th>
      <th>sample_001</th>
      <th>sample_002</th>
      <th>sample_003</th>
      <th>sample_004</th>
      <th>sample_005</th>
      <th>sample_006</th>
      <th>sample_007</th>
      <th>sample_008</th>
      <th>sample_009</th>
      <th>...</th>
      <th>sample_990</th>
      <th>sample_991</th>
      <th>sample_992</th>
      <th>sample_993</th>
      <th>sample_994</th>
      <th>sample_995</th>
      <th>sample_996</th>
      <th>sample_997</th>
      <th>sample_998</th>
      <th>sample_999</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>chirp_000</th>
      <td>0.905353</td>
      <td>0.999836</td>
      <td>0.920138</td>
      <td>0.680144</td>
      <td>0.321662</td>
      <td>-0.092857</td>
      <td>-0.491199</td>
      <td>-0.803969</td>
      <td>-0.976681</td>
      <td>-0.979245</td>
      <td>...</td>
      <td>0.415170</td>
      <td>0.007652</td>
      <td>-0.401199</td>
      <td>-0.740157</td>
      <td>-0.950170</td>
      <td>-0.994652</td>
      <td>-0.865854</td>
      <td>-0.586214</td>
      <td>-0.204448</td>
      <td>0.212936</td>
    </tr>
    <tr>
      <th>chirp_001</th>
      <td>0.976847</td>
      <td>0.804429</td>
      <td>0.491871</td>
      <td>0.093624</td>
      <td>-0.320933</td>
      <td>-0.679581</td>
      <td>-0.919838</td>
      <td>-0.999850</td>
      <td>-0.905677</td>
      <td>-0.653725</td>
      <td>...</td>
      <td>-0.225759</td>
      <td>-0.603753</td>
      <td>-0.876566</td>
      <td>-0.996670</td>
      <td>-0.943140</td>
      <td>-0.725302</td>
      <td>-0.381107</td>
      <td>0.029482</td>
      <td>0.434935</td>
      <td>0.764616</td>
    </tr>
    <tr>
      <th>chirp_002</th>
      <td>0.644981</td>
      <td>0.276856</td>
      <td>-0.139501</td>
      <td>-0.531555</td>
      <td>-0.831006</td>
      <td>-0.985685</td>
      <td>-0.968646</td>
      <td>-0.782858</td>
      <td>-0.460686</td>
      <td>-0.058257</td>
      <td>...</td>
      <td>-0.773011</td>
      <td>-0.964635</td>
      <td>-0.988205</td>
      <td>-0.839616</td>
      <td>-0.544753</td>
      <td>-0.154986</td>
      <td>0.261782</td>
      <td>0.632944</td>
      <td>0.893837</td>
      <td>0.999011</td>
    </tr>
    <tr>
      <th>chirp_003</th>
      <td>0.046790</td>
      <td>-0.365037</td>
      <td>-0.713270</td>
      <td>-0.937241</td>
      <td>-0.997932</td>
      <td>-0.884770</td>
      <td>-0.617469</td>
      <td>-0.242597</td>
      <td>0.174539</td>
      <td>0.561268</td>
      <td>...</td>
      <td>-0.999508</td>
      <td>-0.925246</td>
      <td>-0.689792</td>
      <td>-0.334165</td>
      <td>0.079680</td>
      <td>0.479643</td>
      <td>0.796044</td>
      <td>0.973761</td>
      <td>0.981834</td>
      <td>0.818855</td>
    </tr>
    <tr>
      <th>chirp_004</th>
      <td>-0.570722</td>
      <td>-0.856199</td>
      <td>-0.992513</td>
      <td>-0.955916</td>
      <td>-0.752785</td>
      <td>-0.418507</td>
      <td>-0.011319</td>
      <td>0.397840</td>
      <td>0.737690</td>
      <td>0.949024</td>
      <td>...</td>
      <td>-0.811267</td>
      <td>-0.501932</td>
      <td>-0.105152</td>
      <td>0.309948</td>
      <td>0.671049</td>
      <td>0.915242</td>
      <td>0.999983</td>
      <td>0.910510</td>
      <td>0.662410</td>
      <td>0.298907</td>
    </tr>
    <tr>
      <th>chirp_005</th>
      <td>-0.952571</td>
      <td>-0.993817</td>
      <td>-0.861923</td>
      <td>-0.579868</td>
      <td>-0.196791</td>
      <td>0.220570</td>
      <td>0.599505</td>
      <td>0.873996</td>
      <td>0.996222</td>
      <td>0.944891</td>
      <td>...</td>
      <td>-0.286397</td>
      <td>0.129656</td>
      <td>0.523121</td>
      <td>0.825448</td>
      <td>0.983966</td>
      <td>0.971059</td>
      <td>0.788975</td>
      <td>0.469436</td>
      <td>0.068113</td>
      <td>-0.345076</td>
    </tr>
    <tr>
      <th>chirp_006</th>
      <td>-0.941085</td>
      <td>-0.721066</td>
      <td>-0.375424</td>
      <td>0.035623</td>
      <td>0.440464</td>
      <td>0.768568</td>
      <td>0.962774</td>
      <td>0.989248</td>
      <td>0.843378</td>
      <td>0.550576</td>
      <td>...</td>
      <td>0.357312</td>
      <td>0.707445</td>
      <td>0.934326</td>
      <td>0.998429</td>
      <td>0.888586</td>
      <td>0.623932</td>
      <td>0.250578</td>
      <td>-0.166433</td>
      <td>-0.554448</td>
      <td>-0.845866</td>
    </tr>
    <tr>
      <th>chirp_007</th>
      <td>-0.541005</td>
      <td>-0.150571</td>
      <td>0.266096</td>
      <td>0.636404</td>
      <td>0.895838</td>
      <td>0.999200</td>
      <td>0.928483</td>
      <td>0.696006</td>
      <td>0.342272</td>
      <td>-0.071093</td>
      <td>...</td>
      <td>0.852756</td>
      <td>0.991682</td>
      <td>0.957836</td>
      <td>0.757114</td>
      <td>0.424486</td>
      <td>0.017904</td>
      <td>-0.391798</td>
      <td>-0.733240</td>
      <td>-0.946936</td>
      <td>-0.995655</td>
    </tr>
    <tr>
      <th>chirp_008</th>
      <td>0.082467</td>
      <td>0.482098</td>
      <td>0.797739</td>
      <td>0.974397</td>
      <td>0.981296</td>
      <td>0.817232</td>
      <td>0.510791</td>
      <td>0.115359</td>
      <td>-0.300171</td>
      <td>-0.663405</td>
      <td>...</td>
      <td>0.994356</td>
      <td>0.864427</td>
      <td>0.583895</td>
      <td>0.201635</td>
      <td>-0.215755</td>
      <td>-0.595555</td>
      <td>-0.871595</td>
      <td>-0.995784</td>
      <td>-0.946484</td>
      <td>-0.732285</td>
    </tr>
    <tr>
      <th>chirp_009</th>
      <td>0.671886</td>
      <td>0.915698</td>
      <td>0.999976</td>
      <td>0.910036</td>
      <td>0.661548</td>
      <td>0.297805</td>
      <td>-0.117823</td>
      <td>-0.512924</td>
      <td>-0.818662</td>
      <td>-0.981771</td>
      <td>...</td>
      <td>0.723355</td>
      <td>0.378482</td>
      <td>-0.032331</td>
      <td>-0.437512</td>
      <td>-0.766468</td>
      <td>-0.961886</td>
      <td>-0.989720</td>
      <td>-0.845120</td>
      <td>-0.553281</td>
      <td>-0.165046</td>
    </tr>
  </tbody>
</table>
<p>10 rows × 1000 columns</p>
</div>



Now our table consists all sampled measurement during a cycle. We will use this matrix to calculate range and velocity.

### FFT over the first chirp

For didactic reasons we do the FFT procedure only for the first churp to make an example of the procedure.


```python
chirp0_samples = table_df.iloc[0].values
```


```python
chirp0_magnitude = fft(chirp0_samples)
```


```python
# frequencies found by FFT, will be used later
frequencies = np.arange(0, n_s//2)*f_s/n_s
```


```python
frequencies[:3]
```




    array([     0.,  50000., 100000.])



Each frequency is referred to a range

\begin{align}
f_* = \frac{2m_w}{c}r + \frac{2 f_0}{c}\dot r
\end{align}

Because $m_w$ is so large, the influence of the velocity (second term) is very low:


```python
f_star1 = 2*m_w/c*100
f_star2 = 2*f_0/c*33
print(f_star1, f_star2)
print(f_star2/f_star1*100, "%")
```

    6671281.9039630415 17105.833929951634
    0.25640999999999997 %


That means:

\begin{align}
f_* \approx \frac{2m_w}{c}r
\end{align}

We rearrange:

\begin{align}
r \approx \frac{f_* c}{2 m_w}
\end{align}


```python
def freq_to_range(f):
    return f*c/(2*m_w)
```


```python
ranges = freq_to_range(frequencies)
```


```python
plt.figure(figsize=(10,5))
plt.plot(ranges, 2.0/n_s*np.abs(chirp0_magnitude[0:n_s//2]))
plt.plot(ranges, 2.0/n_s*np.abs(chirp0_magnitude[0:n_s//2]), "k+")
plt.xlabel("range $r$ [m]")
plt.title("derived range (FFT of chirp_0)")
print(freq_to_range(frequencies)[np.argmax(2.0/n_s*np.abs(chirp0_magnitude[0:n_s//2]))])
```

Looks like we got the distance right: $\approx 50m$

# Calculate range bins for each chirp

Now we do the same for each of $n_r$ chirps and save the FFT results in a table of shape $n_r × n_s/2$.


```python
range_table = np.zeros((n_r, n_s//2), dtype=np.csingle)
```


```python
for chirp_nr in range(n_r):
    chirp_ad_values = table_df.iloc[chirp_nr].values
    chirp_fft = fft(chirp_ad_values) # FFT
    range_table[chirp_nr, :] = 2.0/n_s*chirp_fft[:n_s//2]
```

Let's visualize the table


```python
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(15,10), sharex=True, sharey=True)
abs_axes = ax[0, 0]
phi_axes = ax[0, 1]
real_axes = ax[1, 0]
imag_axes = ax[1, 1]

im_asb = abs_axes.imshow(np.abs(range_table), cmap = plt.get_cmap('RdYlBu'))
abs_axes.set_xticks(range(ranges.size)[::50])
abs_axes.set_xticklabels(ranges[::50], rotation=90)
fig.colorbar(im_asb, ax=abs_axes)
abs_axes.set_xlabel("range [m]")
abs_axes.set_ylabel("chirp number")
abs_axes.set_title("$|A(j\omega)|$")

im_phi = phi_axes.imshow(np.angle(range_table)*360/(2*np.pi), cmap = plt.get_cmap('RdYlBu'))
fig.colorbar(im_phi, ax=phi_axes)
phi_axes.set_xlabel("range [m]")
phi_axes.set_ylabel("chirp number")
phi_axes.set_title("$∠ A(j\omega)$")
phi_axes.set_xticks(range(ranges.size)[::50])
phi_axes.set_xticklabels(ranges[::50], rotation=90)

im_real = real_axes.imshow(np.real(range_table), cmap = plt.get_cmap('RdYlBu'))
fig.colorbar(im_real, ax=real_axes)
real_axes.set_xlabel("range [m]")
real_axes.set_ylabel("chirp number")
real_axes.set_title("Real{$A(j\omega)$}")
real_axes.set_xticks(range(ranges.size)[::50])
real_axes.set_xticklabels(ranges[::50], rotation=90)

im_imag = imag_axes.imshow(np.imag(range_table), cmap = plt.get_cmap('RdYlBu'))
fig.colorbar(im_imag, ax=imag_axes)
imag_axes.set_xlabel("range [m]")
imag_axes.set_ylabel("chirp number")
imag_axes.set_title("Imag{$A(j\omega)$}");
imag_axes.set_xticks(range(ranges.size)[::50])
imag_axes.set_xticklabels(ranges[::50], rotation=90);

fig.suptitle("Range FFT table visualized.");
```

## Velocity estimation

We make a second FFT over each range bin (column).

Again, we initialize an empty table. Now the rows will be different velocity ranges.


```python
velocity_table = np.zeros((n_r, range_table.shape[1]), dtype=np.csingle)
```


```python
for r in range(range_table.shape[1]):
    range_bin_magn = range_table[:, r]
    range_bin_fft = fft(range_bin_magn)
    velocity_table[:, r]= 2.0/n_r*range_bin_fft
```

Second FFT on columns returns the phase shift

\begin{align}
... + \frac{4 \pi f_0 }{c} r = ... + \frac{4 \pi f_0 }{c} vt \\
\Rightarrow \frac{4 \pi f_0 }{c}v = \omega_{\text{FFT, columns}}
\end{align}

After rearranging ...


```python
def angle_freq_to_velocity(w):
    return w*c/(4*np.pi*f_0)
```


```python
omega_second = 2*np.pi*np.concatenate((np.arange(0, n_r//2), np.arange(-n_r//2, 0)[::-1]))*f_chirp/n_r
```


```python
velocities = angle_freq_to_velocity(omega_second)
```


```python
plt.figure(figsize=(15,10))
plt.imshow(np.abs(velocity_table), cmap = plt.get_cmap('RdYlBu'))
plt.xticks(range(ranges.size)[::20], ranges[::20], rotation=90);
plt.yticks(range(velocities.size)[::10], velocities[::10]);
plt.xlim([0, 200])
plt.xlabel("range $r$ [m]")
plt.ylabel("velocity $\\dot r = v$ [m/s]");
plt.title("Chirp Sequence Modulation Result - $r, \\dot r$ map")
plt.colorbar();
```

As you can see, we could successfully recognize (blue dot) the target in 50m range with velocity of $10 m/s$.


```python

```

7bb173bd02ae1356b1a8b4a37319026f2e918dd1
ipynb
Jupyter Notebook
RADAR.ipynb
kopytjuk/fmcw
ceba9f71e41d54c3b339c7e40a840a3d8db542d8
MIT
2019-12-23T05:06:19.000Z
2022-02-22T17:19:01.000Z
RADAR.ipynb
kopytjuk/fmcw
ceba9f71e41d54c3b339c7e40a840a3d8db542d8
MIT
RADAR.ipynb
kopytjuk/fmcw
ceba9f71e41d54c3b339c7e40a840a3d8db542d8
MIT
2020-05-06T20:54:58.000Z
2022-02-13T09:42:35.000Z
Qwen/Qwen-72B
1. YES
2. YES