# Interface
# レーダーの信号処理とデータサイエンス
#　演習１：FMCWレーダーのアナログ波形の生成
#　(c) 米本成人
#   July, 2023

import numpy as np
import matplotlib.pyplot as plt

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

N = 1024          # サンプル数
dt = 1/N          # サンプリング周期 [s]
time = np.linspace(0, 1, N)#時間データを作る
#noisesig = norm.rvs(loc=0, scale=1,size=N)  # 標準正規分布
mu, sigma = 0.0, 1 #ノイズの平均値、標準偏差

noisesig=np.random.normal(mu, sigma, N)#ノイズ

ave=np.mean(noisesig) #平均値
vari=np.var(noisesig) #分散
stdev=np.std(noisesig,ddof=0) #標準偏差

#統計データの表示
print("ガウス雑音")#データの種類の表示
print("電圧の平均値=",ave)#平均値
print("電圧の分散 =",vari)#分散
print("電圧の標準偏差=",stdev)#標準偏差

#図の描画（波形）
fig1,ax = plt.subplots()#figureウインドウを表示
ax.plot(time, noisesig, label = "ガウス雑音")#グラフの描画
ax.set_xlabel("時間[sec]")#x軸のラベル
ax.set_ylabel("観測電圧[V]")#y軸のラベル
ax.grid()#グラフの中にグリッド（網）を表示する
plt.legend() #凡例の表示
plt.savefig("./Fig1.png") #図の保存

#図の描画（度数分布）
fig2,ax = plt.subplots()#figureウインドウを表示
plt.hist(noisesig,bins=30, label="ガウス雑音の度数分布")
ax.set_xlabel("観測電圧[V]")#x軸のラベル
ax.set_ylabel("該当するデータの個数、頻度")#y軸のラベル
plt.legend() #凡例の表示
plt.savefig("./Fig2.png") #図の保存

#
# ノイズ信号にFFTを掛けてみる
#
fftdata = np.fft.fft(noisesig)/N #FFTを掛ける
freq = np.fft.fftfreq(N, d=dt) # 周波数
LinFFT = np.abs(fftdata) # 振幅

fig3,ax = plt.subplots()#figureウインドウを表示
ax.plot(freq, LinFFT, label = "ガウス雑音の周波数成分")#グラフの描画
ax.set_xlabel("周波数 [Hz]")#x軸のラベル
ax.set_ylabel("電圧[V]")#y軸のラベル
ax.grid()#グラフの中にグリッド（網）を表示する
plt.legend() #凡例の表示
plt.savefig("./Fig3.png") #図の保存


f1, f2 = 50, 100    # 周波数 [Hz]
p1, p2 = 1, 0.5 #振幅「V]
t = np.arange(0, 1, dt) # 時間 [s]
y = p1*np.sin(2*np.pi*f1*t) + p2*np.sin(2*np.pi*f2*t) +noisesig # さっきのノイズに信号2つ足す

ave=np.mean(y) #平均値
vari=np.var(y) #分散
stdev=np.std(y,ddof=0) #標準偏差

print("\n信号とガウス雑音") #タイトルっ表示
print("電圧の平均値=",ave)#平均値
print("電圧の分散 =",vari)#分散
print("電圧の標準偏差=",stdev)#標準偏差

#図の描画（波形）
fig4,ax = plt.subplots()#figureウインドウを表示
ax.plot(time, y, label = "信号の観測電圧")#グラフの描画
ax.set_xlabel("時間[sec]")#x軸のラベル
ax.set_ylabel("観測電圧[V]")#y軸のラベル
ax.grid()#グラフの中にグリッド（網）を表示する
plt.legend() #凡例の表示
plt.savefig("./Fig4.png") #図の保存

#図の描画（度数分布）
fig5,ax = plt.subplots()#figureウインドウを表示
plt.hist(y,bins=30, label="信号の度数分布")
ax.set_xlabel("観測電圧[V]")#x軸のラベル
ax.set_ylabel("該当するデータの個数、頻度")#y軸のラベル
plt.legend() #凡例の表示
plt.savefig("./Fig5.png") #図の保存


#
# FFTを掛けてみる
#
fftdata = np.fft.fft(y)/N #FFTを掛ける
freq = np.fft.fftfreq(N, d=dt) # 周波数
LinFFT = np.abs(fftdata) # 振幅

#図の描画
fig6,ax = plt.subplots()#figureウインドウを表示
ax.plot(freq, LinFFT, label = "信号の周波数成分")#グラフの描画
ax.set_xlabel("ベースバンド周波数 [Hz]")#x軸のラベル
ax.set_ylabel("振幅[V]")#y軸のラベル
ax.grid()#グラフの中にグリッド（網）を表示する
plt.legend() #凡例の表示
plt.savefig("./Fig6.png") #図の保存

#電力換算
LogAmp = 20*np.log10(LinFFT)-10*np.log10(50)#P＝V^2/Rを対数で計算

#図の描画
fig7,ax = plt.subplots()#figure Secを表示
ax.plot(freq[1:int(N/2)], LogAmp[1:int(N/2)], label = "信号の周波数成分")#グラフの描画
ax.set_xlabel("周波数[Hz]")#x軸のラベル
ax.set_ylabel("振幅dBm]")#y軸のラベル
ax.grid()#グラフの中にグリッド（網）を表示する
plt.legend() #凡例の表示
plt.savefig("./Fig7.png") #図の保存

#図の表示
plt.show()
