diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml
new file mode 100644
index 0000000..d9871b9
--- /dev/null
+++ b/.github/workflows/Documenter.yml
@@ -0,0 +1,24 @@
+name: Documentation
+
+on:
+  push:
+    branches:
+      - master
+    tags: ['*']
+  pull_request:
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: '1'
+      - uses: julia-actions/cache@v2
+      - name: Install dependencies
+        run: julia --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+      - name: Build and deploy
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+        run: julia --project=docs docs/make.jl
diff --git a/.gitignore b/.gitignore
index ba39cc5..5bfcd30 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,2 @@
 Manifest.toml
+docs/build/
diff --git a/Project.toml b/Project.toml
index ac51e47..c6e77d6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,16 +6,19 @@ version = "0.2.0"
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
+Aqua = "0.8"
 DocStringExtensions = "0.8, 0.9"
+LinearAlgebra = "1"
+Test = "1"
 Unitful = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 1.0"
 julia = "1"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-[target]
-test = ["Test"]
+[targets]
+test = ["Aqua", "Test"]
diff --git a/README.md b/README.md
index be48d2f..d1b1ede 100644
--- a/README.md
+++ b/README.md
@@ -1,26 +1,22 @@
-[![Build Status](https://travis-ci.org/JuliaGNSS/TrackingLoopFilters.jl.svg?branch=master)](https://travis-ci.org/JuliaGNSS/TrackingLoopFilters.jl)
-[![Coverage Status](https://coveralls.io/repos/github/JuliaGNSS/TrackingLoopFilters.jl/badge.svg?branch=master)](https://coveralls.io/github/JuliaGNSS/TrackingLoopFilters.jl?branch=master)
+# TrackingLoopFilters.jl
 
-# TrackingLoopFilters
-This implements multiple loop filters for GNSS tracking algorithms.
+[![Build Status](https://github.com/JuliaGNSS/TrackingLoopFilters.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/JuliaGNSS/TrackingLoopFilters.jl/actions/workflows/ci.yml)
+[![codecov](https://codecov.io/gh/JuliaGNSS/TrackingLoopFilters.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaGNSS/TrackingLoopFilters.jl)
+[![Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliagnss.github.io/TrackingLoopFilters.jl/stable)
+[![Documentation](https://img.shields.io/badge/docs-dev-blue.svg)](https://juliagnss.github.io/TrackingLoopFilters.jl/dev)
 
-## Features
+Loop filters for GNSS tracking algorithms.
 
-* First, second and third order loop filters
-* Boxcar and bilinear loop filters
+## Installation
 
-## Getting started
-
-Install:
 ```julia
-julia> ]
 pkg> add TrackingLoopFilters
 ```
 
-## Usage
+## Documentation
 
-```julia
-using TrackingLoopFilters
-carrier_loop_filter = ThirdOrderBilinearLF()
-output, next_carrier_loop_filter = filter_loop(carrier_loop_filter, discriminator_output, Δt, bandwidth)
-```
+See the [documentation](https://juliagnss.github.io/TrackingLoopFilters.jl/stable) for API reference.
+
+## License
+
+MIT
diff --git a/docs/Project.toml b/docs/Project.toml
new file mode 100644
index 0000000..2db9200
--- /dev/null
+++ b/docs/Project.toml
@@ -0,0 +1,6 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+TrackingLoopFilters = "0814aff9-93cb-554c-9fff-9bf1cfdb5efa"
+
+[compat]
+Documenter = "1"
diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000..b24bc98
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,14 @@
+using Documenter
+using TrackingLoopFilters
+
+makedocs(
+    sitename = "TrackingLoopFilters.jl",
+    modules = [TrackingLoopFilters],
+    pages = [
+        "Home" => "index.md",
+    ],
+)
+
+deploydocs(
+    repo = "github.com/JuliaGNSS/TrackingLoopFilters.jl.git",
+)
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000..b79780e
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,28 @@
+# TrackingLoopFilters.jl
+
+Loop filters for GNSS tracking algorithms.
+
+## Installation
+
+```julia
+pkg> add TrackingLoopFilters
+```
+
+## Quick Start
+
+```julia
+using TrackingLoopFilters
+using Unitful: Hz, s
+
+# Create a loop filter
+lf = ThirdOrderBilinearLF()
+
+# Process discriminator output
+output, next_lf = filter_loop(lf, δθ, Δt, bandwidth)
+```
+
+## API Reference
+
+```@autodocs
+Modules = [TrackingLoopFilters]
+```
diff --git a/src/FirstOrderLF.jl b/src/FirstOrderLF.jl
index f287c47..e402764 100644
--- a/src/FirstOrderLF.jl
+++ b/src/FirstOrderLF.jl
@@ -1,12 +1,42 @@
+"""
+$(TYPEDEF)
+
+Abstract base type for first order loop filters.
+"""
 abstract type AbstractFirstOrderLF <: AbstractLoopFilter end
 
+"""
+$(TYPEDEF)
+
+First order loop filter (proportional only).
+
+A stateless filter that provides proportional gain without integration.
+The natural frequency scaling factor is 4.0.
+
+# Example
+```julia
+lf = FirstOrderLF()
+output, next_lf = filter_loop(lf, δθ, Δt, bandwidth)
+```
+"""
 struct FirstOrderLF <: AbstractFirstOrderLF
 end
 
 """
 $(SIGNATURES)
 
-Propagates the state of the loop filter.
+Propagate the first order loop filter state.
+
+Since the first order filter is stateless, this returns the input state unchanged.
+
+# Arguments
+- `state`: Current loop filter state
+- `δθ`: Phase discriminator output
+- `Δt`: Integration time
+- `bandwidth`: Loop bandwidth
+
+# Returns
+The unchanged loop filter state.
 """
 function propagate(state::FirstOrderLF, δθ, Δt, bandwidth)
     state
@@ -15,7 +45,18 @@ end
 """
 $(SIGNATURES)
 
-Calculates the output of the loop filter.
+Calculate the filtered output for the first order loop filter.
+
+Computes `ω₀ * δθ` where `ω₀ = bandwidth * 4.0`.
+
+# Arguments
+- `state`: Current loop filter state
+- `δθ`: Phase discriminator output
+- `Δt`: Integration time
+- `bandwidth`: Loop bandwidth
+
+# Returns
+Filtered frequency estimate.
 """
 function get_filtered_output(state::FirstOrderLF, δθ, Δt, bandwidth)
     ω₀ = bandwidth * 4.0
diff --git a/src/SecondOrderLF.jl b/src/SecondOrderLF.jl
index 488ec38..20c8bce 100644
--- a/src/SecondOrderLF.jl
+++ b/src/SecondOrderLF.jl
@@ -1,17 +1,66 @@
+"""
+$(TYPEDEF)
+
+Abstract base type for second order loop filters.
+"""
 abstract type AbstractSecondOrderLF <: AbstractLoopFilter end
 
+"""
+$(TYPEDEF)
+
+Second order bilinear loop filter.
+
+Uses bilinear transformation for improved frequency response.
+The natural frequency scaling factor is 1.89.
+
+$(TYPEDFIELDS)
+
+# Example
+```julia
+lf = SecondOrderBilinearLF()
+output, next_lf = filter_loop(lf, δθ, Δt, bandwidth)
+```
+"""
 struct SecondOrderBilinearLF{T} <: AbstractSecondOrderLF
+    "Frequency state estimate"
     x::T
 end
 
+"""
+$(TYPEDEF)
+
+Second order boxcar loop filter.
+
+Uses boxcar (rectangular) integration for simpler implementation.
+The natural frequency scaling factor is 1.89.
+
+$(TYPEDFIELDS)
+
+# Example
+```julia
+lf = SecondOrderBoxcarLF()
+output, next_lf = filter_loop(lf, δθ, Δt, bandwidth)
+```
+"""
 struct SecondOrderBoxcarLF{T} <: AbstractSecondOrderLF
+    "Frequency state estimate"
     x::T
 end
 
+"""
+    SecondOrderBilinearLF()
+
+Construct a second order bilinear loop filter with zero initial state.
+"""
 function SecondOrderBilinearLF()
     SecondOrderBilinearLF(0.0Hz)
 end
 
+"""
+    SecondOrderBoxcarLF()
+
+Construct a second order boxcar loop filter with zero initial state.
+"""
 function SecondOrderBoxcarLF()
     SecondOrderBoxcarLF(0.0Hz)
 end
@@ -19,7 +68,18 @@ end
 """
 $(SIGNATURES)
 
-Propagates the state of the loop filter.
+Propagate the second order loop filter state.
+
+Updates the frequency state estimate using `x_next = x + Δt * ω₀² * δθ`.
+
+# Arguments
+- `state`: Current loop filter state
+- `δθ`: Phase discriminator output
+- `Δt`: Integration time
+- `bandwidth`: Loop bandwidth
+
+# Returns
+New loop filter state with updated frequency estimate.
 """
 function propagate(state::T, δθ, Δt, bandwidth) where T <: AbstractSecondOrderLF
     ω₀ = bandwidth * 1.89
@@ -29,7 +89,18 @@ end
 """
 $(SIGNATURES)
 
-Calculates the output of the loop filter.
+Calculate the filtered output for the second order bilinear loop filter.
+
+Computes `x + (√2 * ω₀ + ω₀² * Δt / 2) * δθ`.
+
+# Arguments
+- `state`: Current loop filter state
+- `δθ`: Phase discriminator output
+- `Δt`: Integration time
+- `bandwidth`: Loop bandwidth
+
+# Returns
+Filtered frequency estimate.
 """
 function get_filtered_output(state::SecondOrderBilinearLF, δθ, Δt, bandwidth)
     ω₀ = bandwidth * 1.89
@@ -39,7 +110,18 @@ end
 """
 $(SIGNATURES)
 
-Calculates the output of the loop filter.
+Calculate the filtered output for the second order boxcar loop filter.
+
+Computes `x + √2 * ω₀ * δθ`.
+
+# Arguments
+- `state`: Current loop filter state
+- `δθ`: Phase discriminator output
+- `Δt`: Integration time
+- `bandwidth`: Loop bandwidth
+
+# Returns
+Filtered frequency estimate.
 """
 function get_filtered_output(state::SecondOrderBoxcarLF, δθ, Δt, bandwidth)
     ω₀ = bandwidth * 1.89
diff --git a/src/ThirdOrderLF.jl b/src/ThirdOrderLF.jl
index 9ea4dd1..eb1b870 100644
--- a/src/ThirdOrderLF.jl
+++ b/src/ThirdOrderLF.jl
@@ -1,30 +1,110 @@
+"""
+$(TYPEDEF)
 
+Abstract base type for third order loop filters.
+"""
 abstract type AbstractThirdOrderLF <: AbstractLoopFilter end
+
+"""
+$(TYPEDEF)
+
+Abstract base type for assisted third order loop filters.
+"""
 abstract type AbstractThirdOrderAssistedLF <: AbstractLoopFilter end
 
+"""
+$(TYPEDEF)
+
+Third order bilinear loop filter.
+
+Uses bilinear transformation for improved frequency response.
+The natural frequency scaling factor is 1.2.
+
+$(TYPEDFIELDS)
+
+# Example
+```julia
+lf = ThirdOrderBilinearLF()
+output, next_lf = filter_loop(lf, δθ, Δt, bandwidth)
+```
+"""
 struct ThirdOrderBilinearLF{T1,T2} <: AbstractThirdOrderLF
+    "Frequency state estimate"
     x1::T1
+    "Frequency rate state estimate"
     x2::T2
 end
 
+"""
+$(TYPEDEF)
+
+Third order bilinear loop filter with second order assistance.
+
+Combines a third order bilinear loop with a second order assisted loop
+for improved tracking performance. Accepts a two-element discriminator
+input vector `[δθ_high, δθ_low]`.
+
+$(TYPEDFIELDS)
+
+# Example
+```julia
+lf = ThirdOrderAssistedBilinearLF()
+output, next_lf = filter_loop(lf, [δθ_high, δθ_low], Δt, bandwidth)
+```
+"""
 struct ThirdOrderAssistedBilinearLF{T1,T2} <: AbstractThirdOrderAssistedLF
+    "Frequency state estimate"
     x1::T1
+    "Frequency rate state estimate"
     x2::T2
 end
 
+"""
+$(TYPEDEF)
+
+Third order boxcar loop filter.
+
+Uses boxcar (rectangular) integration for simpler implementation.
+The natural frequency scaling factor is 1.2.
+
+$(TYPEDFIELDS)
+
+# Example
+```julia
+lf = ThirdOrderBoxcarLF()
+output, next_lf = filter_loop(lf, δθ, Δt, bandwidth)
+```
+"""
 struct ThirdOrderBoxcarLF{T1,T2} <: AbstractThirdOrderLF
+    "Frequency state estimate"
     x1::T1
+    "Frequency rate state estimate"
     x2::T2
 end
 
+"""
+    ThirdOrderBilinearLF()
+
+Construct a third order bilinear loop filter with zero initial state.
+"""
 function ThirdOrderBilinearLF()
     ThirdOrderBilinearLF(0.0Hz, 0.0Hz^2)
 end
 
+"""
+    ThirdOrderAssistedBilinearLF()
+
+Construct a third order assisted bilinear loop filter with zero initial state.
+"""
 function ThirdOrderAssistedBilinearLF()
     ThirdOrderAssistedBilinearLF(0.0Hz, 0.0Hz^2)
 end
 
+"""
+    ThirdOrderBoxcarLF()
+
+Construct a third order boxcar loop filter with zero initial state.
+"""
 function ThirdOrderBoxcarLF()
     ThirdOrderBoxcarLF(0.0Hz, 0.0Hz^2)
 end
@@ -32,7 +112,18 @@ end
 """
 $(SIGNATURES)
 
-Propagates the state of the loop filter.
+Propagate the third order loop filter state.
+
+Updates both frequency and frequency rate state estimates.
+
+# Arguments
+- `state`: Current loop filter state
+- `δθ`: Phase discriminator output
+- `Δt`: Integration time
+- `bandwidth`: Loop bandwidth
+
+# Returns
+New loop filter state with updated estimates.
 """
 function propagate(state::T, δθ, Δt, bandwidth) where T <: AbstractThirdOrderLF
     ω₀ = bandwidth * 1.2
@@ -42,7 +133,19 @@ end
 """
 $(SIGNATURES)
 
-Propagates the state of the assisted loop filter.
+Propagate the assisted third order loop filter state.
+
+Updates both frequency and frequency rate state estimates using
+dual discriminator inputs.
+
+# Arguments
+- `state`: Current loop filter state
+- `δθ`: Two-element vector `[δθ_high, δθ_low]` with high and low order discriminator outputs
+- `Δt`: Integration time
+- `bandwidth`: Loop bandwidth
+
+# Returns
+New loop filter state with updated estimates.
 """
 function propagate(state::T, δθ, Δt, bandwidth) where T <: AbstractThirdOrderAssistedLF
     ω₀ = bandwidth * 1.2
@@ -53,7 +156,16 @@ end
 """
 $(SIGNATURES)
 
-Calculates the output of the loop filter.
+Calculate the filtered output for the third order bilinear loop filter.
+
+# Arguments
+- `state`: Current loop filter state
+- `δθ`: Phase discriminator output
+- `Δt`: Integration time
+- `bandwidth`: Loop bandwidth
+
+# Returns
+Filtered frequency estimate.
 """
 function get_filtered_output(state::ThirdOrderBilinearLF, δθ, Δt, bandwidth)
     ω₀= bandwidth * 1.2
@@ -63,7 +175,18 @@ end
 """
 $(SIGNATURES)
 
-Calculates the output of the third order loop filter assisted by a second order loop filter.
+Calculate the filtered output for the assisted third order bilinear loop filter.
+
+Combines outputs from the third order loop and second order assisted loop.
+
+# Arguments
+- `state`: Current loop filter state
+- `δθ`: Two-element vector `[δθ_high, δθ_low]` with high and low order discriminator outputs
+- `Δt`: Integration time
+- `bandwidth`: Loop bandwidth
+
+# Returns
+Filtered frequency estimate.
 """
 function get_filtered_output(state::ThirdOrderAssistedBilinearLF, δθ, Δt, bandwidth)
     ω₀= bandwidth * 1.2
@@ -76,7 +199,16 @@ end
 """
 $(SIGNATURES)
 
-Calculates the output of the loop filter.
+Calculate the filtered output for the third order boxcar loop filter.
+
+# Arguments
+- `state`: Current loop filter state
+- `δθ`: Phase discriminator output
+- `Δt`: Integration time
+- `bandwidth`: Loop bandwidth
+
+# Returns
+Filtered frequency estimate.
 """
 function get_filtered_output(state::ThirdOrderBoxcarLF, δθ, Δt, bandwidth)
     ω₀= bandwidth * 1.2
diff --git a/src/TrackingLoopFilters.jl b/src/TrackingLoopFilters.jl
index c1bf225..448deae 100644
--- a/src/TrackingLoopFilters.jl
+++ b/src/TrackingLoopFilters.jl
@@ -1,3 +1,24 @@
+"""
+    TrackingLoopFilters
+
+A Julia package implementing loop filters for GNSS tracking algorithms.
+
+Provides first, second, and third order loop filters with boxcar and bilinear
+implementations. All filters support Unitful quantities for type-safe calculations.
+
+# Exported Types
+- [`FirstOrderLF`](@ref): First order loop filter (stateless)
+- [`SecondOrderBilinearLF`](@ref): Second order bilinear loop filter
+- [`SecondOrderBoxcarLF`](@ref): Second order boxcar loop filter
+- [`ThirdOrderBilinearLF`](@ref): Third order bilinear loop filter
+- [`ThirdOrderBoxcarLF`](@ref): Third order boxcar loop filter
+- [`ThirdOrderAssistedBilinearLF`](@ref): Third order bilinear loop filter with second order assistance
+
+# Exported Functions
+- [`propagate`](@ref): Propagate the loop filter state
+- [`get_filtered_output`](@ref): Get the filtered output for the current state
+- [`filter_loop`](@ref): Combined propagate and get output in one call
+"""
 module TrackingLoopFilters
 
     using DocStringExtensions
@@ -16,6 +37,11 @@ module TrackingLoopFilters
         ThirdOrderBoxcarLF,
         AbstractLoopFilter
 
+    """
+    $(TYPEDEF)
+
+    Abstract base type for all loop filters.
+    """
     abstract type AbstractLoopFilter end
 
     include("FirstOrderLF.jl")
@@ -25,8 +51,20 @@ module TrackingLoopFilters
     """
     $(SIGNATURES)
 
-    Propagates the state of loop filter and returns the calculated output of the loop filter
-    as well as the propagated loop filter state.
+    Propagate the loop filter state and return both the filtered output and next state.
+
+    This is a convenience function that combines [`propagate`](@ref) and
+    [`get_filtered_output`](@ref) into a single call.
+
+    # Arguments
+    - `state`: Current loop filter state
+    - `δθ`: Phase discriminator output (dimensionless for standard filters, vector for assisted)
+    - `Δt`: Integration time
+    - `bandwidth`: Loop bandwidth
+
+    # Returns
+    - `output`: Filtered frequency estimate
+    - `next_state`: Propagated loop filter state
     """
     function filter_loop(state::T, δθ, Δt, bandwidth) where T <: AbstractLoopFilter
         next_state = propagate(state, δθ, Δt, bandwidth)
diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..3c24aa8
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,5 @@
+[deps]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
diff --git a/test/runtests.jl b/test/runtests.jl
index bec2524..f01b6d1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,9 +1,14 @@
+using Aqua
 using LinearAlgebra
 using Test
 using TrackingLoopFilters
 
 import Unitful: Hz, s
 
+@testset "Aqua" begin
+    Aqua.test_all(TrackingLoopFilters)
+end
+
 @testset "First Order Loop Filter" begin
     bandwidth = 1Hz
     loop_filter = @inferred FirstOrderLF()