diff --git a/src/phasespace.jl b/src/phasespace.jl
index fc0e1a26..1f7e9d67 100644
--- a/src/phasespace.jl
+++ b/src/phasespace.jl
@@ -558,3 +558,34 @@ function _calc_ylm_norm(l, m)
     end
     return n
 end
+
+"""
+    squeeze(b::FockBasis, x)
+
+Squeeze operator on a `FockBasis` described by  
+``S = e^{\\frac{1}{2}(x^* a^2 - x (a^\dagger)^2)}``.
+For ``x = e^{i θ}`` the squeezing is along `θ/2`.
+"""
+function squeeze(b::FockBasis,x)
+    a = destroy(b)
+    ad = create(b)
+    s = exp(conj(x)*dense(a)^2/2 - x*dense(ad)^2/2)
+    return s
+end
+
+"""
+    squeeze(b::SpinBasis, x)
+
+Squeeze operator on a `SpinBasis` described by  
+``S = e^{\\frac{1}{2 N}(x^* J_-^2 - x J_+^2)}``.
+For ``x = e^{i θ}`` the squeezing is along `θ/2`.
+Due to the finiteness of the Hilbert space a too large
+ squeezing ``(|x| > \sqrt{N})`` will create an over-squeezed state. 
+"""
+function squeeze(b::SpinBasis,x)
+    N = Int(b.spinnumber*2)
+    Jm = sigmam(b)/2
+    Jp = sigmap(b)/2
+    s = exp(conj(x)*dense(Jm)^2/2/N - x*dense(Jp)^2/2/N)
+    return s
+end
diff --git a/test/test_phasespace.jl b/test/test_phasespace.jl
index 6ddb21be..f7965211 100644
--- a/test/test_phasespace.jl
+++ b/test/test_phasespace.jl
@@ -207,4 +207,36 @@ end
 t2 = abs(int*2*(pi/res)^2)
 @test isapprox(t2, 0, atol=1e-2)
 
-end # testset
+# squeeze operator tests HO
+co = 100
+b_fock = FockBasis(co)
+a = destroy(b_fock)
+ad = create(b_fock);
+xx = rand() 
+s = squeeze(b_fock,xx*exp(1im*pi));
+
+@test log(real(variance((ad+a),s*fockstate(b_fock,0))))/2 ≈ xx
+
+s = squeeze(b_fock,xx*exp(1im*pi*rand()));
+
+@test log(real(variance((ad+a),s*fockstate(b_fock,0))))/2 < xx
+@test log(real(variance((ad+a),s*fockstate(b_fock,0))))/2 > -xx
+
+# squeeze oerator test SPIN
+N = 500
+b_spin = SpinBasis(N//2)
+ss = squeeze(b_spin,rand()*sqrt(N)*exp(1im*rand()*pi));
+st = ss*spindown(b_spin);
+
+# Heisenberg uncertainty test
+@test abs(variance(sigmax(b_spin)/2,st)*variance(sigmay(b_spin)/2,st)) ≥ abs2(expect(sigmaz(b_spin)/2,st))/4
+
+# small squeezing test
+xx = rand()/10
+ss = squeeze(b_spin,xx*sqrt(N));
+st = ss*spindown(b_spin);
+
+@test isapprox(2*log(real(variance(sigmax(b_spin)/2,st))/N*4) , -xx*sqrt(N), atol=1e-2)
+@test isapprox(2*log(real(variance(sigmay(b_spin)/2,st))/N*4) , xx*sqrt(N), atol=1e-2)
+
+end # testset
\ No newline at end of file