m2cce/internship
m2cce/internship
1
Fork 0
Code Issues Activity

PUV test with multiple reflection coefficients

Browse source
This commit is contained in:
Edgar P. Burkhart 2022-04-04 10:59:17 +02:00
parent 873cc73c84
commit 9e29a30b42
Signed by: edpibu
GPG key ID: 9833D3C5A25BD227
1 changed files with 15 additions and 16 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

31
swash/processing/r_test.py
View file

@ -4,7 +4,7 @@ from numpy import random
import scipy.signal as sgl import scipy.signal as sgl
yi = random.normal(size=2**20) yi = random.normal(size=2**20)
yr = 0.3 * np.roll(yi, -(2**10)) yr = np.roll(yi, -(2**10))
figy, axy = plt.subplots() 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), yi[2**10 : 2**11])
@ -14,8 +14,8 @@ figf, axf = plt.subplots()
axf.plot(*sgl.welch(yi)) axf.plot(*sgl.welch(yi))
axf.plot(*sgl.welch(yr)) axf.plot(*sgl.welch(yr))
eta = yi + yr eta = lambda r: yi + r * yr
u = -yi + yr u = lambda r: -yi + r * yr
def puv(eta, u): def puv(eta, u):
@ -30,21 +30,20 @@ def puv(eta, u):
figr, axr = plt.subplots() figr, axr = plt.subplots()
axr.plot(*puv(eta, u), label="Without noise") for r in np.arange(0, 1.1, 0.1):
axr.plot( axr.plot(*puv(eta(r), u(r)), c="k")
*puv( Rn = puv(
eta + 0.4 * random.normal(size=2**20), u + 0.4 * random.normal(size=2**20) eta(r) + 0.4 * random.normal(size=2**20),
), u(r) + 0.4 * random.normal(size=2**20),
label="With noise" )
) axr.plot(
*Rn,
c="#ff6600",
)
axr.annotate(f"{r=:.1f}", (Rn[0][0], Rn[1][0]), bbox={"boxstyle": "square", "facecolor": "w"})
axr.grid() axr.grid()
axr.autoscale(True, "x", tight=True) axr.autoscale(True, "x", tight=True)
axr.set(ylim=(0, 1), ylabel="R", xlabel="f") axr.set(ylim=(0, 1), ylabel="R", xlabel="f")
axr.legend(loc="lower left") axr.legend(("No noise", "40% noise"), loc="lower left")
axpsd = axr.twinx()
axpsd.plot(*sgl.welch(eta), label=r"PSD ($\eta$)", c="k", alpha=0.2, lw=1)
axpsd.legend(loc="lower right")
axpsd.set(ylabel="PSD")
plt.show() plt.show()