diff --git a/src/DataReduction/SOD.jl b/src/DataReduction/SOD.jl
new file mode 100644
index 0000000..b99fb19
--- /dev/null
+++ b/src/DataReduction/SOD.jl
@@ -0,0 +1,243 @@
+"""
+    SOD <: AbstractDRProblem
+
+Smooth Orthogonal Decomposition (SOD) for extracting linear normal modes and natural
+frequencies from time-series data.
+
+SOD extends Proper Orthogonal Decomposition (POD) by considering temporal smoothness
+in addition to spatial variance. It uses the covariance matrices of displacement and
+velocity responses to form a generalized eigenvalue problem, which yields modes and
+frequencies without requiring knowledge of the system's mass matrix.
+
+# Fields
+- `snapshots`: Displacement snapshot data (matrix or vector of vectors), size (n_dof, n_samples)
+- `velocities`: Velocity snapshot data (same format as snapshots), or `:auto` to compute from snapshots
+- `dt`: Time step for computing velocities when `velocities = :auto`
+- `min_renergy`: Minimum relative energy threshold (0.0 to 1.0)
+- `min_nmodes`: Minimum number of modes to retain
+- `max_nmodes`: Maximum number of modes to retain
+- `nmodes`: Number of modes (computed after reduction)
+- `rbasis`: Reduced basis matrix (smooth orthogonal modes)
+- `renergy`: Relative energy captured by the modes
+- `spectrum`: Smooth orthogonal values (SOVs), which approximate ω² (squared frequencies)
+- `frequencies`: Estimated natural frequencies (√spectrum)
+
+# Constructor
+    SOD(snapshots; velocities=:auto, dt=1.0, min_renergy=1.0, min_nmodes=1, max_nmodes=size(snapshots,1))
+    SOD(snapshots, nmodes; velocities=:auto, dt=1.0)
+
+# References
+- Chelidze D, Zhou W. Smooth orthogonal decomposition-based vibration mode identification.
+  Journal of Sound and Vibration. 2006; 292(3-5): 461-473.
+  https://www.sciencedirect.com/science/article/abs/pii/S0022460X05005948
+"""
+mutable struct SOD <: AbstractDRProblem
+    # specified
+    snapshots::Any
+    velocities::Any
+    dt::Float64
+    min_renergy::Any
+    min_nmodes::Int
+    max_nmodes::Int
+    # computed
+    nmodes::Int
+    rbasis::Any
+    renergy::Any
+    spectrum::Any
+    frequencies::Any
+    # constructors
+    function SOD(
+            snaps;
+            velocities = :auto,
+            dt::Real = 1.0,
+            min_renergy = 1.0,
+            min_nmodes::Int = 1,
+            max_nmodes::Int = _get_ndof(snaps)
+        )
+        nmodes = min_nmodes
+        errorhandle_sod(snaps, velocities, dt, nmodes, min_renergy, min_nmodes, max_nmodes)
+        return new(
+            snaps, velocities, Float64(dt), min_renergy, min_nmodes, max_nmodes,
+            nmodes, missing, 1.0, missing, missing
+        )
+    end
+    function SOD(snaps, nmodes::Int; velocities = :auto, dt::Real = 1.0)
+        errorhandle_sod(snaps, velocities, dt, nmodes, 0.0, nmodes, nmodes)
+        return new(
+            snaps, velocities, Float64(dt), 0.0, nmodes, nmodes,
+            nmodes, missing, 1.0, missing, missing
+        )
+    end
+end
+
+# Helper to get number of degrees of freedom from snapshots
+function _get_ndof(snaps::AbstractMatrix)
+    return size(snaps, 1)
+end
+
+function _get_ndof(snaps::AbstractVector{<:AbstractVector})
+    return length(snaps[1])
+end
+
+# Error handling for SOD
+function errorhandle_sod(
+        data::AbstractMatrix, velocities, dt, nmodes::Int, max_energy,
+        min_nmodes, max_nmodes
+    )
+    @assert size(data, 1) > 1 "State vector is expected to be vector valued."
+    s = minimum(size(data))
+    @assert 0 < nmodes <= s "Number of modes should be in {1,2,...,$s}."
+    @assert min_nmodes <= max_nmodes "Minimum number of modes must lie below maximum number of modes"
+    @assert 0.0 <= max_energy <= 1.0 "Maximum relative energy must be in [0,1]"
+    @assert dt > 0 "Time step dt must be positive"
+    if velocities !== :auto
+        @assert size(velocities) == size(data) "Velocities must have same size as snapshots"
+    end
+    return nothing
+end
+
+function errorhandle_sod(
+        data::AbstractVector{T}, velocities, dt, nmodes::Int, max_energy, min_nmodes,
+        max_nmodes
+    ) where {T <: AbstractVector}
+    @assert size(data[1], 1) > 1 "State vector is expected to be vector valued."
+    s = min(size(data, 1), size(data[1], 1))
+    @assert 0 < nmodes <= s "Number of modes should be in {1,2,...,$s}."
+    @assert min_nmodes <= max_nmodes "Minimum number of modes must lie below maximum number of modes"
+    @assert 0.0 <= max_energy <= 1.0 "Maximum relative energy must be in [0,1]"
+    @assert dt > 0 "Time step dt must be positive"
+    if velocities !== :auto
+        @assert length(velocities) == length(data) "Velocities must have same length as snapshots"
+    end
+    return nothing
+end
+
+# Compute velocity from displacement using finite differences
+function _compute_velocity(X::AbstractMatrix, dt::Real)
+    # Central difference for interior points, forward/backward for boundaries
+    n_dof, n_samples = size(X)
+    V = similar(X)
+
+    if n_samples == 1
+        V .= 0
+        return V
+    elseif n_samples == 2
+        # Forward difference only
+        V[:, 1] = (X[:, 2] - X[:, 1]) / dt
+        V[:, 2] = V[:, 1]
+        return V
+    end
+
+    # Forward difference for first point
+    V[:, 1] = (X[:, 2] - X[:, 1]) / dt
+
+    # Central difference for interior points
+    for i in 2:(n_samples - 1)
+        V[:, i] = (X[:, i + 1] - X[:, i - 1]) / (2 * dt)
+    end
+
+    # Backward difference for last point
+    V[:, n_samples] = (X[:, n_samples] - X[:, n_samples - 1]) / dt
+
+    return V
+end
+
+function _compute_velocity(X::AbstractVector{<:AbstractVector}, dt::Real)
+    X_mat = matricize(X)
+    V_mat = _compute_velocity(X_mat, dt)
+    return V_mat
+end
+
+# Convert velocities to matrix form
+function _get_velocity_matrix(sod::SOD)
+    if sod.velocities === :auto
+        X = sod.snapshots isa AbstractMatrix ? sod.snapshots : matricize(sod.snapshots)
+        return _compute_velocity(X, sod.dt)
+    else
+        return sod.velocities isa AbstractMatrix ? sod.velocities : matricize(sod.velocities)
+    end
+end
+
+# Get displacement matrix
+function _get_displacement_matrix(sod::SOD)
+    return sod.snapshots isa AbstractMatrix ? sod.snapshots : matricize(sod.snapshots)
+end
+
+"""
+    reduce!(sod::SOD)
+
+Perform Smooth Orthogonal Decomposition on the snapshot data.
+
+This solves the generalized eigenvalue problem:
+    Σ_xx Ψ = λ Σ_vv Ψ
+
+where Σ_xx = X X^T / N (displacement covariance) and Σ_vv = V V^T / N (velocity covariance).
+
+The eigenvalues λ (smooth orthogonal values, SOVs) approximate ω² (squared natural frequencies).
+The eigenvectors Ψ are the smooth orthogonal modes (SOMs).
+"""
+function reduce!(sod::SOD)
+    X = _get_displacement_matrix(sod)
+    V = _get_velocity_matrix(sod)
+
+    n_dof, n_samples = size(X)
+
+    # Compute covariance matrices
+    # Σ_xx = X * X' / n_samples (displacement covariance)
+    # Σ_vv = V * V' / n_samples (velocity covariance)
+    Σ_xx = X * X' / n_samples
+    Σ_vv = V * V' / n_samples
+
+    # Add small regularization to Σ_vv if it's singular
+    # This handles cases where velocity variance is very small
+    eps_reg = 1.0e-10 * maximum(abs, Σ_vv)
+    if eps_reg < 1.0e-14
+        eps_reg = 1.0e-14
+    end
+    Σ_vv_reg = Σ_vv + eps_reg * I
+
+    # Solve generalized eigenvalue problem: Σ_xx Ψ = λ Σ_vv Ψ
+    # This is equivalent to finding eigenvalues of Σ_vv^(-1) Σ_xx
+    # For numerical stability, we use eigen on the symmetric pencil
+    eig_result = eigen(Symmetric(Σ_xx), Symmetric(Σ_vv_reg))
+
+    # Eigenvalues are the SOVs (smooth orthogonal values)
+    # They should be sorted in decreasing order (largest eigenvalue = highest energy)
+    sovs = real.(eig_result.values)
+    soms = real.(eig_result.vectors)
+
+    # Sort by eigenvalue magnitude (largest first)
+    perm = sortperm(sovs, rev = true)
+    sovs = sovs[perm]
+    soms = soms[:, perm]
+
+    # Determine truncation based on energy criterion
+    sod.nmodes, sod.renergy = determine_truncation(
+        sovs, sod.min_nmodes, sod.min_renergy, sod.max_nmodes
+    )
+
+    # Store results
+    sod.rbasis = soms[:, 1:(sod.nmodes)]
+    sod.spectrum = sovs
+
+    # Compute natural frequencies from SOVs
+    # SOVs approximate ω², so frequencies are √SOVs
+    # Only take sqrt of positive values
+    sod.frequencies = sqrt.(max.(sovs[1:(sod.nmodes)], 0.0))
+
+    return nothing
+end
+
+function Base.show(io::IO, sod::SOD)
+    print(io, "SOD \n")
+    print(io, "Reduction Order = ", sod.nmodes, "\n")
+    n_dof = _get_ndof(sod.snapshots)
+    n_samples = sod.snapshots isa AbstractMatrix ? size(sod.snapshots, 2) :
+        length(sod.snapshots)
+    print(io, "Snapshot size = (", n_dof, ",", n_samples, ")\n")
+    print(io, "Relative Energy = ", sod.renergy, "\n")
+    if !ismissing(sod.frequencies)
+        print(io, "Estimated frequencies = ", sod.frequencies, "\n")
+    end
+    return nothing
+end
diff --git a/src/ModelOrderReduction.jl b/src/ModelOrderReduction.jl
index 270084e..d06eae2 100644
--- a/src/ModelOrderReduction.jl
+++ b/src/ModelOrderReduction.jl
@@ -5,7 +5,7 @@ using DocStringExtensions: DocStringExtensions, FUNCTIONNAME, SIGNATURES, TYPEDS
 using ModelingToolkit: ModelingToolkit, @variables, Differential, Equation, Num, ODESystem,
     SymbolicUtils, Symbolics, arguments, build_function, complete,
     expand, substitute, tearing_substitution
-using LinearAlgebra: LinearAlgebra, /, \, mul!, qr, svd
+using LinearAlgebra: LinearAlgebra, /, \, I, Symmetric, eigen, mul!, qr, svd
 
 using Setfield: Setfield, @set!
 
@@ -15,8 +15,9 @@ include("Types.jl")
 include("ErrorHandle.jl")
 
 include("DataReduction/POD.jl")
+include("DataReduction/SOD.jl")
 export SVD, TSVD, RSVD
-export POD, reduce!
+export POD, SOD, reduce!
 
 include("deim.jl")
 export deim
diff --git a/test/SOD.jl b/test/SOD.jl
new file mode 100644
index 0000000..a4e2c54
--- /dev/null
+++ b/test/SOD.jl
@@ -0,0 +1,139 @@
+using Test, ModelOrderReduction
+using OrdinaryDiffEq
+
+# Simple harmonic oscillator system for testing SOD
+# This system has known natural frequencies that SOD should recover
+function harmonic_oscillator_prob()
+    # 2-DOF mass-spring system with known frequencies
+    # m1*x1'' + k1*x1 + k2*(x1-x2) = 0
+    # m2*x2'' + k2*(x2-x1) = 0
+    # With m1=m2=1, k1=1, k2=1, natural frequencies are ω₁ ≈ 0.618, ω₂ ≈ 1.618
+
+    function harmonic!(du, u, p, t)
+        # u = [x1, v1, x2, v2]
+        k1, k2, m1, m2 = p
+        x1, v1, x2, v2 = u
+        du[1] = v1
+        du[2] = (-k1 * x1 - k2 * (x1 - x2)) / m1
+        du[3] = v2
+        return du[4] = (-k2 * (x2 - x1)) / m2
+    end
+
+    u0 = [1.0, 0.0, 0.5, 0.0]  # Initial displacement, zero velocity
+    p = [1.0, 1.0, 1.0, 1.0]   # k1, k2, m1, m2
+    tspan = (0.0, 50.0)
+    prob = ODEProblem(harmonic!, u0, tspan, p)
+    sol = solve(prob, Tsit5(), saveat = 0.01)
+    return sol
+end
+
+# Use Lorenz system for comparison with POD (same test as POD uses)
+function lorenz_prob()
+    function lorenz!(du, u, p, t)
+        du[1] = p[1] * (u[2] - u[1])
+        du[2] = u[1] * (p[2] - u[3]) - u[2]
+        return du[3] = u[1] * u[2] - p[3] * u[3]
+    end
+
+    u0 = [1, 0, 0]
+    p = [10, 28, 8 / 3]
+    tspan = (0, 100)
+    prob = ODEProblem(lorenz!, u0, tspan, p)
+    sol = solve(prob, Tsit5(), saveat = 0.1)
+    return sol
+end
+
+@testset "SOD - Basic Functionality" begin
+    sol = lorenz_prob()
+    solution = Array(sol)
+    dt = 0.1
+
+    # Test with matrix input
+    order = 2
+    matrix_reducer = SOD(solution, order; dt = dt)
+    reduce!(matrix_reducer)
+
+    @test size(matrix_reducer.rbasis, 2) == matrix_reducer.nmodes
+    @test size(matrix_reducer.rbasis, 1) == size(solution, 1)
+    @test !ismissing(matrix_reducer.spectrum)
+    @test !ismissing(matrix_reducer.frequencies)
+    @test length(matrix_reducer.frequencies) == order
+
+    # Test with vector of vectors input
+    snapshot_reducer = SOD(sol.u, order; dt = dt)
+    reduce!(snapshot_reducer)
+
+    @test size(snapshot_reducer.rbasis, 2) == snapshot_reducer.nmodes
+    @test !ismissing(snapshot_reducer.frequencies)
+end
+
+@testset "SOD - Energy Truncation" begin
+    sol = lorenz_prob()
+    solution = Array(sol)
+    dt = 0.1
+
+    # Test with energy-based truncation
+    reducer = SOD(solution; dt = dt, min_nmodes = 1, max_nmodes = 3, min_renergy = 0.1)
+    reduce!(reducer)
+    @test reducer.renergy >= 0.1
+    @test 1 <= reducer.nmodes <= 3
+end
+
+@testset "SOD - Provided Velocities" begin
+    sol = lorenz_prob()
+    solution = Array(sol)
+    dt = 0.1
+
+    # Manually compute velocities (central difference)
+    n_dof, n_samples = size(solution)
+    velocities = similar(solution)
+    velocities[:, 1] = (solution[:, 2] - solution[:, 1]) / dt
+    for i in 2:(n_samples - 1)
+        velocities[:, i] = (solution[:, i + 1] - solution[:, i - 1]) / (2 * dt)
+    end
+    velocities[:, n_samples] = (solution[:, n_samples] - solution[:, n_samples - 1]) / dt
+
+    # Test with provided velocities
+    order = 2
+    reducer = SOD(solution, order; velocities = velocities, dt = dt)
+    reduce!(reducer)
+
+    @test size(reducer.rbasis, 2) == reducer.nmodes
+    @test !ismissing(reducer.frequencies)
+end
+
+@testset "SOD - Harmonic Oscillator Frequency Recovery" begin
+    sol = harmonic_oscillator_prob()
+
+    # Extract position data (x1 and x2, not velocities)
+    positions = vcat(sol[1, :]', sol[3, :]')  # Shape: (2, n_samples)
+    dt = 0.01
+
+    # SOD should find the modes of this oscillating system
+    reducer = SOD(positions, 2; dt = dt)
+    reduce!(reducer)
+
+    @test reducer.nmodes == 2
+    @test !ismissing(reducer.frequencies)
+
+    # Check that we get reasonable positive frequencies
+    # (The exact values depend on the sampling and system dynamics)
+    @test all(reducer.frequencies .>= 0)
+end
+
+@testset "SOD - Show Method" begin
+    sol = lorenz_prob()
+    solution = Array(sol)
+    dt = 0.1
+
+    reducer = SOD(solution, 2; dt = dt)
+    reduce!(reducer)
+
+    # Test that show doesn't error
+    io = IOBuffer()
+    show(io, reducer)
+    str = String(take!(io))
+    @test occursin("SOD", str)
+    @test occursin("Reduction Order", str)
+    @test occursin("frequencies", str)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 6d1c1af..0e6e978 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,6 +9,9 @@ end
 @safetestset "POD" begin
     include("DataReduction.jl")
 end
+@safetestset "SOD" begin
+    include("SOD.jl")
+end
 @safetestset "utils" begin
     include("utils.jl")
 end