Emulator readout module modification

readout_example/readout_example.py

@@ -0,0 +1,88 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+from qibolab.instruments.emulator.readout import ReadoutSimulator
+from qibolab.pulses import ReadoutPulse
+from qibolab.qubits import Qubit
+
+bare_resonator_frequency = 5.045e9
+nshots = 100
+
+SNR = 50  # dB
+READOUT_AMPLITUDE = 1
+NOISE_AMP = np.power(10, -SNR / 20)
+AWGN = lambda t: np.random.normal(loc=0, scale=NOISE_AMP, size=len(t))
+
+qb = Qubit(
+    0,
+    bare_resonator_frequency=bare_resonator_frequency,
+    drive_frequency=3.99e9,
+    anharmonicity=-263e6,
+)
+readout = ReadoutSimulator(
+    qubit=qb,
+    g=10e6,
+    noise_model=AWGN,
+    internal_Q=2.5e6,
+    coupling_Q=6e4,
+    sampling_rate=1966.08e6,
+)
+
+
+# demonstrates effect of state dependent dispersive shift on amplitude of reflected microwave from resonator;
+span = 1e6
+center_frequency = bare_resonator_frequency
+freq_sweep = np.linspace(center_frequency - span / 2, center_frequency + span / 2, 1000)
+y_gnd = np.abs(readout.ground_s21(freq_sweep))
+y_exc = np.abs(readout.excited_s21(freq_sweep))
+
+freq_sweep /= 1e9
+plt.plot(freq_sweep, y_gnd, label=r"$|0\rangle$")
+plt.plot(freq_sweep, y_exc, label=r"$|1\rangle$")
+plt.ylabel("|S21| [arb. units]")
+plt.xlabel("Readout Frequency [GHz]")
+plt.legend()
+plt.grid()
+plt.tight_layout()
+# plt.savefig("S21.png", dpi=300)
+plt.show()
+
+freq_sweep *= 1e9
+
+# demonstrates effect of state dependent dispersive shift on phase of reflected microwave from resonator; (for domain of [-pi/2,pi/2])
+y_gnd1 = np.angle(readout.ground_s21(freq_sweep))
+y_exc1 = np.angle(readout.excited_s21(freq_sweep))
+freq_sweep /= 1e9
+plt.plot(freq_sweep, y_gnd1, label=r"$|0\rangle$")
+plt.plot(freq_sweep, y_exc1, label=r"$|1\rangle$")
+plt.ylabel(r"$\theta$ phase [rad]")
+plt.xlabel("Readout Frequency [GHz]")
+plt.legend()
+plt.grid()
+plt.tight_layout()
+# plt.savefig("phase.png", dpi=300)
+plt.show()
+
+# demonstrates the separation of V_I/V_Q data of reflected microwave on V_I/V_Q plane
+ro_frequency = 5.0450e9 + readout.lambshift
+ro_pulse = ReadoutPulse(
+    start=0,
+    duration=1000,
+    amplitude=READOUT_AMPLITUDE,
+    frequency=ro_frequency,
+    shape="Rectangular()",
+    relative_phase=0,
+)
+
+rgnd = [readout.simulate_ground_state_iq(ro_pulse) for k in range(nshots)]
+rexc = [readout.simulate_excited_state_iq(ro_pulse) for k in range(nshots)]
+plt.scatter(np.real(rgnd), np.imag(rgnd), label=r"$|0\rangle$")
+plt.scatter(np.real(rexc), np.imag(rexc), label=r"$|1\rangle$")
+# plt.xlim([0,1])
+# plt.ylim([-1,1])
+plt.xlabel(r"$V_I$")
+plt.ylabel(r"$V_Q$")
+plt.legend()
+plt.tight_layout()
+# plt.savefig("IQ.png", dpi=300)
+plt.show()

src/qibolab/instruments/emulator/readout.py

@@ -21,8 +21,8 @@ def lamb_shift(g, delta):
 
 
 def dispersive_shift(g, delta, alpha):
-    """Calculates the dispersive shift of the readout resonator for the first excited state of the qubit
-    @see https://arxiv.org/abs/2106.06173, equation 35
+    """Calculates the dispersive shift of the readout resonator for depending on state of the qubit
+    @see https://arxiv.org/pdf/1904.06560, equation 146, we omit the negative sign as we included it in our definition of delta
 
     Args:
         g (float): Coupling strength between readout resonator and qubit in Hz.
@@ -32,7 +32,7 @@ def dispersive_shift(g, delta, alpha):
     Returns:
         chi (float): Dispersive shift in Hz.
     """
-    return 2 * alpha * g * g / delta / delta
+    return g * g / delta * (1 / (1 + (delta / alpha)))
 
 
 def s21_function(resonator_frequency, total_Q, coupling_Q):
@@ -65,12 +65,18 @@ def __init__(
             coupling_Q (float): Coupling/external Q factor of the readout resonator.
             sampling_rate (float): Sampling rate of the ADC/digitizer.
         """
-
-        delta = qubit.bare_resonator_frequency - qubit.drive_frequency
-        ground_state_frequency = qubit.bare_resonator_frequency + lamb_shift(g, delta)
-        excited_state_frequency = ground_state_frequency + dispersive_shift(
+        # maintaining the definition of |0> = |e> = (1 0) with the JC Hamiltonian model
+        delta = qubit.drive_frequency - qubit.bare_resonator_frequency
+        dressed_frequency = qubit.bare_resonator_frequency - lamb_shift(g, delta)
+        ground_state_frequency = dressed_frequency - dispersive_shift(
+            g, delta, qubit.anharmonicity
+        )
+        excited_state_frequency = dressed_frequency + dispersive_shift(
             g, delta, qubit.anharmonicity
         )
+        # print(f"ground_state_frequency= {ground_state_frequency}; excited_state_frequency= {excited_state_frequency}")
+
+        self.lambshift = -lamb_shift(g, delta)
         self.noise_model = noise_model
         self.sampling_rate = sampling_rate
 
@@ -99,8 +105,8 @@ def simulate_excited_state_iq(self, pulse: ReadoutPulse):
         return self.simulate_and_demodulate(s21, pulse)
 
     def simulate_and_demodulate(self, s21: complex, pulse: ReadoutPulse):
-        """Simulates the readout pulse for a given S21-parameter and
-        demodulates it.
+        """Simulates the readout pulse for a given S21-parameter and homodyne
+        demodulation/2nd stage of heterodyne demodulation.
 
         Args:
             s21 (complex): Complex S21 parameter.
@@ -109,34 +115,28 @@ def simulate_and_demodulate(self, s21: complex, pulse: ReadoutPulse):
         Returns:
             IQ (complex): IQ data for a shot.
         """
-        logmag = np.abs(s21)
-        phase = np.angle(s21)
+        reflected_amplitude = np.abs(s21)
+        reflected_phase = np.angle(s21)
 
-        env_I, env_Q = pulse.envelope_waveforms(self.sampling_rate / 1e9)
+        env_I, env_Q = pulse.envelope_waveforms(
+            self.sampling_rate / 1e9
+        )  # Gigasample per second
 
-        start = int(pulse.start * 1e-9 * self.sampling_rate)
-        t = np.arange(start, start + len(env_I)) / self.sampling_rate
-        reference_signal_I = np.sin(
-            2 * np.pi * pulse.frequency * t + pulse.relative_phase
-        )
-        reference_signal_Q = np.cos(
-            2 * np.pi * pulse.frequency * t + pulse.relative_phase
-        )
+        start = int(pulse.start * 1e-9 * self.sampling_rate)  # n = gigasample index
+        t = np.arange(start, start + len(env_I)) / self.sampling_rate  # t_n
 
-        waveform = (
-            logmag
-            * pulse.amplitude
-            * (
-                env_I.data
-                * np.cos(pulse.relative_phase + phase)
-                * np.sin(2 * np.pi * pulse.frequency * t)
-                + env_Q.data
-                * np.sin(pulse.relative_phase + phase)
-                * np.cos(2 * np.pi * pulse.frequency * t)
-                + self.noise_model(t)
-            )
-        )
+        # Low-pass filtered I-component (with intermediate frequency = carrier frequency)
+        I_filtered = reflected_amplitude * np.cos(
+            2 * np.pi * t * pulse.frequency + pulse.relative_phase + reflected_phase
+        ) + self.noise_model(t)
 
-        return np.dot(waveform, reference_signal_I) + 1j * np.dot(
-            waveform, reference_signal_Q
-        )
+        # Low-pass filtered Q-component
+        Q_filtered = reflected_amplitude * np.sin(
+            2 * np.pi * t * pulse.frequency + pulse.relative_phase + reflected_phase
+        ) + self.noise_model(t)
+
+        z = I_filtered + 1j * Q_filtered
+        z *= np.exp(-1j * 2 * np.pi * t * pulse.frequency)
+        z = np.sum(z) / len(t)
+
+        return z