diff --git a/Project.toml b/Project.toml
index 21e3466..77871ab 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,7 +4,7 @@ keywords = ["AeroAcoustics"]
 license = "MIT"
 desc = "A package for post-processing of microphone array measurements"
 authors = ["Oliver Lylloff <ollyl@dtu.dk>"]
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
@@ -15,6 +15,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
 
 [compat]
 DSP = "0.6.7"
@@ -23,6 +24,7 @@ FFTW = "1.1"
 LazyArrays = "0.16"
 NLsolve = "4.4"
 Parameters = "0.12.3"
+ThreadsX = "0.1.10"
 julia = "1.6"
 
 [extras]
diff --git a/docs/src/index.md b/docs/src/index.md
index ffe802e..91238bd 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -50,7 +50,8 @@ E = Environment(
     Nx = 21,
     Ny = 21,
     xlim=(-0.5,0.5),
-    ylim=(-0.5,0.5)
+    ylim=(-0.5,0.5),
+    multi_thread = true # multi-threading can be enabled globally like this
     )
 ```
 Where the measurement distance `z0`, the microphone geometry `micgeom`, and the csm `CSM` are required variables. 
diff --git a/examples/Quick_start.ipynb b/examples/Quick_start.ipynb
index 327874c..9c252d6 100644
--- a/examples/Quick_start.ipynb
+++ b/examples/Quick_start.ipynb
@@ -2,18 +2,9 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "┌ Warning: PyPlot is using tkagg backend, which is known to cause crashes on MacOS (#410); use the MPLBACKEND environment variable to request a different backend.\n",
-      "└ @ PyPlot /Users/oliver/.julia/packages/PyPlot/XaELc/src/init.jl:198\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "using HDF5, AeroAcoustics, PyPlot"
    ]
@@ -27,7 +18,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
@@ -36,7 +27,7 @@
        "1.6"
       ]
      },
-     "execution_count": 2,
+     "execution_count": 4,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -53,7 +44,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -62,7 +53,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 21,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -74,13 +65,14 @@
     "    Nx = 21,\n",
     "    Ny = 21,\n",
     "    xlim=(-0.5,0.5),\n",
-    "    ylim=(-0.5,0.5)\n",
+    "    ylim=(-0.5,0.5),\n",
+    "    multi_thread = true, # multi-threading can be enabled globally like this\n",
     "    );"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -89,7 +81,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -98,7 +90,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [
     {
@@ -125,7 +117,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
@@ -134,7 +126,7 @@
        "(0.1, -0.05)"
       ]
      },
-     "execution_count": 8,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -153,14 +145,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "  4.742135 seconds (6.09 M allocations: 1.055 GiB, 5.72% gc time, 56.05% compilation time)\n"
+      "  4.762169 seconds (6.85 M allocations: 1.066 GiB, 6.01% gc time, 58.86% compilation time)\n"
      ]
     }
    ],
@@ -170,7 +162,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [
     {
@@ -204,14 +196,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "  2.603157 seconds (1.39 M allocations: 1.265 GiB, 11.42% gc time, 24.62% compilation time)\n"
+      "  2.377261 seconds (1.35 M allocations: 1.255 GiB, 12.46% gc time, 24.53% compilation time)\n"
      ]
     }
    ],
@@ -221,7 +213,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [
     {
@@ -255,18 +247,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "breaking at iter = 42...\n",
-      "breaking at iter = 47...\n",
-      "breaking at iter = 71...\n",
-      "breaking at iter = 55...\n",
-      "  2.753007 seconds (5.40 M allocations: 345.352 MiB, 3.79% gc time)\n"
+      "  2.979043 seconds (5.47 M allocations: 328.361 MiB, 3.92% gc time, 98.86% compilation time)\n"
      ]
     }
    ],
@@ -277,7 +265,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [
     {
@@ -311,7 +299,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 17,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -325,7 +313,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [
     {
@@ -339,7 +327,7 @@
        " 60.905           54.7118       54.8137    54.5308"
       ]
      },
-     "execution_count": 21,
+     "execution_count": 18,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -357,7 +345,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -366,7 +354,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 20,
    "metadata": {},
    "outputs": [
     {
@@ -380,7 +368,7 @@
        " 60.905           54.1891  54.7118       54.8137    54.5308"
       ]
      },
-     "execution_count": 23,
+     "execution_count": 20,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -395,26 +383,19 @@
    "source": [
     "The four methods `SPI`, `Clean-SC`, `DAMAS`, and `FISTA` now return very similar results of the source power for the 4 frequency bins defined in the environment (3000-4000 Hz)."
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Julia 1.6.2",
+   "display_name": "Julia 1.7.2",
    "language": "julia",
-   "name": "julia-1.6"
+   "name": "julia-1.7"
   },
   "language_info": {
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.6.2"
+   "version": "1.7.2"
   }
  },
  "nbformat": 4,
diff --git a/src/AeroAcoustics.jl b/src/AeroAcoustics.jl
index a5f4fe9..3a931b3 100644
--- a/src/AeroAcoustics.jl
+++ b/src/AeroAcoustics.jl
@@ -7,6 +7,7 @@ using Parameters
 using LazyArrays
 using FFTW
 import DSP
+using ThreadsX
 
 import Base.length,
        Base.push!,
@@ -104,4 +105,9 @@ function DR(CSM::FreqArray)
     return FreqArray(CSMd,CSM.fc)
 end
 
-end
\ No newline at end of file
+function check_multithread(multi_thread)
+    _foreach = multi_thread ? ThreadsX.foreach : Base.foreach
+    return _foreach
+end
+
+end # Module
\ No newline at end of file
diff --git a/src/beamforming.jl b/src/beamforming.jl
index 6e485ac..6c70144 100644
--- a/src/beamforming.jl
+++ b/src/beamforming.jl
@@ -25,10 +25,11 @@ end
 Calculate frequency-domain beamforming using cross-spectral matrix `csm` and steering vector `v`
 stored in Environment struct
 """
-function beamforming(E::Environment)
+function beamforming(E::Environment;multi_thread=E.multi_thread)
+    _foreach = AeroAcoustics.check_multithread(multi_thread)
     @unpack CSM_s,w,steeringvec,M,N,Nf,fn = E
     b = Array{Float64}(undef, N, Nf)
-    @views @inbounds for j in 1:Nf
+    @views @inbounds _foreach(1:Nf) do j
         st = E.steeringvec.arr[:,:,j]
         @. st = complex(w[:,j].*real(st), imag(st))
         bf_col!(b[:,j],st,selectdim(CSM_s.arr, 3, j))
diff --git a/src/cleansc.jl b/src/cleansc.jl
index 2883e02..7fd7d1b 100644
--- a/src/cleansc.jl
+++ b/src/cleansc.jl
@@ -8,7 +8,8 @@ CLEAN-SC algorithm for source identification and quantification optionally setti
 
 With inspiration from https://github.com/acoular/acoular/blob/66cba3cffb3bc72602c869f99347be76798f4ac1/acoular/fbeamform.py#L1496
 """
-function cleanSC(E;maxiter=50,ϕ=0.5,stopn=10,peak_removal=false,trust_region=nothing)
+function cleanSC(E;maxiter=50,ϕ=0.5,stopn=10,peak_removal=false,trust_region=nothing,multi_thread=E.multi_thread)
+    _foreach = AeroAcoustics.check_multithread(multi_thread)
     @unpack CSM_s,steeringvec,M,N,Nf,fn = E
     x = zeros(N,Nf)
     
@@ -20,7 +21,7 @@ function cleanSC(E;maxiter=50,ϕ=0.5,stopn=10,peak_removal=false,trust_region=no
         trust_indices = LinearIndices((E.Nx,E.Ny))
     end
     
-    @views @inbounds for j = 1:Nf
+    @views @inbounds _foreach(1:Nf) do j
         _cleanSC!(x[:,j],steeringvec.arr[:,:,j],CSM_s.arr[:,:,j],maxiter,ϕ,stopn,peak_removal,trust_indices)
     end
     return FreqArray(x,E.fn)
diff --git a/src/fista.jl b/src/fista.jl
index 78f5335..7c0803f 100644
--- a/src/fista.jl
+++ b/src/fista.jl
@@ -50,7 +50,7 @@ function _fista!(x,b,psf,gam;maxit=1000,tol=1e-6)
         x .= max.(real.(y),0.0)
         # stopping criterion
         if norm(x_ex-x)/gam <= tol*(1+norm(x))
-            println("breaking at iter = $(it)...")
+            #println("breaking at iter = $(it)...")
             break
         end
     end
diff --git a/src/psf.jl b/src/psf.jl
index 77449b7..12e4552 100644
--- a/src/psf.jl
+++ b/src/psf.jl
@@ -4,10 +4,11 @@
 Calculate frequency-domain point spread function using the Environment struct to access
 steeringvectors. Optionally, supply the index where the psf is centered, default is (N/2)+1.
 """
-function psf(E::Environment,cent::Int64=floor(Int,E.N/2)+1)
+function psf(E::Environment,cent::Int64=floor(Int,E.N/2)+1;multi_thread=false)
+    _foreach = AeroAcoustics.check_multithread(multi_thread)
     @unpack steeringvec,M,N,Nf,fn = E
     p = Array{Float64, 2}(undef, N, Nf)
-    @views @inbounds for j in 1:Nf
+    @views @inbounds _foreach(1:Nf) do j
         AeroAcoustics.psf_col!(p[:,j],steeringvec.arr[:,:,j],cent)
     end
     return FreqArray(p,fn)
diff --git a/src/steeringvectors.jl b/src/steeringvectors.jl
index 94f2bbc..5903d4e 100644
--- a/src/steeringvectors.jl
+++ b/src/steeringvectors.jl
@@ -81,17 +81,19 @@ end
 
 Pre-compute steeringvectors for beamforming using an `Environment` with the needed parameters.
 """
-function steeringvectors(E::Environment)
+function steeringvectors(E::Environment;multi_thread=E.multi_thread)
+    _foreach = AeroAcoustics.check_multithread(multi_thread)
     @unpack fn,c,M,N,Nf,D,D0 = E
     kw = 2pi*fn/c
     vi = Array{ComplexF64,3}(undef,N,M,Nf)
-    for j in 1:Nf
+    @views @inbounds _foreach(1:Nf) do j
         vi[:,:,j] .= 1 ./(D.*D0.*sum(1 ./D.^2,dims=2)).*exp.(-im.*kw[j].*(D.-D0))
     end
     steeringvec = FreqArray(permutedims(vi,(2,1,3)),fn)
 end
 
-function steeringvectors!(E::Environment)
+function steeringvectors!(E::Environment;multi_thread=E.multi_thread)
+    _foreach = AeroAcoustics.check_multithread(multi_thread)
     @unpack fn,c,M,N,Nf,D,D0,Ma,h = E
     vi = Array{ComplexF64,3}(undef,N,M,Nf)
 
@@ -99,7 +101,7 @@ function steeringvectors!(E::Environment)
         w = 2pi*fn
         dx,pc_pm = refraction_correction(E,h)
         ta = dx./c
-        for j in 1:Nf
+        @views @inbounds _foreach(1:Nf) do j
             vi[:,:,j] .= 1 ./(ta.*c.*D0.*sum(1 ./(ta.*c).^2,dims=2)).*exp.(-im.*w[j].*ta)
         end
         if E.ampcorr
@@ -109,7 +111,7 @@ function steeringvectors!(E::Environment)
         end
     else
         kw = 2pi*fn/c
-        for j in 1:Nf
+        @views @inbounds _foreach(1:Nf) do j
             vi[:,:,j] .= 1 ./(D.*D0.*sum(1 ./D.^2,dims=2)).*exp.(-im.*kw[j].*(D.-D0))
         end
     end
diff --git a/src/types.jl b/src/types.jl
index f098b98..21de779 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -43,13 +43,15 @@ setup, constants, and stores the relevant data together. The microphone array is
 - `c::Real=343.`: Speed of sound [m/s].
 - `Ma::Real=0.0`: Mach number (sign determines flow direction)
 - `h::Real=0.0`: Distance from array center to shear layer (Amiet correction) should be supplied when `Ma != 0`.
-
+- `multi_thread::Bool = false`: Enable multi-threading via the ThreadsX package.
 """
 @with_kw mutable struct Environment <: AeroAcousticType
+    # Required:
     micgeom::Matrix{<:AbstractFloat}
     z0::Real
     CSM::FreqArray
     flim::NTuple{2,Real} = extrema(CSM.fc)
+    # Optional:
     Nx::Int = 21
     Ny::Int = 21
     xlim::NTuple{2,Real} = (-1.,1.)
@@ -59,6 +61,10 @@ setup, constants, and stores the relevant data together. The microphone array is
     c::Real = 343.0 # Speed of sound
     Ma::Real = 0.0 # Flow Mach speed (in positive x-direction) TODO: Generalize to NTuple{3,Real}
     h::Real = 0.0 # Distance from array center to shear layer
+    multi_thread = false
+    #if multi_thread #&& Threads.nthreads() == 1
+    #    @warn "Multi-threading is enabled but julia was started with only 1 thread. This could lead to errornous results. Restart julia with `julia -t auto` to enable multi-threading."
+    #end
     ### Compute extra parameters
     dr::Bool = false
     Cinds = (CSM.fc.>=flim[1]) .& (CSM.fc.<=flim[2])