Python による実開口レーダーの実装

 

まずは，ライブラリをインポートします．

import numpy as np

import matplotlib.pyplot as plt

from tqdm import tqdm

続いて，周波数，時間，波長，サンプル数などのパラメータを設定します．

# Examples Parameters

freqency = 2 * (1E3 ** 2)

time_sec = 1E-8

wavelength = 0.025

num_sample = 1024 * 3

print(f'- Frequency: {freqency} Hz \n- Time: {time_sec} s \n- Wavelength: {wavelength} m')

times = np.linspace(0, time_sec, num=num_sample)

描画してみます．

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('Time')

plt.plot(times, label='Time')

plt.legend(loc='upper left',)

plt.xlabel('Sample')

plt.ylabel('Time (s)')

plt.savefig(f'section_03_001_Time.png', bbox_inches='tight')

plt.show()

図-1 経過時間

図-1の横軸はサンプリング，縦軸はサンプリングした時の経過時間です．続いて，送信する矩形の送信信号を生成して描画してみます．

pulse_sample = 128

pulse_times = np.linspace(0, pulse_sample, pulse_sample)

pulse = np.zeros_like(pulse_times, dtype=np.complex64)

pulse[16:128-16] = 1

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('Pulse Pattarn')

plt.plot(pulse.real, label='Real', alpha=0.75)

plt.plot(pulse.imag, label='Imaginary', alpha=0.75)

plt.legend()

plt.xlabel('Sample')

plt.ylabel('Amplitude')

plt.savefig(f'section_03_002_pulse_pattarn.png', bbox_inches='tight')

plt.show()

図-2  送信する信号

図-2の横軸はサンプリング，縦軸は信号の大きさです．この送信信号から受信信号を作成してみます．

points_pulse = [

256,

1024, 1024+256,

1024+1024, 1024+1024+128, 1024+1024+256,

]

respons_signal = np.zeros_like(times, dtype=np.complex64)

for point in points_pulse:

respons_signal[point:point + pulse_sample] += pulse

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('Recieved Signal')

plt.plot(respons_signal.real, label='Real', alpha=0.75)

plt.plot(respons_signal.imag, label='Imaginary', alpha=0.75)

plt.legend()

plt.xlabel('Sample')

plt.ylabel('Amplitude')

plt.savefig(f'section_03_003_recieved_signal.png', bbox_inches='tight')

plt.show()

図-3 受信信号

def impulse_response(respons_signal, pulse):

length_pulse, length_signal = len(pulse), len(respons_signal)

# padding

_respons_signal = np.zeros((length_pulse+length_signal), dtype=np.complex64)

_respons_signal[:length_signal] = respons_signal

# process the signal corration

corr = np.zeros_like(_respons_signal)

for i in tqdm(range(len(respons_signal))):

corr[i] = np.sum(np.dot(_respons_signal[i:i+length_pulse] ,np.conj(pulse)))

corr = np.abs(corr)

return corr[:length_signal]

impulse_response_signal = impulse_response(respons_signal, pulse)

plt.figure(figsize=(16, 4), dpi=200, facecolor='w', edgecolor='k')

plt.title('impulse Response')

plt.plot(impulse_response_signal ,)

plt.axhline(80, color='r', linestyle='--', label='Threshold')

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_004_impulse_response.png', bbox_inches='tight')

plt.show()

ノイズ

noises = np.random.normal(0, 1, len(impulse_response_signal))

noises = np.abs(noises)

plt.figure(figsize=(16, 4), dpi=200, facecolor='w', edgecolor='k')

plt.title('Random Noise')

plt.plot(noises,)

plt.xlabel('Sample')

plt.ylabel('Amplitude')

plt.savefig(f'section_03_005_noise.png', bbox_inches='tight')

plt.show()

impulse_response_signal_noise = impulse_response(respons_signal+ noises, pulse)

plt.figure(figsize=(16, 4), dpi=200, facecolor='w', edgecolor='k')

plt.title('impulse Response with Noise')

plt.plot(impulse_response_signal_noise,)

plt.axhline(150, color='r', linestyle='--', label='Threshold')

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_006_impulse_response_with_noise.png', bbox_inches='tight')

plt.show()

noises = np.random.normal(0, 1, len(impulse_response_signal))

noises = np.abs(noises)

plt.figure(figsize=(16, 4), dpi=200, facecolor='w', edgecolor='k')

plt.title('Random Noise')

plt.plot(noises,)

plt.xlabel('Sample')

plt.ylabel('Amplitude')

plt.savefig(f'section_03_005_noise.png', bbox_inches='tight')

plt.show()

impulse_response_signal_noise = impulse_response(respons_signal+ noises, pulse)

plt.figure(figsize=(16, 4), dpi=200, facecolor='w', edgecolor='k')

plt.title('impulse Response with Noise')

plt.plot(impulse_response_signal_noise,)

plt.axhline(150, color='r', linestyle='--', label='Threshold')

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_006_impulse_response_with_noise.png', bbox_inches='tight')

plt.show()

noises = np.random.normal(0, 1, len(impulse_response_signal))

ratio = 0.9 # SNR Ratio

impulse_response_signal_noise = impulse_response(respons_signal *ratio + (1 - ratio) * noises, pulse)

plt.figure(figsize=(16, 4), dpi=200, facecolor='w', edgecolor='k')

plt.title(f'impulse Response with Noise {1 - ratio:.2f}')

plt.plot(impulse_response_signal_noise,)

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_007_impulse_response_with_noise_ratio{ratio:.2f}.png', bbox_inches='tight')

plt.show()

noises = np.random.normal(0, 1, len(impulse_response_signal))

ratio = 0.5 # SNR Ratio

impulse_response_signal_noise = impulse_response(respons_signal *ratio + (1 - ratio) * noises, pulse)

plt.figure(figsize=(16, 4), dpi=200, facecolor='w', edgecolor='k')

plt.title(f'impulse Response with Noise {1 - ratio:.2f}')

plt.plot(impulse_response_signal_noise,)

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_007_impulse_response_with_noise_ratio{ratio:.2f}.png', bbox_inches='tight')

plt.show()

noises = np.random.normal(0, 1, len(impulse_response_signal))

ratio = 0.1 # SNR Ratio

impulse_response_signal_noise = impulse_response(respons_signal *ratio + (1 - ratio) * noises, pulse)

plt.figure(figsize=(16, 4), dpi=200, facecolor='w', edgecolor='k')

plt.title(f'impulse Response with Noise {1 - ratio:.2f}')

plt.plot(impulse_response_signal_noise,)

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_007_impulse_response_with_noise_ratio{ratio:.2f}.png', bbox_inches='tight')

plt.show()

Chirp Pulse

ratio = 0.25 # SNR Ratio

pulse_sample = 128

pulse_times = np.linspace(0, pulse_sample, pulse_sample)

pulse = np.exp(2j * np.pi * pulse_times)

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('Sin Pulse')

plt.plot(pulse.real, label='Real', alpha=0.75)

plt.plot(pulse.imag, label='Imaginary', alpha=0.75)

plt.legend()

plt.xlabel('Sample')

plt.ylabel('Amplitude')

plt.savefig(f'section_03_008_sin_pulse.png', bbox_inches='tight')

plt.show()

respons_signal = np.zeros_like(times, dtype=np.complex64)

for point in points_pulse:

respons_signal[point:point + pulse_sample] += pulse

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('Recieved Signal Sin')

plt.plot(respons_signal.real, label='Real', alpha=0.75)

plt.plot(respons_signal.imag, label='Imaginary', alpha=0.75)

plt.legend()

plt.xlabel('Sample')

plt.ylabel('Amplitude')

plt.savefig(f'section_03_009_recieved_signal_sin.png', bbox_inches='tight')

plt.show()

impulse_response_signal = impulse_response(respons_signal, pulse)

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('impulse Response')

plt.plot(impulse_response_signal,)

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_010_impulse_response_sin.png', bbox_inches='tight')

plt.show()

noises = np.random.normal(0, 1, len(impulse_response_signal))

impulse_response_signal_noise = impulse_response(respons_signal *ratio + (1 - ratio) * noises, pulse)

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title(f'impulse Response with Noise {1 - ratio:.2f}')

plt.plot(impulse_response_signal_noise,)

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_011_impulse_response_sin_with_noise_ratio{ratio:.2f}.png', bbox_inches='tight')

plt.show()

pulse_sample = 128

pulse_times = np.linspace(0, pulse_sample, pulse_sample)

pulse = np.exp(1j * np.pi * pulse_times)

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('Chirp Pulse')

plt.plot(pulse.real, label='Real', alpha=0.75)

plt.plot(pulse.imag, label='Imaginary', alpha=0.75)

plt.legend()

plt.xlabel('Sample')

plt.ylabel('Amplitude')

plt.savefig(f'section_03_012_chirp_pulse.png', bbox_inches='tight')

plt.show()

respons_signal = np.zeros_like(times, dtype=np.complex64)

for point in points_pulse:

respons_signal[point:point + pulse_sample] += pulse

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('Recieved Signal Chirp')

plt.plot(respons_signal.real, label='Real', alpha=0.75)

plt.plot(respons_signal.imag, label='Imaginary', alpha=0.75)

plt.legend()

plt.xlabel('Sample')

plt.ylabel('Amplitude')

plt.savefig(f'section_03_013_recieved_signal_chirp.png', bbox_inches='tight')

plt.show()

impulse_response_signal = impulse_response(respons_signal, pulse)

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('impulse Response')

plt.plot(impulse_response_signal,)

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_014_impulse_response_chirp.png', bbox_inches='tight')

plt.show()

impulse_response_signal_noise = impulse_response(respons_signal *ratio + (1 - ratio) * noises, pulse)

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title(f'impulse Response with Noise {1 - ratio:.2f}')

plt.plot(impulse_response_signal_noise,)

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_015_impulse_response_chirp_with_noise_ratio{ratio:.2f}.png', bbox_inches='tight')

plt.show()

respons_signal = np.zeros_like(times, dtype=np.complex64)

for point in points_pulse:

respons_signal[point:point + pulse_sample] += pulse

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('Recieved Signal Chirp')

plt.plot(respons_signal.real, label='Real', alpha=0.75)

plt.plot(respons_signal.imag, label='Imaginary', alpha=0.75)

plt.legend()

plt.xlabel('Sample')

plt.ylabel('Amplitude')

plt.savefig(f'section_03_017_recieved_signal_chirp.png', bbox_inches='tight')

plt.show()

impulse_response_signal = impulse_response(respons_signal, pulse)

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title('impulse Response')

plt.plot(impulse_response_signal,)

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_018_impulse_response_chirp.png', bbox_inches='tight')

plt.show()

impulse_response_signal_noise = impulse_response(respons_signal *ratio + (1 - ratio) * noises, pulse)

plt.figure(figsize=(16, 3), dpi=200, facecolor='w', edgecolor='k')

plt.title(f'impulse Response with Noise {1 - ratio:.2f}')

plt.plot(impulse_response_signal_noise,)

plt.xlabel('Sample')

plt.ylabel('Correlation')

plt.savefig(f'section_03_019_impulse_response_chirp_with_noise_ratio{ratio:.2f}.png', bbox_inches='tight')

plt.show()

Appendix

Numpy implimentation `impulse Response`

def impulse_response(respons_signal, pulse):

return np.convolve(respons_signal, pulse)

Movie

import os

import warnings, gc

from scipy.signal import chirp, spectrogram, correlate

warnings.simplefilter('ignore')

skip_step = 10

t_sec = 60

t_sec_chirp_start = 30

tau_p_sec = 10

fs_sampling = 1E2

K_r = 1.01

noise_top = 1.

SNR = 0.25

assert 0 < SNR < 1

tau = np.arange(-tau_p_sec/2, tau_p_sec/2, 1/fs_sampling)

phase = 1j*np.pi*K_r*tau**2

ra_chirp_temp = np.exp(phase)

slice_chirp = slice(int(t_sec_chirp_start*fs_sampling), int((t_sec_chirp_start+tau_p_sec)*fs_sampling))

T = np.arange(0, t_sec, 1/fs_sampling)

noise = np.random.uniform(low=-noise_top, high=noise_top, size=T.shape)

reflect_noise = noise.copy()

noise_tau_sec = np.random.uniform(low=-noise_top, high=noise_top, size=tau.shape)

corr_T = []

for i, t1 in tqdm(enumerate(T), total=len(T)):

# shortcut plot

if i % skip_step != 0:

continue

# pulse width time

if t1 > T.max() - tau_p_sec:

break

_slice_chirp = slice(round(t1*fs_sampling), round((t1+tau_p_sec)*fs_sampling))

# reference pulse

kernel = np.zeros(T.shape)

kernel[_slice_chirp] = ra_chirp_temp

# correlation

corr = correlate(kernel.real[_slice_chirp], reflect_noise.real[_slice_chirp], mode='same')

corr_center = corr[int(len(corr)/2)]

if i == 0:

corr_T.append(corr_center)

else:

corr_T.extend([corr_center]*skip_step)

# visualization

plt.figure(figsize=(26, 18), dpi=200, facecolor='w', edgecolor='k')

plt.subplot(6, 1, 1)

plt.grid(True,)

plt.plot(tau, corr, alpha=0.8, color='green', label='Correlation Value')

plt.scatter([0], corr_center, label='Target Center', color='orange')

plt.scatter([tau[np.argmax(corr)]], corr[np.argmax(corr)], label='Maximum Peak Correlation', color='black')

plt.legend()

plt.ylim(-50, 400)

plt.title(f'Local Correlation Time:{t1: .2f}')

plt.ylabel('Correlation')

plt.subplot(6, 1, 2)

plt.plot(T, kernel, color='orange', label='Kernel')

plt.grid(True,)

plt.title('Reference Pulse Chirp Signal')

plt.xlabel('Time axis [sec]')

plt.ylabel('Amplitude')

reflect_noise = noise.copy()

reflect_noise[slice_chirp] = ((1 - SNR)*ra_chirp_temp + SNR*noise_tau_sec)

plt.subplot(6, 1, 3)

plt.grid(True,)

plt.plot(T, reflect_noise.real, 'c')

plt.xlabel('Time axis [sec]')

plt.ylabel('Real')

plt.subplot(6, 1, 4)

plt.grid(True,)

plt.plot(T, kernel, 'navy', alpha=0.7, label='Kernel')

plt.plot(T, reflect_noise.real, 'magenta', alpha=0.4, label='Noise Reflect')

plt.legend()

plt.ylabel('Amplitude')

plt.subplot(6, 1, 5)

plt.grid(True,)

plt.plot(range(i+1), np.array(corr_T), 'b', alpha=0.8, label='Peak')

plt.scatter([i], corr_center, label='Center Correlation', color='red')

plt.legend()

plt.xlim(0, int(len(T) - len(ra_chirp_temp)))

plt.ylim(-50, 400)

plt.title('Correlation Value')

plt.xlabel('Sample [n]')

plt.ylabel('Correlation')

f, t, Sxx = spectrogram(reflect_noise, fs_sampling)

plt.subplot(6, 1, 6)

plt.grid(True,)

plt.pcolormesh(t, f, Sxx, shading='gouraud', label='spectrum', vmax=1e-1, cmap='inferno')

plt.vlines(t1 + round(tau_p_sec/2), ymin=0, ymax=50, color='red', label='position')

plt.legend()

plt.ylabel('Spectrum')

plt.xlabel('Time [sec]')

plt.tight_layout()

os.makedirs('frames', exist_ok=True)

plt.savefig(os.path.join('frames', f'impulse_respose_frame_{str(i).zfill(5)}.png'), format='png')

plt.clf()

plt.close()