diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 0025978cc..13b2bc4f6 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -2,6 +2,18 @@
 
 <h3>New features since last release</h3>
 
+* Rectangular GKP states are now supported.
+  [(#668)](https://github.com/XanaduAI/strawberryfields/pull/668)
+
+  For example,
+  
+  ```python
+   prog = sf.Program(1)
+   
+   with prog.context as q:
+       sf.ops.GKP(shape="rectangular", alpha=2) | q[0]
+  ```
+
 <h3>Breaking Changes</h3>
 
 <h3>Bug fixes</h3>
@@ -12,6 +24,8 @@
 
 This release contains contributions from (in alphabetical order):
 
+J. Eli Bourassa
+
 # Release 0.21.0 (current release)
 
 <h3>New features since last release</h3>
diff --git a/strawberryfields/backends/bosonicbackend/backend.py b/strawberryfields/backends/bosonicbackend/backend.py
index c276c1797..f5a20983c 100644
--- a/strawberryfields/backends/bosonicbackend/backend.py
+++ b/strawberryfields/backends/bosonicbackend/backend.py
@@ -540,7 +540,9 @@ def prepare_cat_real_rep(self, a, theta, p, ampl_cutoff, D):
 
         return weights, means, cov
 
-    def prepare_gkp(self, state, epsilon, ampl_cutoff, representation="real", shape="square"):
+    def prepare_gkp(
+        self, state, epsilon, ampl_cutoff, representation="real", shape="square", alpha=1
+    ):
         r"""Prepares the arrays of weights, means and covs for a finite energy GKP state.
 
         GKP states are qubits, with the qubit state defined by:
@@ -555,6 +557,7 @@ def prepare_gkp(self, state, epsilon, ampl_cutoff, representation="real", shape=
             ampl_cutoff (float): this determines how many terms to keep
             representation (str): ``'real'`` or ``'complex'`` reprsentation
             shape (str): shape of the lattice; default 'square'
+            alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
 
         Returns:
             tuple: arrays of the weights, means and covariances for the state
@@ -566,16 +569,24 @@ def prepare_gkp(self, state, epsilon, ampl_cutoff, representation="real", shape=
         if representation == "complex":
             raise NotImplementedError("The complex description of GKP is not implemented")
 
-        if shape != "square":
-            raise NotImplementedError("Only square GKP are implemented for now")
+        if shape not in ["square", "rectangular"]:
+            raise NotImplementedError("Only square and rectangular GKP are implemented.")
+
+        if shape == "square":
+            if alpha != 1:
+                raise ValueError(
+                    "For square GKPs, alpha must be 1. For alpha not equal to "
+                    + "1, use shape='rectangular'."
+                )
 
         theta, phi = state[0], state[1]
 
-        def coeff(peak_loc):
+        def coeff(peak_loc, alpha):
             """Returns the value of the weight for a given peak.
 
             Args:
                 peak_loc (array): location of the ideal peak in phase space
+                alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
 
             Returns:
                 float: weight of the peak
@@ -610,7 +621,7 @@ def coeff(peak_loc):
             prefactor = np.exp(
                 -np.pi
                 * 0.25
-                * (l ** 2 + m ** 2)
+                * ((l * np.sqrt(alpha)) ** 2 + (m / np.sqrt(alpha)) ** 2)
                 * (1 - np.exp(-2 * epsilon))
                 / (1 + np.exp(-2 * epsilon))
             )
@@ -649,7 +660,7 @@ def coeff(peak_loc):
         )
 
         # Calculate the weights for each peak
-        weights = coeff(means)
+        weights = coeff(means, alpha)
         filt = abs(weights) > ampl_cutoff
         weights = weights[filt]
 
@@ -657,7 +668,11 @@ def coeff(peak_loc):
         # Apply finite energy effect to means
         means = means[filt]
 
+        means[:, 0] *= np.sqrt(alpha)
+        means[:, 1] /= np.sqrt(alpha)
+
         means *= 0.5 * damping * np.sqrt(np.pi * self.circuit.hbar)
+
         # Covariances all the same
         covs = (
             0.5
diff --git a/strawberryfields/backends/fockbackend/backend.py b/strawberryfields/backends/fockbackend/backend.py
index befa6bd84..71d64da72 100644
--- a/strawberryfields/backends/fockbackend/backend.py
+++ b/strawberryfields/backends/fockbackend/backend.py
@@ -295,7 +295,7 @@ def measure_fock(self, modes, shots=1, select=None, **kwargs):
         return self.circuit.measure_fock(self._remap_modes(modes), select=select)
 
     def prepare_gkp(
-        self, state, epsilon, ampl_cutoff, representation="real", shape="square", mode=None
+        self, state, epsilon, ampl_cutoff, representation="real", shape="square", alpha=1, mode=None
     ):
         r"""Prepares the Fock representation of a finite energy GKP state.
 
@@ -311,6 +311,7 @@ def prepare_gkp(
             amplcutoff (float): this determines how many terms to keep
             representation (str): ``'real'`` or ``'complex'`` reprsentation
             shape (str): shape of the lattice; default 'square'
+            alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
 
         Returns:
             tuple: arrays of the weights, means and covariances for the state
@@ -322,8 +323,15 @@ def prepare_gkp(
         if representation == "complex":
             raise NotImplementedError("The complex description of GKP is not implemented")
 
-        if shape != "square":
-            raise NotImplementedError("Only square GKP are implemented for now")
+        if shape not in ["square", "rectangular"]:
+            raise NotImplementedError("Only square and rectangular GKP are implemented.")
+
+        if shape == "square":
+            if alpha != 1:
+                raise ValueError(
+                    "For square GKPs, alpha must be 1. For alpha not equal to "
+                    + "1, use shape='rectangular'."
+                )
 
         theta, phi = state[0], state[1]
-        self.circuit.prepare_gkp(theta, phi, epsilon, ampl_cutoff, self._remap_modes(mode))
+        self.circuit.prepare_gkp(theta, phi, epsilon, ampl_cutoff, alpha, self._remap_modes(mode))
diff --git a/strawberryfields/backends/fockbackend/circuit.py b/strawberryfields/backends/fockbackend/circuit.py
index 41868fb2d..b5962b5ae 100644
--- a/strawberryfields/backends/fockbackend/circuit.py
+++ b/strawberryfields/backends/fockbackend/circuit.py
@@ -799,13 +799,15 @@ def measure_homodyne(self, phi, mode, select=None, **kwargs):
         # `homodyne_sample` will always be a single value since multiple shots is not supported
         return np.array([[homodyne_sample]])
 
-    def prepare_gkp(self, theta, phi, epsilon, ampl_cutoff, mode):
+    def prepare_gkp(self, theta, phi, epsilon, ampl_cutoff, alpha, mode):
         """
         Prepares a mode in a GKP state.
         """
 
         if self._pure:
-            self.prepare(ops.square_gkp_state(theta, phi, epsilon, ampl_cutoff, self._trunc), mode)
+            self.prepare(
+                ops.rect_gkp_state(theta, phi, epsilon, ampl_cutoff, alpha, self._trunc), mode
+            )
         else:
-            st = ops.square_gkp_state(theta, phi, epsilon, ampl_cutoff, self._trunc)
+            st = ops.rect_gkp_state(theta, phi, epsilon, ampl_cutoff, alpha, self._trunc)
             self.prepare(np.outer(st, st.conjugate()), mode)
diff --git a/strawberryfields/backends/fockbackend/ops.py b/strawberryfields/backends/fockbackend/ops.py
index 476792eef..9d004bb52 100644
--- a/strawberryfields/backends/fockbackend/ops.py
+++ b/strawberryfields/backends/fockbackend/ops.py
@@ -515,7 +515,7 @@ def hermiteVals(q_mag, num_bins, m_omega_over_hbar, trunc):
     return q_tensor, Hvals
 
 
-def gkp_displacements(t, k, epsilon):
+def gkp_displacements(t, k, epsilon, alpha):
     """
     Helper function to generate the displacements parameters associated with the teeth of
     GKP computational basis state k.
@@ -524,14 +524,15 @@ def gkp_displacements(t, k, epsilon):
         t (array): the teeth of GKP computational basis
         k (int): a computational basis state label, can be either 0 or 1
         epsilon (float): finite energy parameter of the state
+        alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
 
     Returns:
         array: the displacements
     """
-    return np.sqrt(0.5 * np.pi) * (2 * t + k) / np.cosh(epsilon)
+    return np.sqrt(0.5 * np.pi * alpha) * (2 * t + k) / np.cosh(epsilon)
 
 
-def gkp_coeffs(t, k, epsilon):
+def gkp_coeffs(t, k, epsilon, alpha):
     """
     Helper function to generate the coefficient parameters associated with the teeth of
     GKP computational basis state k.
@@ -540,22 +541,24 @@ def gkp_coeffs(t, k, epsilon):
         t (array): the teeth of GKP computational basis
         k (int): a computational basis state label, can be either 0 or 1
         epsilon (float): finite energy parameter of the state
+        alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
 
     Returns:
         array: the coefficients
     """
-    return np.exp(-0.5 * np.pi * np.tanh(epsilon) * (k + 2 * t) ** 2)
+    return np.exp(-0.5 * np.pi * alpha * np.tanh(epsilon) * (k + 2 * t) ** 2)
 
 
 @functools.lru_cache()
-def square_gkp_basis_state(i, epsilon, ampl_cutoff, cutoff):
+def rect_gkp_basis_state(i, epsilon, ampl_cutoff, alpha, cutoff):
     """
-    Generate the Fock expansion of a (subnormalized) computational GKP basis state. Normalization occurs in the ``square_gkp_state`` method.
+    Generate the Fock expansion of a (subnormalized) computational GKP basis state. Normalization occurs in the ``rect_gkp_state`` method.
 
     Args:
         i (int): a computational basis state label, can be either 0 or 1
         epsilon (float): finite energy parameter of the state
         ampl_cutoff (float): this determines how many terms to keep in the Hilbert space expansion
+        alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
         cutoff (int): Fock space truncation
 
     Returns
@@ -563,16 +566,16 @@ def square_gkp_basis_state(i, epsilon, ampl_cutoff, cutoff):
 
     """
     z_max = int(np.ceil(np.sqrt(-0.25 / np.pi * np.log(ampl_cutoff) / np.tanh(epsilon))))
-    coeffs = [gkp_coeffs(t, i, epsilon) for t in range(-z_max, z_max + 1)]
+    coeffs = [gkp_coeffs(t, i, epsilon, alpha) for t in range(-z_max, z_max + 1)]
     r = -0.5 * np.log(np.tanh(epsilon))
-    alphas = [gkp_displacements(t, i, epsilon) for t in range(-z_max, z_max + 1)]
-    num_kets = len(alphas)
-    ket = [coeffs[j] * displacedSqueezed(alphas[j], 0, r, 0, cutoff) for j in range(num_kets)]
+    disps = [gkp_displacements(t, i, epsilon, alpha) for t in range(-z_max, z_max + 1)]
+    num_kets = len(disps)
+    ket = [coeffs[j] * displacedSqueezed(disps[j], 0, r, 0, cutoff) for j in range(num_kets)]
     return sum(ket)
 
 
 @functools.lru_cache()
-def square_gkp_state(theta, phi, epsilon, ampl_cutoff, cutoff):
+def rect_gkp_state(theta, phi, epsilon, ampl_cutoff, alpha, cutoff):
     r"""
     Generate the Fock expansion of an abitrary GKP state parametrized as
     :math:`|\psi\rangle = \cos{\tfrac{\theta}{2}} \vert 0 \rangle_{\rm gkp} + e^{-i \phi} \sin{\tfrac{\theta}{2}} \vert 1 \rangle_{\rm gkp}`.
@@ -582,6 +585,7 @@ def square_gkp_state(theta, phi, epsilon, ampl_cutoff, cutoff):
         phi (float): the longitude with respect to the x-axis in the Bloch sphere
         epsilon (float): finite energy parameter of the state
         ampl_cutoff (float): this determines how many terms to keep
+        alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
         cutoff (int): Fock space truncation
 
     Returns:
@@ -589,8 +593,8 @@ def square_gkp_state(theta, phi, epsilon, ampl_cutoff, cutoff):
     """
     qubit_coeff0 = np.cos(theta / 2)
     qubit_coeff1 = np.sin(theta / 2) * np.exp(-1j * phi)
-    ket0 = square_gkp_basis_state(0, epsilon, ampl_cutoff, cutoff)
-    ket1 = square_gkp_basis_state(1, epsilon, ampl_cutoff, cutoff)
+    ket0 = rect_gkp_basis_state(0, epsilon, ampl_cutoff, alpha, cutoff)
+    ket1 = rect_gkp_basis_state(1, epsilon, ampl_cutoff, alpha, cutoff)
     ket = qubit_coeff0 * ket0 + qubit_coeff1 * ket1
     ket /= np.linalg.norm(ket)
     return ket
diff --git a/strawberryfields/ops.py b/strawberryfields/ops.py
index a8c563558..8d8389cab 100644
--- a/strawberryfields/ops.py
+++ b/strawberryfields/ops.py
@@ -949,14 +949,21 @@ class GKP(Preparation):
         cutoff (float): this determines how many terms to keep
         representation (str): ``'real'`` or ``'complex'`` reprsentation
         shape (str): shape of the lattice; default ``'square'``
+        alpha (float): peak spacing in q is given by sqrt(alpha * pi * hbar)
     """
 
     def __init__(
-        self, state=None, epsilon=0.2, ampl_cutoff=1e-12, representation="real", shape="square"
+        self,
+        state=None,
+        epsilon=0.2,
+        ampl_cutoff=1e-12,
+        representation="real",
+        shape="square",
+        alpha=1,
     ):
         if state is None:
             state = [0, 0]
-        super().__init__([state, epsilon, ampl_cutoff, representation, shape])
+        super().__init__([state, epsilon, ampl_cutoff, representation, shape, alpha])
 
     def _apply(self, reg, backend, **kwargs):
         backend.prepare_gkp(*self.p, mode=reg[0])
diff --git a/tests/backend/test_bosonic_backend.py b/tests/backend/test_bosonic_backend.py
index d70f95ba4..b46abc716 100644
--- a/tests/backend/test_bosonic_backend.py
+++ b/tests/backend/test_bosonic_backend.py
@@ -474,6 +474,35 @@ def test_gkp_complex(self):
         with pytest.raises(NotImplementedError):
             backend.run_prog(prog)
 
+    @pytest.mark.parametrize("eps", EPS_VALS)
+    def test_sensor_gkp(self, eps):
+        r"""Checks that the GKP sensor state is invariant under Fourier gate."""
+        x = np.linspace(-3 * np.sqrt(np.pi), 3 * np.sqrt(np.pi), 40)
+        p = np.copy(x)
+
+        # Prepare GKP 0 and apply Hadamard
+        prog_sensor = sf.Program(1)
+        with prog_sensor.context as q:
+            sf.ops.GKP(state=[np.pi / 2, 0], epsilon=eps, shape="rectangular", alpha=2) | q[0]
+
+        backend_sensor = bosonic.BosonicBackend()
+        backend_sensor.run_prog(prog_sensor)
+        state_sensor = backend_sensor.state()
+        wigner_sensor = state_sensor.wigner(0, x, p)
+
+        # Prepare GKP +
+        prog_sensor_F = sf.Program(1)
+        with prog_sensor_F.context as qF:
+            sf.ops.GKP(state=[np.pi / 2, 0], epsilon=eps, shape="rectangular", alpha=2) | qF[0]
+            sf.ops.Rgate(np.pi / 2) | qF[0]
+
+        backend_sensor_F = bosonic.BosonicBackend()
+        backend_sensor_F.run_prog(prog_sensor_F)
+        state_sensor_F = backend_sensor_F.state()
+        wigner_sensor_F = state_sensor_F.wigner(0, x, p)
+
+        assert np.allclose(wigner_sensor, wigner_sensor_F)
+
     @pytest.mark.parametrize("fock_eps", [0.1, 0.2, 0.3])
     def test_gkp_wigner_compare_backends(self, fock_eps):
         """Test correct Wigner function are generated by comparing fock and bosonic outputs.
@@ -501,6 +530,34 @@ def test_gkp_wigner_compare_backends(self, fock_eps):
         W_bosonic = gkp_bosonic.wigner(0, xvec, xvec)
         assert np.allclose(W_fock, W_bosonic)
 
+    @pytest.mark.parametrize("fock_eps", [0.1, 0.2, 0.3])
+    @pytest.mark.parametrize("rect_ratio", [0.8, 1, 1.5])
+    def test_rect_gkp_wigner_compare_backends(self, fock_eps, rect_ratio):
+        """Test that correct Wigner function is generated by comparing fock and bosonic outputs.
+        Since the fock backend cannot simulate very small epsilon, we restrict to fock_eps > 0.1."""
+        theta = np.pi / 8
+        phi = np.pi / 6
+        ket = [theta, phi]
+        xvec = np.linspace(-1, 1, 200)
+        prog_fock = sf.Program(1)
+        with prog_fock.context as q:
+            sf.ops.GKP(epsilon=fock_eps, state=ket, shape="rectangular", alpha=rect_ratio) | q
+        cutoff = 100
+        hbar = 2
+        eng_fock = sf.Engine("fock", backend_options={"cutoff_dim": cutoff, "hbar": hbar})
+        gkp_fock = eng_fock.run(prog_fock).state
+        W_fock = gkp_fock.wigner(0, xvec, xvec)
+
+        prog_bosonic = sf.Program(1)
+
+        with prog_bosonic.context as q:
+            sf.ops.GKP(epsilon=fock_eps, state=ket, shape="rectangular", alpha=rect_ratio) | q
+        hbar = 2
+        eng_bosonic = sf.Engine("bosonic", backend_options={"hbar": hbar})
+        gkp_bosonic = eng_bosonic.run(prog_bosonic).state
+        W_bosonic = gkp_bosonic.wigner(0, xvec, xvec)
+        assert np.allclose(W_fock, W_bosonic)
+
 
 class TestBosonicUserSpecifiedState:
     """Checks the Bosonic preparation method."""
diff --git a/tests/frontend/test_ops_gate.py b/tests/frontend/test_ops_gate.py
index 422e8fd9b..99a4872dd 100644
--- a/tests/frontend/test_ops_gate.py
+++ b/tests/frontend/test_ops_gate.py
@@ -174,12 +174,22 @@ class TestGKPBasics:
     @pytest.mark.parametrize("ampl", [0, 0.001])
     @pytest.mark.parametrize("eps", [0, 0.001])
     @pytest.mark.parametrize("r", ["real", "complex"])
-    @pytest.mark.parametrize("s", ["square"])
-    def test_gkp_str_representation(self, state, ampl, eps, r, s):
+    @pytest.mark.parametrize("s", ["square", "rectangular"])
+    @pytest.mark.parametrize("alph", [0.8, 1.1])
+    def test_gkp_str_representation(self, state, ampl, eps, r, s, alph):
         """Test the string representation of the GKP operation"""
         assert (
-            str(ops.GKP(state=state, ampl_cutoff=ampl, epsilon=eps, representation=r, shape=s))
-            == f"GKP({str(state)}, {str(eps)}, {str(ampl)}, {r}, {s})"
+            str(
+                ops.GKP(
+                    state=state,
+                    ampl_cutoff=ampl,
+                    epsilon=eps,
+                    representation=r,
+                    shape=s,
+                    alpha=alph,
+                )
+            )
+            == f"GKP({str(state)}, {str(eps)}, {str(ampl)}, {r}, {s}, {str(alph)})"
         )