import
numpy as np
import
matplotlib.pyplot as plt
def
gen_test_data(periods, noise
=
0
, rotary
=
False
, npts
=
1000
, dt
=
1.0
/
24
):
t
=
np.arange(npts, dtype
=
float
)
*
dt
if
rotary:
x
=
noise
*
(np.random.randn(npts)
+
1j
*
np.random.randn(npts))
else
:
x
=
noise
*
np.random.randn(npts)
for
p
in
periods:
if
rotary:
x
+
=
np.exp(
2j
*
np.pi
*
t
/
p)
else
:
x
+
=
np.cos(
2
*
np.pi
*
t
/
p)
return
t, x
def
complex_demodulation(t, x, central_period, hwidth
=
2
):
rotary
=
x.dtype.kind
=
=
'c'
c
=
np.exp(
-
1j
*
2
*
np.pi
*
t
/
central_period)
product
=
x
*
c
dt
=
t[
1
]
-
t[
0
]
hwpts
=
int
(
round
(hwidth
*
abs
(central_period)
/
dt))
nf
=
hwpts
*
2
+
1
x1
=
np.linspace(
-
1
,
1
, nf, endpoint
=
True
)
x1
=
x1[
1
:
-
1
]
w1
=
0.42
+
0.5
*
np.cos(x1
*
np.pi)
+
0.08
*
np.cos(x1
*
2
*
np.pi)
ytop
=
np.convolve(product, w1, mode
=
'same'
)
ybot
=
np.convolve(np.ones_like(product), w1, mode
=
'same'
)
demod
=
ytop
/
ybot
if
not
rotary:
demod
*
=
2
reconstructed
=
(demod
*
np.conj(c))
if
not
rotary:
reconstructed
=
reconstructed.real
if
np.sign(central_period) <
0
:
demod
=
np.conj(demod)
d
=
{
't'
:t,
'signal'
: x,
'hwpts'
: hwpts,
'demod'
: demod,
'reconstructed'
: reconstructed}
return
d
def
plot_demodulation(dm):
fig, axs
=
plt.subplots(
3
, sharex
=
True
)
resid
=
dm.get(
'signal'
)
-
dm.get(
'reconstructed'
)
if
dm.get(
'signal'
).dtype.kind
=
=
'c'
:
axs[
0
].plot(dm.get(
't'
), dm.get(
'signal'
).real, label
=
'signal.real'
)
axs[
0
].plot(dm.get(
't'
), dm.get(
'signal'
).imag, label
=
'signal.imag'
)
axs[
0
].plot(dm.get(
't'
), resid.real, label
=
'difference real'
)
axs[
0
].plot(dm.get(
't'
), resid.imag, label
=
'difference imag'
)
else
:
axs[
0
].plot(dm.get(
't'
), dm.get(
'signal'
), label
=
'signal'
)
axs[
0
].plot(dm.get(
't'
), dm.get(
'reconstructed'
), label
=
'reconstructed'
)
axs[
0
].plot(dm.get(
't'
), dm.get(
'signal'
)
-
dm.get(
'reconstructed'
), label
=
'difference'
)
axs[
0
].legend(loc
=
'upper right'
, fontsize
=
'small'
)
axs[
1
].plot(dm.get(
't'
), np.
abs
(dm.get(
'demod'
)), label
=
'amplitude'
, color
=
'C3'
)
axs[
1
].legend(loc
=
'upper right'
, fontsize
=
'small'
)
axs[
2
].plot(dm.get(
't'
), np.angle(dm.get(
'demod'
), deg
=
True
),
'.'
, label
=
'phase'
,
color
=
'C4'
)
axs[
2
].set_ylim(
-
180
,
180
)
axs[
2
].legend(loc
=
'upper right'
, fontsize
=
'small'
)
for
ax
in
axs:
ax.locator_params(axis
=
'y'
, nbins
=
5
)
return
fig, axs
def
test_demodulation(periods, central_period,
noise
=
0
,
rotary
=
False
,
hwidth
=
1
,
npts
=
1000
,
dt
=
1.0
/
24
):
t, x
=
gen_test_data(periods, noise
=
noise, rotary
=
rotary,
npts
=
npts, dt
=
dt)
dm
=
complex_demodulation(t, x, central_period, hwidth
=
hwidth)
fig, axs
=
plot_demodulation(dm)
return
fig, axs, dm
test_demodulation([
12.0
/
24
],
12.0
/
24
);
test_demodulation([
11.0
/
24
],
12.0
/
24
)