import numpy as np
import matplotlib.pyplot as plt

# レーダーパラメータ
c = 3e8                # 光速 [m/s]
fc = 10e9              # キャリア周波数 [Hz] (Xバンド)
lam = c / fc           # 波長 [m]
B = 500e6              # 帯域幅 [Hz]
delta_R = c / (2 * B)  # レンジ分解能 [m]

# パルスパラメータ
PRF = 500              # パルス繰り返し周波数 [Hz]
N_pulses = 256         # パルス数
T_obs = N_pulses / PRF # 観測時間 [s]

# 目標の回転
omega = np.deg2rad(3) / T_obs  # 観測時間中に3°回転
delta_theta = omega * T_obs     # 総回転角度 [rad]

print(f"波長: {lam*100:.1f} cm")
print(f"レンジ分解能: {delta_R:.2f} m")
print(f"クロスレンジ分解能: {lam/(2*delta_theta):.2f} m")
print(f"観測時間: {T_obs:.2f} s")
print(f"回転角度: {np.rad2deg(delta_theta):.1f}°")

# 目標モデル（船舶状の散乱点配置）
# 各散乱点: (x, y, amplitude)
# x: レンジ方向 [m], y: クロスレンジ方向 [m]
scatterers = [
    # 船体中心線
    (-30, 0, 1.0),
    (-20, 0, 0.8),
    (-10, 0, 0.6),
    (0, 0, 1.5),     # 船橋
    (10, 0, 0.7),
    (20, 0, 0.5),
    (30, 0, 0.4),
    # 上部構造
    (-5, 3, 0.9),
    (5, 3, 0.8),
    (-5, -3, 0.9),
    (5, -3, 0.8),
    (0, 5, 1.2),     # マスト
    (0, -5, 1.2),
    # 煙突
    (15, 2, 1.0),
    (15, -2, 1.0),
]

scatterers = np.array(scatterers)
x_s = scatterers[:, 0]  # レンジ方向位置
y_s = scatterers[:, 1]  # クロスレンジ方向位置
amp_s = scatterers[:, 2]  # 散乱振幅

# 基準距離
R0 = 50e3  # 50 km

# レンジビンの設定
N_range = 256
range_bins = np.arange(N_range) * delta_R - N_range / 2 * delta_R

# 信号シミュレーション
# 各パルス時刻
t_slow = np.arange(N_pulses) / PRF  # slow-time

# データ行列 (パルス × レンジビン)
signal = np.zeros((N_pulses, N_range), dtype=complex)

for m in range(N_pulses):
    t = t_slow[m]
    theta = omega * t  # 現在の回転角度

    for i in range(len(x_s)):
        # 回転後の散乱点位置
        x_rot = x_s[i] * np.cos(theta) - y_s[i] * np.sin(theta)
        y_rot = x_s[i] * np.sin(theta) + y_s[i] * np.cos(theta)

        # レーダーからの距離（レンジ方向の射影）
        R_p = R0 + x_rot

        # パルス圧縮後の信号（sinc関数で近似）
        range_diff = range_bins - x_rot
        # パルス圧縮応答
        pc_response = amp_s[i] * np.sinc(range_diff / delta_R)

        # 位相項（ドップラー情報を含む）
        phase = -4 * np.pi * R_p / lam
        signal[m, :] += pc_response * np.exp(1j * phase)

# 雑音追加
noise_level = 0.05
signal += noise_level * (np.random.randn(N_pulses, N_range) +
                          1j * np.random.randn(N_pulses, N_range))

print(f"\n信号行列サイズ: {signal.shape}")

import numpy as np
import matplotlib.pyplot as plt

# （上のコードの続き — パラメータ再掲）
c = 3e8
fc = 10e9
lam = c / fc
B = 500e6
delta_R = c / (2 * B)
PRF = 500
N_pulses = 256
T_obs = N_pulses / PRF
omega = np.deg2rad(3) / T_obs
delta_theta = omega * T_obs
N_range = 256
range_bins = np.arange(N_range) * delta_R - N_range / 2 * delta_R
t_slow = np.arange(N_pulses) / PRF

# 目標散乱点（船舶モデル）
scatterers_data = [
    (-30, 0, 1.0), (-20, 0, 0.8), (-10, 0, 0.6),
    (0, 0, 1.5), (10, 0, 0.7), (20, 0, 0.5), (30, 0, 0.4),
    (-5, 3, 0.9), (5, 3, 0.8), (-5, -3, 0.9), (5, -3, 0.8),
    (0, 5, 1.2), (0, -5, 1.2), (15, 2, 1.0), (15, -2, 1.0),
]
scatterers_arr = np.array(scatterers_data)
x_s, y_s, amp_s = scatterers_arr[:, 0], scatterers_arr[:, 1], scatterers_arr[:, 2]

R0 = 50e3
signal = np.zeros((N_pulses, N_range), dtype=complex)

for m in range(N_pulses):
    t = t_slow[m]
    theta = omega * t
    for i in range(len(x_s)):
        x_rot = x_s[i] * np.cos(theta) - y_s[i] * np.sin(theta)
        y_rot = x_s[i] * np.sin(theta) + y_s[i] * np.cos(theta)
        R_p = R0 + x_rot
        range_diff = range_bins - x_rot
        pc_response = amp_s[i] * np.sinc(range_diff / delta_R)
        phase = -4 * np.pi * R_p / lam
        signal[m, :] += pc_response * np.exp(1j * phase)

noise_level = 0.05
np.random.seed(0)
signal += noise_level * (np.random.randn(N_pulses, N_range) +
                          1j * np.random.randn(N_pulses, N_range))

# ---- 運動補償 ----

# 1. 包絡線アライメント（相互相関によるレンジシフト補正）
reference = signal[0, :]  # 基準パルス
aligned_signal = np.zeros_like(signal)
aligned_signal[0, :] = signal[0, :]

for m in range(1, N_pulses):
    # 相互相関（FFTベース）
    corr = np.fft.ifft(np.fft.fft(reference) * np.conj(np.fft.fft(signal[m, :])))
    shift = np.argmax(np.abs(corr))
    if shift > N_range // 2:
        shift -= N_range
    aligned_signal[m, :] = np.roll(signal[m, :], -shift)

# 2. 位相補償（支配散乱点法）
# 最も強い散乱点のレンジビンを特定
range_profile_mean = np.mean(np.abs(aligned_signal), axis=0)
dominant_bin = np.argmax(range_profile_mean)

# 支配散乱点の位相履歴を抽出
phase_history = np.angle(aligned_signal[:, dominant_bin])
# 位相補正
for m in range(N_pulses):
    aligned_signal[m, :] *= np.exp(-1j * phase_history[m])

# ---- ISAR画像生成（レンジ-ドップラー処理） ----

# Hamming窓の適用（クロスレンジ方向）
window = np.hamming(N_pulses)
windowed_signal = aligned_signal * window[:, np.newaxis]

# クロスレンジFFT
N_cross = 512  # ゼロパディング
isar_image = np.fft.fftshift(
    np.fft.fft(windowed_signal, n=N_cross, axis=0), axes=0
)

# 画像のパワー（dBスケール）
isar_power = 20 * np.log10(np.abs(isar_image) + 1e-12)
isar_power -= isar_power.max()

# 軸の計算
cross_range_res = lam / (2 * delta_theta)
doppler_axis = np.fft.fftshift(np.fft.fftfreq(N_cross, d=1/PRF))
cross_range_axis = doppler_axis * lam / (2 * omega)

# 表示範囲の設定
range_display = range_bins
cross_display = cross_range_axis

# ISAR画像の表示
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# 1. 散乱点の配置（真値）
axes[0].scatter(x_s, y_s, s=amp_s * 100, c='red', edgecolors='black', zorder=5)
for i in range(len(x_s)):
    axes[0].annotate(f'{amp_s[i]:.1f}', (x_s[i] + 1, y_s[i] + 0.5), fontsize=8)
axes[0].set_xlabel('Range [m]', fontsize=12)
axes[0].set_ylabel('Cross-Range [m]', fontsize=12)
axes[0].set_title('True Scatterer Locations', fontsize=13)
axes[0].set_xlim([-40, 40])
axes[0].set_ylim([-15, 15])
axes[0].set_aspect('equal')
axes[0].grid(True, alpha=0.3)

# 2. ISAR画像
extent = [range_bins[0], range_bins[-1],
          cross_range_axis[0], cross_range_axis[-1]]
im = axes[1].imshow(
    isar_power, aspect='auto', extent=extent,
    origin='lower', cmap='jet', vmin=-30, vmax=0
)
axes[1].set_xlabel('Range [m]', fontsize=12)
axes[1].set_ylabel('Cross-Range [m]', fontsize=12)
axes[1].set_title('ISAR Image', fontsize=13)
axes[1].set_xlim([-40, 40])
axes[1].set_ylim([-15, 15])
plt.colorbar(im, ax=axes[1], label='Power [dB]')

# 3. レンジプロファイル
range_profile = 20 * np.log10(np.mean(np.abs(aligned_signal), axis=0) + 1e-12)
range_profile -= range_profile.max()
axes[2].plot(range_bins, range_profile, 'b-', linewidth=1.5)
axes[2].set_xlabel('Range [m]', fontsize=12)
axes[2].set_ylabel('Power [dB]', fontsize=12)
axes[2].set_title('Average Range Profile', fontsize=13)
axes[2].set_xlim([-40, 40])
axes[2].set_ylim([-30, 5])
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"レンジ分解能: {delta_R:.2f} m")
print(f"クロスレンジ分解能: {cross_range_res:.2f} m")

import numpy as np
import matplotlib.pyplot as plt

c = 3e8
fc = 10e9
lam = c / fc
B = 500e6
delta_R = c / (2 * B)
PRF = 500
N_range = 256
range_bins = np.arange(N_range) * delta_R - N_range / 2 * delta_R
R0 = 50e3

# 散乱点（簡略化した目標）
scat = np.array([
    (-15, 0, 1.0), (0, 0, 1.5), (15, 0, 1.0),
    (0, 5, 1.2), (0, -5, 1.2),
    (-10, 3, 0.8), (10, 3, 0.8),
    (-10, -3, 0.8), (10, -3, 0.8),
])
x_s, y_s, amp_s = scat[:, 0], scat[:, 1], scat[:, 2]

# 異なる回転角度で比較
rotation_angles = [1.0, 3.0, 6.0, 10.0]  # [deg]
N_pulses = 256

fig, axes = plt.subplots(1, 4, figsize=(20, 5))

for idx, rot_deg in enumerate(rotation_angles):
    rot_rad = np.deg2rad(rot_deg)
    T_obs = N_pulses / PRF
    omega_val = rot_rad / T_obs
    t_slow = np.arange(N_pulses) / PRF

    # 信号生成
    sig = np.zeros((N_pulses, N_range), dtype=complex)
    for m in range(N_pulses):
        theta = omega_val * t_slow[m]
        for i in range(len(x_s)):
            x_rot = x_s[i] * np.cos(theta) - y_s[i] * np.sin(theta)
            R_p = R0 + x_rot
            range_diff = range_bins - x_rot
            pc_resp = amp_s[i] * np.sinc(range_diff / delta_R)
            phase = -4 * np.pi * R_p / lam
            sig[m, :] += pc_resp * np.exp(1j * phase)

    # 雑音
    np.random.seed(42)
    sig += 0.05 * (np.random.randn(N_pulses, N_range) +
                    1j * np.random.randn(N_pulses, N_range))

    # 運動補償（簡略版 - 位相補償のみ）
    rp_mean = np.mean(np.abs(sig), axis=0)
    dom_bin = np.argmax(rp_mean)
    ph = np.angle(sig[:, dom_bin])
    for m in range(N_pulses):
        sig[m, :] *= np.exp(-1j * ph[m])

    # ISAR画像
    window = np.hamming(N_pulses)
    sig_w = sig * window[:, np.newaxis]
    N_cross = 512
    isar = np.fft.fftshift(np.fft.fft(sig_w, n=N_cross, axis=0), axes=0)
    isar_dB = 20 * np.log10(np.abs(isar) + 1e-12)
    isar_dB -= isar_dB.max()

    # クロスレンジ軸
    doppler = np.fft.fftshift(np.fft.fftfreq(N_cross, d=1/PRF))
    cr_axis = doppler * lam / (2 * omega_val)

    extent = [range_bins[0], range_bins[-1], cr_axis[0], cr_axis[-1]]
    axes[idx].imshow(isar_dB, aspect='auto', extent=extent,
                     origin='lower', cmap='jet', vmin=-25, vmax=0)
    cr_res = lam / (2 * rot_rad)
    axes[idx].set_title(f'Δθ={rot_deg}°, Δy={cr_res:.2f}m', fontsize=12)
    axes[idx].set_xlabel('Range [m]', fontsize=11)
    axes[idx].set_ylabel('Cross-Range [m]', fontsize=11)
    axes[idx].set_xlim([-25, 25])
    axes[idx].set_ylim([-15, 15])

plt.suptitle('ISAR Image Quality vs Rotation Angle', fontsize=14)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

c = 3e8
fc = 10e9
lam = c / fc
B = 500e6
delta_R = c / (2 * B)
PRF = 500
N_pulses = 256
N_range = 256
range_bins = np.arange(N_range) * delta_R - N_range / 2 * delta_R
R0 = 50e3

# 2つの散乱点（クロスレンジ方向に接近）
scat = np.array([
    (0, 0, 1.0),
    (0, 2, 0.3),  # 近くに弱い散乱点
])
x_s, y_s, amp_s = scat[:, 0], scat[:, 1], scat[:, 2]

rot_deg = 3.0
rot_rad = np.deg2rad(rot_deg)
T_obs = N_pulses / PRF
omega_val = rot_rad / T_obs
t_slow = np.arange(N_pulses) / PRF

# 信号生成
sig = np.zeros((N_pulses, N_range), dtype=complex)
for m in range(N_pulses):
    theta = omega_val * t_slow[m]
    for i in range(len(x_s)):
        x_rot = x_s[i] * np.cos(theta) - y_s[i] * np.sin(theta)
        R_p = R0 + x_rot
        range_diff = range_bins - x_rot
        pc_resp = amp_s[i] * np.sinc(range_diff / delta_R)
        phase = -4 * np.pi * R_p / lam
        sig[m, :] += pc_resp * np.exp(1j * phase)

np.random.seed(0)
sig += 0.02 * (np.random.randn(N_pulses, N_range) +
                1j * np.random.randn(N_pulses, N_range))

# 位相補償
rp_mean = np.mean(np.abs(sig), axis=0)
dom_bin = np.argmax(rp_mean)
ph = np.angle(sig[:, dom_bin])
for m in range(N_pulses):
    sig[m, :] *= np.exp(-1j * ph[m])

# 窓関数の比較
windows = {
    'Rectangular': np.ones(N_pulses),
    'Hamming': np.hamming(N_pulses),
    'Hanning': np.hanning(N_pulses),
    'Blackman': np.blackman(N_pulses),
}

N_cross = 1024
doppler = np.fft.fftshift(np.fft.fftfreq(N_cross, d=1/PRF))
cr_axis = doppler * lam / (2 * omega_val)

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

# クロスレンジプロファイル（レンジビン = 支配散乱点）
for name, win in windows.items():
    sig_w = sig * win[:, np.newaxis]
    isar = np.fft.fftshift(np.fft.fft(sig_w, n=N_cross, axis=0), axes=0)
    cross_profile = np.abs(isar[:, dom_bin])
    cross_dB = 20 * np.log10(cross_profile / cross_profile.max() + 1e-12)
    axes[0].plot(cr_axis, cross_dB, linewidth=1.5, label=name)

axes[0].set_xlabel('Cross-Range [m]', fontsize=13)
axes[0].set_ylabel('Power [dB]', fontsize=13)
axes[0].set_title('Cross-Range Profile at Dominant Range Bin', fontsize=13)
axes[0].set_xlim([-10, 10])
axes[0].set_ylim([-60, 5])
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# 真の散乱点位置をマーク
axes[0].axvline(x=0, color='gray', linestyle=':', alpha=0.5)
axes[0].axvline(x=2, color='gray', linestyle=':', alpha=0.5)

# Hamming窓でのISAR画像
sig_h = sig * np.hamming(N_pulses)[:, np.newaxis]
isar_h = np.fft.fftshift(np.fft.fft(sig_h, n=N_cross, axis=0), axes=0)
isar_dB = 20 * np.log10(np.abs(isar_h) + 1e-12)
isar_dB -= isar_dB.max()

extent = [range_bins[0], range_bins[-1], cr_axis[0], cr_axis[-1]]
im = axes[1].imshow(isar_dB, aspect='auto', extent=extent,
                     origin='lower', cmap='jet', vmin=-30, vmax=0)
axes[1].set_xlabel('Range [m]', fontsize=12)
axes[1].set_ylabel('Cross-Range [m]', fontsize=12)
axes[1].set_title('ISAR Image (Hamming Window)', fontsize=13)
axes[1].set_xlim([-10, 10])
axes[1].set_ylim([-10, 10])
plt.colorbar(im, ax=axes[1], label='Power [dB]')

plt.tight_layout()
plt.show()