52 lines
1.2 KiB
Python
52 lines
1.2 KiB
Python
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
from numpy import random
|
|
import scipy.signal as sgl
|
|
|
|
yi = random.normal(size=2**20)
|
|
yr = np.roll(yi, -(2**10))
|
|
|
|
figy, axy = plt.subplots()
|
|
axy.plot(np.arange(2**10, 2**11), yi[2**10 : 2**11])
|
|
axy.plot(np.arange(2**10, 2**11), yr[2**10 : 2**11])
|
|
|
|
figf, axf = plt.subplots()
|
|
axf.plot(*sgl.welch(yi))
|
|
axf.plot(*sgl.welch(yr))
|
|
|
|
eta = lambda r: yi + r * yr
|
|
u = lambda r: -yi + r * yr
|
|
|
|
|
|
def puv(eta, u):
|
|
f, phi_eta = sgl.welch(eta)
|
|
phi_u = sgl.welch(u)[1]
|
|
phi_eta_u = np.abs(sgl.csd(eta, u)[1].real)
|
|
|
|
return f, np.sqrt(
|
|
(phi_eta + phi_u - 2 * phi_eta_u) / (phi_eta + phi_u + 2 * phi_eta_u)
|
|
)
|
|
|
|
|
|
figr, axr = plt.subplots()
|
|
|
|
for r in np.arange(0, 1.1, 0.1):
|
|
axr.plot(*puv(eta(r), u(r)), c="k")
|
|
Rn = puv(
|
|
eta(r) + 0.4 * random.normal(size=2**20),
|
|
u(r) + 0.4 * random.normal(size=2**20),
|
|
)
|
|
axr.plot(
|
|
*Rn,
|
|
c="#ff6600",
|
|
)
|
|
axr.annotate(
|
|
f"{r=:.1f}", (Rn[0][0], Rn[1][0]), bbox={"boxstyle": "square", "facecolor": "w"}
|
|
)
|
|
axr.grid()
|
|
axr.autoscale(True, "x", tight=True)
|
|
axr.set(ylim=(0, 1), ylabel="R", xlabel="f")
|
|
axr.legend(("No noise", "40% noise"), loc="lower left")
|
|
|
|
figr.savefig("out_r_test.pdf")
|
|
plt.show()
|