Note
Click here to download the full example code
Compute a spectrogramΒΆ
Compute the spectrogram of a set of SFTs. This is useful to produce visualizations of the Doppler modulation of a CW signal.
9 import os
10
11 import matplotlib.pyplot as plt
12
13 import pyfstat
14
15 # not github-action compatible
16 # plt.rcParams["font.family"] = "serif"
17 # plt.rcParams["font.size"] = 18
18 # plt.rcParams["text.usetex"] = True
19
20 # workaround deprecation warning
21 # see https://github.com/matplotlib/matplotlib/issues/21723
22 plt.rcParams["axes.grid"] = False
23
24 label = "PyFstatExampleSpectrogram"
25 outdir = os.path.join("PyFstat_example_data", label)
26 logger = pyfstat.set_up_logger(label=label, outdir=outdir)
27
28 depth = 5
29
30 data_parameters = {
31 "sqrtSX": 1e-23,
32 "tstart": 1000000000,
33 "duration": 2 * 365 * 86400,
34 "detectors": "H1",
35 "Tsft": 1800,
36 }
37
38 signal_parameters = {
39 "F0": 100.0,
40 "F1": 0,
41 "F2": 0,
42 "Alpha": 0.0,
43 "Delta": 0.5,
44 "tp": data_parameters["tstart"],
45 "asini": 25.0,
46 "period": 50 * 86400,
47 "tref": data_parameters["tstart"],
48 "h0": data_parameters["sqrtSX"] / depth,
49 "cosi": 1.0,
50 }
51
52 # making data
53 data = pyfstat.BinaryModulatedWriter(
54 label=label, outdir=outdir, **data_parameters, **signal_parameters
55 )
56 data.make_data()
57
58 logger.info("Loading SFT data and computing normalized power...")
59 freqs, times, sft_data = pyfstat.utils.get_sft_as_arrays(data.sftfilepath)
60 sft_power = sft_data["H1"].real ** 2 + sft_data["H1"].imag ** 2
61 normalized_power = (
62 2 * sft_power / (data_parameters["Tsft"] * data_parameters["sqrtSX"] ** 2)
63 )
64
65 plotfile = os.path.join(outdir, label + ".png")
66 logger.info(f"Plotting to file: {plotfile}")
67 fig, ax = plt.subplots(figsize=(0.8 * 16, 0.8 * 9))
68 ax.set(xlabel="Time [days]", ylabel="Frequency [Hz]", ylim=(99.98, 100.02))
69 c = ax.pcolormesh(
70 (times["H1"] - times["H1"][0]) / 86400,
71 freqs,
72 normalized_power,
73 cmap="inferno_r",
74 shading="nearest",
75 )
76 fig.colorbar(c, label="Normalized Power")
77 plt.tight_layout()
78 fig.savefig(plotfile)
Total running time of the script: ( 0 minutes 0.000 seconds)