# Interface
# レーダーの信号処理とデータサイエンス
#　演習２：回転型レーダーのレーダー画像の生成
#　(c) 米本成人
#   July, 2023

import numpy as np
import matplotlib.pyplot as plt

#タイトルで漢字が使えるようフォントを設定
plt.rcParams['font.family'] = 'Meiryo'


#レーダーのパラメータ設定
MaxRange = 100#最大距離
ViewAngle = 30#m
N=MaxRange*2+2#サンプル数1メモリ1mとして、0からMaxRangeまで
dt = 1/(N-1)#サンプリング周期
M = 16 #アンテナ数
mu, sigma = 0.0, 1 #ノイズの平均値、標準偏差

antenna = np.linspace(1,M,M)#アンテナ番号
angle = np.pi*np.linspace(-ViewAngle/2, ViewAngle/2,M)/180#アンテナの画角設定（radian）
distance = np.linspace(0,N,N)#距離(m)
t = np.linspace(0, 1, N) # 時間 [s]

# 3D描画用のmeshgrid を作成
range_mesh, antenna_mesh = np.meshgrid(distance, antenna)#アンテナ番号と距離に合わせたメッシュデータを作る
range_mesh, angle_mesh = np.meshgrid(distance, angle)#角度と距離に合わせたメッシュデータを作る。
x_mesh = range_mesh * np.sin(angle_mesh)#レーダーからみて横方向の座標をX軸とする
y_mesh = range_mesh * np.cos(angle_mesh)#レーダーから見て距離方向の座標をY軸とする
x1, x2 = np.array_split(x_mesh, 2, 1)#メッシュデータを分割して必要なところだけ使う
y1, y2 = np.array_split(y_mesh, 2, 1)#メッシュデータを分割して必要なところだけ使う

#データ作成
#とりあえず空の行列を作る
noisedata = 0*range_mesh#ノイズデータ用変数
signaldata= 0*range_mesh#信号データ用変数
fftdata = 0j*range_mesh#FFT計算後のデータを入れる変数

#２つの物体からの反射波の周波数、振幅、および到来角を設定。
f1, f2 = 50, 90    # 周波数 [kHz]
a1, a2 = 10, 10        # 振幅
p1, p2 = -6/16, 2/16    # 位相
for m in range(M):
    noisedata[m]=np.random.normal(mu, sigma, N)#ノイズ

signaldata[2]=a1*np.cos(2*np.pi*(f1*t+p1*m)) # 信号1
signaldata[8]= a2*np.sin(2*np.pi*(f2*t+p2*m)) # 信号2

timedomaindata = signaldata+noisedata#M個のアンテナの時間波形

#データ描画
fig8 = plt.figure(8)
ax = fig8.add_subplot(projection='3d')
ax.set_title("測定時の時間領域波形") #グラフタイトルを設定
ax.set_ylabel("データ番号",color="black")#軸ラベル
ax.set_xlabel("サンプリング番号",color="black")#軸ラベル
ax.set_zlabel("信号強度",color="black")#軸ラベル
ax.plot_surface(antenna_mesh,range_mesh,timedomaindata,cmap="jet")#3次元的な曲面を描く
plt.savefig("./Fig8.png") #図の保存

#
#   各種信号処理
#
#   1DFFT　
for m in range(M):
     fftdata[m] = np.fft.fft(timedomaindata[m])/N #時間領域のデータを#1行ずつ1次元FFT

absdata = abs(fftdata)#振幅計算
magdata = 20*np.log10(absdata)#dBに変換
posdata, negadata = np.array_split(magdata, 2, 1)#ベースバンドを正の周波数と負の周波数に分割
mapdata = posdata#画像には正のでーだだけを仕様

fig9 = plt.figure(9)
ax = fig9.add_subplot(projection='3d')
ax.plot_surface(x1,y1,mapdata,cmap="jet")#レーダ画像
ax.set_title("測定値の周波数分析") #グラフタイトルを設定
ax.set_xlabel("x座標(m)",color="black")#軸ラベル
ax.set_ylabel("y座標(m)",color="black")#軸ラベル
ax.set_zlabel("信号強度",color="black")#軸ラベル
plt.savefig("./Fig9.png") #図の保存

#フィルタ作成
theta = np.linspace(0,1,M)#アンテナ毎のポイント
#filan = 0.5-0.5*np.cos(2*np.pi*theta)#空間方向にHanning window
filan = np.linspace(1,1,M)#時間方向の振幅今回は方形窓、
fildis = 0.5-0.5*np.cos(2*np.pi*distance/N)#時間方向にHanning window
f, a = np.meshgrid(fildis, filan) #フィルタの座標作成
filterdata = f*a#フィルタの係数計算

#フィルタの断面
fig10 = plt.figure(10)
plt.plot(fildis)
plt.savefig("./Fig10.png") #図の保存

#データに掛ける係数
fig11 = plt.figure(11)
ax = fig11.add_subplot(projection='3d')
ax.plot_surface(antenna_mesh,range_mesh,filterdata,cmap="jet")
ax.set_title("フィルタ形状")#グラフタイトルを設定
ax.set_ylabel("アンテナ番号",color="black")#軸ラベル
ax.set_xlabel("レンジ軸",color="black")#軸ラベル
ax.set_zlabel("信号強度",color="black")#軸ラベル
plt.savefig("./Fig11.png") #図の保存

#空間フィルタリング
filtereddata= filterdata*timedomaindata#フィルタを掛ける

fig12 = plt.figure(12)
ax = fig12.add_subplot(projection='3d')
ax.plot_surface(antenna_mesh,range_mesh,filtereddata,cmap="jet")
ax.set_title("フィルタリングした時間波形")#グラフタイトルを設定
ax.set_ylabel("アンテナの向き",color="black")#軸ラベル
ax.set_xlabel("サンプリング番号",color="black")#軸ラベル
ax.set_zlabel("信号強度",color="black")#軸ラベル
plt.savefig("./Fig12.png") #図の保存

for m in range(M):
     fftdata[m] = np.fft.fft(filtereddata[m])#1行ずつ1次元FFT
absdata = abs(fftdata)
magdata = 20*np.log10(absdata)
posdata, negadata = np.array_split(magdata, 2, 1)
mapdata = posdata


fig13 = plt.figure(13)
ax = fig13.add_subplot(projection='3d')
#ax.set_aspect('equal')
#曲面を描画
ax.plot_surface(x1,y1,mapdata,cmap="jet")
ax.set_title("生成したレーダー画像")#グラフタイトルを設定
ax.set_xlabel("x座標(m)",color="black")#軸ラベル
ax.set_ylabel("y座標(m)",color="black")#軸ラベル
ax.set_zlabel("信号強度",color="black")#軸ラベル
plt.savefig("./Fig13.png") #図の保存

#描画
plt.show()

