パルスドップラーレーダーの原理を解説して実装する

import numpy as np
import matplotlib.pyplot as plt

# =============================================================
# レーダーパラメータ
# =============================================================
c = 3e8            # 光速 [m/s]
f0 = 10e9          # 搬送波周波数 [Hz] (Xバンド)
lam = c / f0       # 波長 [m]
B = 5e6            # LFM帯域幅 [Hz]
Tp = 10e-6         # パルス幅 [s]
PRF = 1000         # PRF [Hz]
T_PRI = 1 / PRF    # PRI [s]
fs_fast = 20e6     # Fast-timeサンプリング周波数 [Hz]
M = 64             # パルス数
Pt = 1e3           # 送信電力 [W]（シミュレーション上の概念値）

# 導出されるパラメータ
delta_R = c / (2 * B)  # 距離分解能
R_ua = c * T_PRI / 2   # 最大無曖昧距離
v_ua = lam * PRF / 4   # 最大無曖昧速度
delta_v = lam * PRF / (2 * M)  # 速度分解能

print('=== Radar Parameters ===')
print(f'Carrier frequency: {f0/1e9:.1f} GHz')
print(f'Wavelength: {lam*100:.2f} cm')
print(f'LFM bandwidth: {B/1e6:.1f} MHz')
print(f'Pulse width: {Tp*1e6:.1f} μs')
print(f'PRF: {PRF} Hz')
print(f'Number of pulses: {M}')
print(f'Range resolution: {delta_R:.1f} m')
print(f'Max unambiguous range: {R_ua/1e3:.1f} km')
print(f'Max unambiguous velocity: {v_ua:.1f} m/s')
print(f'Velocity resolution: {delta_v:.2f} m/s')

# =============================================================
# 目標の定義
# =============================================================
targets = [
    {'range': 30e3, 'velocity': 50, 'rcs': 10},    # 目標1: 30km, 50m/s
    {'range': 50e3, 'velocity': -100, 'rcs': 5},    # 目標2: 50km, -100m/s
    {'range': 30e3, 'velocity': -30, 'rcs': 3},     # 目標3: 30km, -30m/s（目標1と同距離）
    {'range': 70e3, 'velocity': 0, 'rcs': 20},      # 目標4: 70km, 静止（クラッタ）
]

# =============================================================
# LFMパルスの生成
# =============================================================
t_pulse = np.arange(0, Tp, 1/fs_fast)  # パルス内の時間
N_pulse = len(t_pulse)

# LFMチャープ信号（ベースバンド）
chirp_rate = B / Tp
tx_pulse = np.exp(1j * np.pi * chirp_rate * (t_pulse - Tp/2)**2)

# マッチドフィルタ（送信パルスの時間反転複素共役）
matched_filter = np.conj(tx_pulse[::-1])

# =============================================================
# Fast-timeサンプル数（レンジ方向）
# =============================================================
R_max_sim = 100e3  # シミュレーション最大距離
tau_max = 2 * R_max_sim / c
N_fast = int(np.ceil((tau_max + Tp) * fs_fast))

# =============================================================
# 受信信号の生成
# =============================================================
np.random.seed(42)
rx_data = np.zeros((M, N_fast), dtype=complex)

for m in range(M):
    for tgt in targets:
        R = tgt['range'] + tgt['velocity'] * m * T_PRI
        tau = 2 * R / c  # 遅延時間
        fd = 2 * tgt['velocity'] / lam  # ドップラー周波数

        # 遅延サンプル数
        delay_samples = int(np.round(tau * fs_fast))

        if delay_samples + N_pulse <= N_fast and delay_samples >= 0:
            # 受信パルス = 遅延・ドップラーシフトされた送信パルス
            amplitude = np.sqrt(tgt['rcs'])
            phase_shift = np.exp(1j * 2 * np.pi * fd * m * T_PRI)
            rx_data[m, delay_samples:delay_samples+N_pulse] += (
                amplitude * phase_shift * tx_pulse
            )

    # 受信ノイズ
    noise_power = 0.1
    rx_data[m] += np.sqrt(noise_power/2) * (
        np.random.randn(N_fast) + 1j * np.random.randn(N_fast)
    )

# =============================================================
# パルス圧縮（各パルスにマッチドフィルタ適用）
# =============================================================
N_conv = N_fast + N_pulse - 1
pc_data = np.zeros((M, N_conv), dtype=complex)

for m in range(M):
    pc_data[m] = np.convolve(rx_data[m], matched_filter, mode='full')

# パルス圧縮後のデータをレンジ方向にトリム
pc_data = pc_data[:, N_pulse-1:N_fast+N_pulse-1]

# レンジ軸
range_axis = np.arange(pc_data.shape[1]) * c / (2 * fs_fast)

# =============================================================
# ドップラー処理（各レンジゲートでFFT）
# =============================================================
# 窓関数を適用（サイドローブ低減）
window = np.hamming(M)
rd_map = np.zeros_like(pc_data)

for k in range(pc_data.shape[1]):
    rd_map[:, k] = np.fft.fftshift(np.fft.fft(pc_data[:, k] * window))

# ドップラー軸
doppler_axis = np.fft.fftshift(np.fft.fftfreq(M, T_PRI))
velocity_axis = doppler_axis * lam / 2

# =============================================================
# Range-Dopplerマップの表示
# =============================================================
rd_power = np.abs(rd_map)**2
rd_dB = 10 * np.log10(rd_power / np.max(rd_power) + 1e-12)

fig, axes = plt.subplots(2, 2, figsize=(16, 12))

# Range-Dopplerマップ
im = axes[0, 0].pcolormesh(
    range_axis / 1e3, velocity_axis,
    rd_dB, cmap='jet', vmin=-40, vmax=0, shading='auto'
)
plt.colorbar(im, ax=axes[0, 0], label='Power (dB)')
for tgt in targets:
    axes[0, 0].plot(tgt['range']/1e3, tgt['velocity'], 'wo',
                    markersize=10, markeredgecolor='k', markeredgewidth=2)
axes[0, 0].set_xlabel('Range (km)')
axes[0, 0].set_ylabel('Velocity (m/s)')
axes[0, 0].set_title('Range-Doppler Map')
axes[0, 0].set_xlim([0, R_max_sim/1e3])

# パルス圧縮結果（1パルス目）
axes[0, 1].plot(range_axis / 1e3,
                20*np.log10(np.abs(pc_data[0])/np.max(np.abs(pc_data[0]))+1e-12),
                'b', linewidth=1)
for tgt in targets:
    axes[0, 1].axvline(x=tgt['range']/1e3, color='r', linestyle='--', alpha=0.5)
axes[0, 1].set_xlabel('Range (km)')
axes[0, 1].set_ylabel('Amplitude (dB)')
axes[0, 1].set_title('Pulse Compression Output (Pulse #0)')
axes[0, 1].set_xlim([0, R_max_sim/1e3])
axes[0, 1].set_ylim([-40, 5])
axes[0, 1].grid(True)

# ドップラースペクトル（目標1のレンジゲート）
tgt1_range_bin = int(np.round(targets[0]['range'] * 2 / c * fs_fast))
doppler_spectrum = 20 * np.log10(
    np.abs(rd_map[:, tgt1_range_bin])
    / np.max(np.abs(rd_map[:, tgt1_range_bin])) + 1e-12
)
axes[1, 0].plot(velocity_axis, doppler_spectrum, 'b', linewidth=2)
axes[1, 0].axvline(x=targets[0]['velocity'], color='r', linestyle='--',
                    label=f"Target 1: {targets[0]['velocity']} m/s")
axes[1, 0].axvline(x=targets[2]['velocity'], color='g', linestyle='--',
                    label=f"Target 3: {targets[2]['velocity']} m/s")
axes[1, 0].set_xlabel('Velocity (m/s)')
axes[1, 0].set_ylabel('Amplitude (dB)')
axes[1, 0].set_title(f"Doppler Spectrum at R={targets[0]['range']/1e3:.0f} km")
axes[1, 0].legend()
axes[1, 0].grid(True)
axes[1, 0].set_ylim([-40, 5])

# MTIの効果
# 簡易2パルスMTIフィルタ
mti_data = pc_data[1:] - pc_data[:-1]
mti_power = np.mean(np.abs(mti_data)**2, axis=0)
no_mti_power = np.mean(np.abs(pc_data)**2, axis=0)

axes[1, 1].plot(range_axis/1e3,
                10*np.log10(no_mti_power/np.max(no_mti_power)+1e-12),
                'b', alpha=0.5, label='Without MTI')
axes[1, 1].plot(range_axis/1e3,
                10*np.log10(mti_power/np.max(no_mti_power)+1e-12),
                'r', linewidth=2, label='With MTI')
for tgt in targets:
    if tgt['velocity'] != 0:
        axes[1, 1].axvline(x=tgt['range']/1e3, color='g', linestyle=':', alpha=0.5)
axes[1, 1].set_xlabel('Range (km)')
axes[1, 1].set_ylabel('Power (dB)')
axes[1, 1].set_title('MTI Effect: Clutter Suppression')
axes[1, 1].set_xlim([0, R_max_sim/1e3])
axes[1, 1].set_ylim([-50, 5])
axes[1, 1].legend()
axes[1, 1].grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

c = 3e8
f0 = 10e9
lam = c / f0

# PRF範囲
PRF_range = np.logspace(2, 5, 500)  # 100 Hz ~ 100 kHz

# 最大無曖昧距離と速度
R_ua = c / (2 * PRF_range)
v_ua = lam * PRF_range / 4

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# R_ua vs PRF
axes[0].loglog(PRF_range, R_ua / 1e3, 'b', linewidth=2)
axes[0].set_xlabel('PRF (Hz)')
axes[0].set_ylabel('R_ua (km)')
axes[0].set_title('Max Unambiguous Range vs PRF')
axes[0].grid(True, which='both', alpha=0.3)
axes[0].axhspan(0, 10, alpha=0.1, color='red', label='Short range (<10 km)')
axes[0].axhspan(10, 100, alpha=0.1, color='yellow', label='Medium range')
axes[0].axhspan(100, 2000, alpha=0.1, color='green', label='Long range')
axes[0].legend(fontsize=8)

# v_ua vs PRF
axes[1].loglog(PRF_range, v_ua, 'r', linewidth=2)
axes[1].set_xlabel('PRF (Hz)')
axes[1].set_ylabel('v_ua (m/s)')
axes[1].set_title('Max Unambiguous Velocity vs PRF')
axes[1].grid(True, which='both', alpha=0.3)

# R_ua × v_ua = const のプロット
axes[2].loglog(R_ua/1e3, v_ua, 'g', linewidth=2)
product = c * lam / 8
axes[2].set_xlabel('R_ua (km)')
axes[2].set_ylabel('v_ua (m/s)')
axes[2].set_title(f'R_ua × v_ua = cλ/8 = {product:.0f} m²/s (X-band)')
axes[2].grid(True, which='both', alpha=0.3)

# PRFの代表的な値をプロット
for prf_val, label in [(500, 'Low'), (5000, 'Medium'), (50000, 'High')]:
    r = c / (2*prf_val) / 1e3
    v = lam * prf_val / 4
    axes[2].plot(r, v, 'ko', markersize=10)
    axes[2].annotate(f'{label} PRF\n({prf_val} Hz)',
                     xy=(r, v), xytext=(r*1.5, v*1.5),
                     arrowprops=dict(arrowstyle='->', color='black'),
                     fontsize=9)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

PRF = 1000
f0 = 10e9
lam = 3e8 / f0

# 正規化周波数
f_norm = np.linspace(-0.5, 0.5, 1000)
f = f_norm * PRF  # 実際のドップラー周波数
v = f * lam / 2   # 対応する速度

# 1段MTI: H(z) = 1 - z^{-1}
H1 = 2 * np.abs(np.sin(np.pi * f_norm))

# 2段MTI: H(z) = (1 - z^{-1})^2
H2 = 4 * np.sin(np.pi * f_norm)**2

# 3段MTI: H(z) = (1 - z^{-1})^3
H3 = 8 * np.abs(np.sin(np.pi * f_norm))**3

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 振幅特性（リニア）
axes[0].plot(v, H1/np.max(H1), 'b', linewidth=2, label='Single canceller')
axes[0].plot(v, H2/np.max(H2), 'r', linewidth=2, label='Double canceller')
axes[0].plot(v, H3/np.max(H3), 'g', linewidth=2, label='Triple canceller')
axes[0].set_xlabel('Velocity (m/s)')
axes[0].set_ylabel('Normalized Amplitude')
axes[0].set_title('MTI Filter Frequency Response')
axes[0].legend()
axes[0].grid(True)

# 振幅特性（dB）
H1_dB = 20 * np.log10(H1/np.max(H1) + 1e-12)
H2_dB = 20 * np.log10(H2/np.max(H2) + 1e-12)
H3_dB = 20 * np.log10(H3/np.max(H3) + 1e-12)

axes[1].plot(v, H1_dB, 'b', linewidth=2, label='Single canceller')
axes[1].plot(v, H2_dB, 'r', linewidth=2, label='Double canceller')
axes[1].plot(v, H3_dB, 'g', linewidth=2, label='Triple canceller')
axes[1].set_xlabel('Velocity (m/s)')
axes[1].set_ylabel('Amplitude (dB)')
axes[1].set_title('MTI Filter Response (dB)')
axes[1].set_ylim([-60, 5])
axes[1].legend()
axes[1].grid(True)

# ブラインドスピードをマーク
v_blind = lam * PRF / 2
for ax in axes:
    for n in range(1, 3):
        ax.axvline(x=n*v_blind/2, color='orange', linestyle=':', alpha=0.5)
        ax.axvline(x=-n*v_blind/2, color='orange', linestyle=':', alpha=0.5)

plt.tight_layout()
plt.show()

print(f'Blind speed (1st): {v_blind:.1f} m/s = {v_blind*3.6:.1f} km/h')