diff --git a/Project.toml b/Project.toml
index 54248df..d5e2d64 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,13 +4,12 @@ keywords = ["AeroAcoustics"]
 license = "MIT"
 desc = "A package for post-processing of microphone array measurements"
 authors = ["Oliver Lylloff <ollyl@dtu.dk>"]
-version = "0.2.4"
+version = "0.2.5"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
@@ -21,9 +20,9 @@ ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
 DSP = "0.6.7, 0.7"
 Distances = "0.10"
 FFTW = "1.1"
-LazyArrays = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 1"
 NLsolve = "4.4"
 Parameters = "0.12.3"
+Statistics = "1"
 ThreadsX = "0.1.10"
 julia = "1.6"
 
diff --git a/src/AeroAcoustics.jl b/src/AeroAcoustics.jl
index 3a931b3..a9b1afe 100644
--- a/src/AeroAcoustics.jl
+++ b/src/AeroAcoustics.jl
@@ -4,7 +4,6 @@ using LinearAlgebra
 using Statistics
 using NLsolve
 using Parameters
-using LazyArrays
 using FFTW
 import DSP
 using ThreadsX
diff --git a/src/csm.jl b/src/csm.jl
index 4f13f72..3ba8948 100644
--- a/src/csm.jl
+++ b/src/csm.jl
@@ -25,14 +25,10 @@ end
 
 # flat_t is a helper function to reshape a S x M matrix to a vector of vectors
 function flat_t(t::AbstractArray{T,N}) where {T <: AbstractFloat, N}
-    if size(t,1) < size(t,2)
+    if size(t, 1) < size(t, 2)
         t = permutedims(t)
     end
-    ta = Array{Vector{T}}(undef,size(t,2))
-    for i in axes(t,2)
-        ta[i] = view(t,:,i)
-    end
-    return ta
+    return map(i -> view(t, :, i), axes(t, 2))
 end
 
 """
@@ -45,7 +41,8 @@ function csm(t::AbstractArray{T};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n
     csm(flat_t(t);n=n,noverlap=noverlap,fs=fs,win=win,scaling=scaling)
 end
 
-function csm(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat
+function csm(t::Vector{<:AbstractVector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum",multi_thread=true) where T <: AbstractFloat
+    _foreach = AeroAcoustics.check_multithread(multi_thread)
     M = length(t)
     Nf = div(n,2)+1
     nout = div((length(t[1]) - n), n - noverlap)+1
@@ -53,9 +50,12 @@ function csm(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(
     C = Array{Complex{T}}(undef,M,M,Nf)
     S = DSP.stft.(t, n, noverlap; fs=fs, window=win, onesided=true)
     Sc = conj.(S)
-    for m in 1:M
-        for j in m:M
-            C[j,m,:] .= dropdims(mean(LazyArray(@~ Sc[m].*S[j]),dims=2);dims=2)
+    @views @inbounds _foreach(1:Nf) do ω
+        for m in 1:M
+            for j in m:M
+                C[j, m, ω] = mean(Sc[m][ω,:] .* S[j][ω,:])
+                C[m, j, ω] = conj(C[j, m, ω])  
+            end
         end
     end
 
@@ -68,9 +68,6 @@ function csm(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(
     C .*= scale
     C[:,:,2:end-1] .*= 2
 
-    for ω in 1:Nf
-        C[:,:,ω] = Hermitian(C[:,:,ω],:L)
-    end
     return FreqArray(C,fc)
 end
 
@@ -85,7 +82,7 @@ function csm_slow(t::AbstractArray{T};n=1024,noverlap=div(n,2),fs=1,win=DSP.hann
     csm_slow(flat_t(t);n=n,noverlap=noverlap,fs=fs,win=win,scaling=scaling)
 end
 
-function csm_slow(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat
+function csm_slow(t::Vector{<:AbstractVector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat
     M = length(t)
     Nf = div(n,2)+1
     nout = div((length(t[1]) - n), n - noverlap)+1
@@ -94,14 +91,19 @@ function csm_slow(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.han
     C = Array{Complex{T}}(undef,M,M,Nf)
     S = Array{Complex{T}}(undef,Nf,nout)
     Sc = similar(S)
+    
     for m in 1:M
-        stft!(Sc,t[m],n,noverlap; fs=fs, window=win, onesided=true)
-        conj!(Sc)
+        stft!(Sc, t[m], n, noverlap; fs=fs, window=win, onesided=true)  
+        conj!(Sc) 
+    
         for j in m:M
-            stft!(S,t[j],n,noverlap; fs=fs, window=win, onesided=true)
-            C[j,m,:] .= dropdims(mean(LazyArray(@~ Sc.*S),dims=2);dims=2)
+            stft!(S, t[j], n, noverlap; fs=fs, window=win, onesided=true) 
+            product = Sc .* S
+            mean_product = mean(product, dims=2)
+            C[j, m, :] .= dropdims(mean_product; dims=2)  # Drop the singleton dimension
         end
     end
+    
 
     if scaling == "density"
         scale = 1/(fs*sum(abs2,win))