diff --git a/Project.toml b/Project.toml
index c7706bb..aabfe9c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,6 +14,7 @@ PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
+UnitfulEquivalences = "da9c4bc3-91c8-4f02-8a40-6b990d2a7e0c"
 
 [compat]
 DustExtinction = ">=0.11.7"
@@ -23,4 +24,5 @@ PhysicalConstants = "0.2"
 RecipesBase = "1"
 Unitful = "1"
 UnitfulAstro = "1"
+UnitfulEquivalences = "0.2.1"
 julia = "1.10"
diff --git a/docs/src/index.md b/docs/src/index.md
index 1df9a08..f5bfb19 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -37,7 +37,7 @@ spec = spectrum(wave, flux)
 plot(spec)
 ```
 
-For constructing higher dimensional spectra, e.g., for echelle or IFU spectra, see the docstrings for [EchelleSpectrum](@ref) and [IFUSpectrum](@ref), respectively.
+For constructing higher dimensional spectra, e.g., for echelle or IFU spectra, see the docstrings for [`EchelleSpectrum`](@ref) and [`IFUSpectrum`](@ref), respectively.
 
 ## Citation
 
diff --git a/docs/src/spectrum.md b/docs/src/spectrum.md
index 315b621..297ef38 100644
--- a/docs/src/spectrum.md
+++ b/docs/src/spectrum.md
@@ -12,6 +12,7 @@ SpectrumBase.Spectrum
 SpectrumBase.SingleSpectrum
 SpectrumBase.EchelleSpectrum
 SpectrumBase.IFUSpectrum
+SpectrumBase.BinnedSpectrum
 ```
 
 ## Constructors
diff --git a/docs/src/transforms.md b/docs/src/transforms.md
index dbcdd4c..ee1c248 100644
--- a/docs/src/transforms.md
+++ b/docs/src/transforms.md
@@ -43,7 +43,7 @@ julia> spectral_axis(shifted) ≈ spectral_axis(spec) .* 1.5
 true
 ```
 
-**Doppler shift** shifts by a radial velocity `v`. Pass a `Unitful` velocity or a plain number (interpreted as m/s). Set `relativistic=true` for the full relativistic formula:
+**Doppler shift** shifts by a radial velocity `v`. Pass a `Unitful` velocity or a plain number (interpreted as km/s). Set `relativistic=true` for the full relativistic formula:
 
 ```jldoctest
 julia> using SpectrumBase, Unitful
@@ -57,6 +57,22 @@ julia> shifted_rel = doppler_shift(spec, 100u"km/s"; relativistic=true);
 
 Both `redshift` and `doppler_shift` return a new spectrum. In-place variants `redshift!` and `doppler_shift!` are also available.
 
+Binned spectra are supported as well: both bin edges are shifted and the bin values are carried unchanged. An energy axis is compressed rather than stretched:
+
+```jldoctest
+julia> using SpectrumBase, Unitful
+
+julia> spec = spectrum([1.0 2.0; 2.0 3.0; 3.0 4.0]u"keV", [10.0, 20.0, 30.0]u"erg/s/cm^2");
+
+julia> shifted = redshift(spec, 0.5);
+
+julia> spectral_axis(shifted) ≈ spectral_axis(spec) ./ 1.5
+true
+
+julia> flux_axis(shifted) == flux_axis(spec)
+true
+```
+
 ### API/Reference
 
 ```@docs
@@ -66,6 +82,32 @@ doppler_shift
 doppler_shift!
 ```
 
+## Axis conversions
+
+Powered by [UnitfulEquivalences.jl](https://github.com/sostock/UnitfulEquivalences.jl), [`uconvert`](@ref) from [Unitful.jl](https://github.com/JuliaPhysics/Unitful.jl) is extended to convert the spectral axis between wavelength, frequency, and photon energy, along with its associated flux density:
+
+```jldoctest
+julia> using SpectrumBase, Unitful, UnitfulAstro
+
+julia> spec = spectrum([1.0, 1.5, 2.0]u"μm", [1.0, 2.0, 3.0]u"W/m^2/μm");
+
+julia> uconvert(u"THz", spec);
+
+julia> νspec = uconvert((u"THz", u"Jy"), spec);
+
+julia> issorted(spectral_axis(νspec); rev = true)
+true
+
+julia> flux_axis(νspec)[1] ≈ uconvert(u"Jy", 1.0u"W/m^2/μm" * (1.0u"μm")^2 / Unitful.c0)
+true
+```
+
+### API/Reference
+
+```@docs
+Unitful.uconvert
+```
+
 ## Resampling
 
 External interpolators, e.g., from [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) or [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl), can be used to resample spectra onto a given wavelength grid. Starting with a sample spectrum `spec`, we first create a [`SpectrumResampler`](@ref) object `resampler` which stores the initial spectrum and interpolator `interp` together. We then apply this object to the wavelength grid of our choice to produce the resampled spectrum. We show example usage in the docstring below:
diff --git a/src/SpectrumBase.jl b/src/SpectrumBase.jl
index 061259f..94777d6 100644
--- a/src/SpectrumBase.jl
+++ b/src/SpectrumBase.jl
@@ -4,7 +4,7 @@ module SpectrumBase
 export AbstractSpectrum, Spectrum, spectrum, spectral_axis, flux_axis
 
 # AbstractSpectrum types
-export SingleSpectrum, IFUSpectrum, EchelleSpectrum
+export SingleSpectrum, IFUSpectrum, EchelleSpectrum, BinnedSpectrum
 
 # Transforms
 export SpectrumResampler, redden, redden!, deredden, deredden!
@@ -54,9 +54,7 @@ mutable struct Spectrum{S<:Number, F<:Number, M, N, A<:AbstractArray{S, M}, B<:A
             - SingleSpectrum: spectral axis (M-length vector), flux axis (M-length vector)
             - EchelleSpectrum: spectral axis (M x N matrix), flux axis (M x N matrix)
             - IFUSpectrum: spectral axis (M-length vector), flux axis (M x N x K matrix)
-            - TODO: BinnedSpectrum (final name(s) tbd):
-                - energy (M × 2 matrix), flux (N-length vector)
-                - others?
+            - BinnedSpectrum: spectral axis (M x 2 matrix of bin edges), flux axis (M-length vector)
 
             See the documentation for each spectrum type for more.
             """))
@@ -68,6 +66,19 @@ mutable struct Spectrum{S<:Number, F<:Number, M, N, A<:AbstractArray{S, M}, B<:A
             throw(ArgumentError("Spectral axis must be strictly increasing or decreasing."))
         end
 
+        # BinnedSpectrum invariants: m × 2 bin edges, each row strictly ordered in the axis direction
+        if M == 2 && N == 1
+            if size(s, 2) != 2
+                throw(ArgumentError("BinnedSpectrum spectral axis must be an m × 2 matrix of bin edges; got size $(size(s))"))
+            end
+            rows = eachrow(s)
+            ordered = all(issorted, spec_ax) ? all(r -> first(r) < last(r), rows) :
+                all(r -> first(r) > last(r), rows)
+            if !ordered
+                throw(ArgumentError("BinnedSpectrum bin edges must be strictly ordered within each row, in the direction of the spectral axis."))
+            end
+        end
+
         return new{S, F, M, N, A, B}(s, f, meta)
     end
 end
@@ -208,9 +219,9 @@ Unitful.unit(spec::AbstractSpectrum) = unit(eltype(spectral_axis(spec))), unit(e
 
 # Spectrum types and basic arithmetic
 include("spectrum_single.jl")
+include("spectrum_binned.jl")
 include("spectrum_echelle.jl")
 include("spectrum_ifu.jl")
-#include("spectrum_binned.jl")
 
 """
     spectrum(spectral_axis, flux_axis, [meta])
@@ -290,13 +301,25 @@ function spectrum(spectral_axis::AbstractMatrix{<:Real}, flux_axis::AbstractMatr
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
+function spectrum(spectral_axis::AbstractMatrix{<:Real}, flux_axis::AbstractVector{<:Real}; kwds...)
+    Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
+end
+
+# Recognized spectral-axis dimensions: wavelength, frequency, photon energy.
+const SPECTRAL_DIMENSIONS = (u"𝐋", u"𝐓^-1", u"𝐋^2 * 𝐌 * 𝐓^-2")
+
 function spectrum(spectral_axis::AbstractVector{<:Unitful.Quantity}, flux_axis::AbstractVector{<:Unitful.Quantity}; kwds...)
-    @assert dimension(eltype(spectral_axis)) ∈ (u"𝐋", u"𝐋^2 * 𝐌 * 𝐓^-2") "spectral_axis not recognized as having dimensions of wavelength or energy."
+    @assert dimension(eltype(spectral_axis)) ∈ SPECTRAL_DIMENSIONS "spectral_axis not recognized as having dimensions of wavelength, frequency, or energy."
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
 function spectrum(spectral_axis::AbstractMatrix{<:Unitful.Quantity}, flux_axis::AbstractMatrix{<:Unitful.Quantity}; kwds...)
-    @assert dimension(eltype(spectral_axis)) ∈ (u"𝐋", u"𝐋^2 * 𝐌 * 𝐓^-2") "spectral_axis not recognized as having dimensions of wavelength or energy."
+    @assert dimension(eltype(spectral_axis)) ∈ SPECTRAL_DIMENSIONS "spectral_axis not recognized as having dimensions of wavelength, frequency, or energy."
+    Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
+end
+
+function spectrum(spectral_axis::AbstractMatrix{<:Unitful.Quantity}, flux_axis::AbstractVector{<:Unitful.Quantity}; kwds...)
+    @assert dimension(eltype(spectral_axis)) ∈ SPECTRAL_DIMENSIONS "spectral_axis not recognized as having dimensions of wavelength, frequency, or energy."
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
diff --git a/src/spectrum_binned.jl b/src/spectrum_binned.jl
new file mode 100644
index 0000000..17978ea
--- /dev/null
+++ b/src/spectrum_binned.jl
@@ -0,0 +1,31 @@
+"""
+    BinnedSpectrum <: AbstractSpectrum
+
+An instance of [`Spectrum`](@ref) where the spectral axis is an ``m \\times 2`` matrix of bin edges, and the flux axis is a 1D array.
+"""
+const BinnedSpectrum = Spectrum{S, F, 2, 1} where {S, F}
+
+Base.getindex(spec::BinnedSpectrum, i::Integer) = spec[i:i]
+
+function Base.getindex(spec::BinnedSpectrum, inds)
+    return Spectrum(spectral_axis(spec)[inds, :], flux_axis(spec)[inds], meta(spec))
+end
+
+Base.firstindex(spec::BinnedSpectrum) = firstindex(flux_axis(spec))
+Base.lastindex(spec::BinnedSpectrum) = lastindex(flux_axis(spec))
+
+function Base.show(io::IO, spec::BinnedSpectrum)
+    w = spectral_axis(spec)
+    f = flux_axis(spec)
+    println(io, "BinnedSpectrum($(eltype(w)), $(eltype(f)))")
+    if isempty(f)
+        println(io, "  spectral axis $(size(w)): empty")
+        println(io, "  flux axis $(size(f)): empty")
+    else
+        println(io, "  spectral axis $(size(w)): ", first(w), " .. ", last(w))
+        println(io, "  flux axis $(size(f)): ", first(f), " .. ", last(f))
+    end
+    print(io, "  meta: ", meta(spec))
+end
+
+# rebinning -> important in the future, make sure rmf gets rebinned correctly or it could have different size
diff --git a/src/transforms/redden.jl b/src/transforms/redden.jl
index a5cf031..b34f910 100644
--- a/src/transforms/redden.jl
+++ b/src/transforms/redden.jl
@@ -1,5 +1,10 @@
 import DustExtinction
 
+# Extinction laws evaluate at a single wavelength per sample, which bin edges cannot provide.
+const _BINNED_REDDEN_MSG = "reddening is not defined for binned spectra: applying an extinction law to bin edges requires a bin-integration convention"
+redden!(::BinnedSpectrum, Av; kwargs...) = throw(ArgumentError(_BINNED_REDDEN_MSG))
+deredden!(::BinnedSpectrum, Av; kwargs...) = throw(ArgumentError(_BINNED_REDDEN_MSG))
+
 """
     redden!(::AbstractSpectrum, Av; Rv = 3.1, law = DustExtinction.CCM89)
 
diff --git a/src/transforms/redshift.jl b/src/transforms/redshift.jl
index e2265fe..24629a2 100644
--- a/src/transforms/redshift.jl
+++ b/src/transforms/redshift.jl
@@ -11,12 +11,19 @@ function _doppler_factor(v; relativistic)
     return sqrt((1 + β) / (1 - β))
 end
 
+# Cosmological redshift
+function _redshift_factor(z)
+    z > -1 || throw(DomainError(z, "Cosmological redshift requires z > -1."))
+    return 1 + z
+end
+
 # Resolve the (1 + z) stretch factor to the scalar actually multiplied onto a spectral axis of element type `S`.
 _spectral_axis_factor(::Type{<:Real}, factor) = factor # Unitless: assume wavelength
 _spectral_axis_factor(::Type{<:Unitful.Length}, factor) = factor # Wavelength: stretches
+_spectral_axis_factor(::Type{<:Unitful.Frequency}, factor) = inv(factor) # Frequency: compresses
 _spectral_axis_factor(::Type{<:Unitful.Energy}, factor) = inv(factor) # Energy: compresses
 _spectral_axis_factor(::Type{S}, factor) where {S <: Unitful.Quantity} =
-    throw(ArgumentError("cannot shift a spectral axis with dimension $(dimension(S)); expected a wavelength or energy"))
+    throw(ArgumentError("cannot shift a spectral axis with dimension $(dimension(S)); expected a wavelength, frequency, or energy"))
 _spectral_axis_factor(spec::AbstractSpectrum, factor) = _spectral_axis_factor(eltype(spectral_axis(spec)), factor)
 
 # Rebuilding through the constructor lets the element type promote and re-runs the invariant checks.
@@ -35,7 +42,7 @@ Reassigns the spectral axis, so its element and array type cannot change.
 Use [`redshift`](@ref) for spectra with an integer spectral axis.
 """
 function redshift!(spec::AbstractSpectrum, z::Real)
-    spec.spectral_axis = spectral_axis(spec) .* _spectral_axis_factor(spec, 1 + z)
+    spec.spectral_axis = spectral_axis(spec) .* _spectral_axis_factor(spec, _redshift_factor(z))
     return spec
 end
 
@@ -52,6 +59,7 @@ The observed wavelength is related to the rest-frame wavelength by:
 
 where ``z`` is the cosmological redshift parameter. Only the spectral axis is transformed.
 Flux density values are not corrected for the stretching of the wavelength bins.
+For a [`BinnedSpectrum`](@ref), both bin edges are shifted and the bin values are carried unchanged.
 
 A wavelength axis is stretched by ``1 + z``. An energy axis (distinguished by its `Unitful` dimension)
 is compressed by the same factor. A unitless axis is assumed to be a wavelength.
@@ -59,7 +67,7 @@ is compressed by the same factor. A unitless axis is assumed to be a wavelength.
 # Arguments
 
 - `spec`: The input spectrum.
-- `z`: Redshift parameter. Positive values redshift (longer wavelengths),
+- `z`: Redshift parameter, must satisfy `z > -1`. Positive values redshift (longer wavelengths),
   negative values blueshift (shorter wavelengths).
 
 # Examples
@@ -75,7 +83,7 @@ true
 
 See also [`doppler_shift`](@ref) for velocity-based Doppler shifting.
 """
-redshift(spec::AbstractSpectrum, z::Real) = _scale_spectral_axis(spec, 1 + z)
+redshift(spec::AbstractSpectrum, z::Real) = _scale_spectral_axis(spec, _redshift_factor(z))
 
 """
     doppler_shift!(spec::AbstractSpectrum, v; relativistic=false)
@@ -108,6 +116,7 @@ Apply a Doppler shift to a spectrum, returning a new spectrum with shifted wavel
 ```
 
 Only the spectral axis is transformed; flux density values are preserved as-is.
+For a [`BinnedSpectrum`](@ref), both bin edges are shifted and the bin values are carried unchanged.
 
 # Arguments
 
diff --git a/src/transforms/spectral_axis.jl b/src/transforms/spectral_axis.jl
new file mode 100644
index 0000000..f24d00f
--- /dev/null
+++ b/src/transforms/spectral_axis.jl
@@ -0,0 +1,65 @@
+using UnitfulEquivalences: Spectral, SpectralDensity
+
+# Same-dimension targets are a plain unit conversion, keeping conversions idempotent.
+_spectral_value(spectral_unit, x) =
+    dimension(spectral_unit) == dimension(x) ? Unitful.uconvert(spectral_unit, x) :
+    Unitful.uconvert(spectral_unit, x, Spectral())
+
+_flux_value(flux_unit, x, at) =
+    dimension(flux_unit) == dimension(x) ? Unitful.uconvert(flux_unit, x) :
+    Unitful.uconvert(flux_unit, x, SpectralDensity(at))
+
+_convert_spectral_axis(spec::Spectrum, spectral_unit) = _spectral_value.(spectral_unit, spectral_axis(spec))
+
+_convert_flux_axis(spec::Spectrum, ::Nothing) = copy(flux_axis(spec))
+function _convert_flux_axis(spec::Spectrum{<:Unitful.Quantity}, flux_unit::Unitful.Units)
+    spec isa BinnedSpectrum && throw(ArgumentError(
+        "flux density conversion is not defined for binned spectra: bin values are integrated quantities, so converting them requires a bin-integration convention"))
+    return _flux_value.(flux_unit, flux_axis(spec), spectral_axis(spec))
+end
+
+function _convert_spectrum(spec::Spectrum{<:Unitful.Quantity}, spectral_unit::Unitful.Units, flux_unit::Union{Nothing, Unitful.Units})
+    dimension(spectral_unit) ∈ SPECTRAL_DIMENSIONS || throw(ArgumentError(
+        "spectral axis unit must be a wavelength, frequency, or energy unit; got $spectral_unit"))
+    return Spectrum(_convert_spectral_axis(spec, spectral_unit), _convert_flux_axis(spec, flux_unit), deepcopy(meta(spec)))
+end
+
+@doc raw"""
+    Unitful.uconvert(spectral_unit::Unitful.Units, spec::Spectrum)
+    Unitful.uconvert(units::Tuple{Unitful.Units, Unitful.Units}, spec::Spectrum)
+
+Convert the spectral axis of `spec` to `spectral_unit`, which may be any wavelength, frequency, or
+photon-energy unit. Coordinates are related by the photon equivalence ``E = hν = hc/λ``.
+A target of the same dimension as the current axis is a plain unit conversion.
+
+The conversion is applied elementwise, so a reciprocal conversion (e.g., wavelength -->
+frequency) yields an axis with the opposite direction of monotonicity.
+
+Given a single unit, flux values are copied unchanged. Beware that a flux density then
+remains "per" its original coordinate (e.g., a per-wavelength density against a frequency
+axis). Pass a `(spectral_unit, flux_unit)` tuple, i.e., the shape returned by `unit(spec)`, to
+also convert flux density values into the matching convention via [`UnitfulEquivalences.SpectralDensity`](https://github.com/sostock/UnitfulEquivalences.jl),
+which preserves integrals: ``F_ν = F_λ λ^2 / c``.
+
+# Examples
+
+```jldoctest
+julia> using Unitful, UnitfulAstro
+
+julia> spec = spectrum([1.0, 1.5, 2.0]u"μm", [1.0, 2.0, 3.0]u"W/m^2/μm");
+
+julia> νspec = uconvert((u"THz", u"Jy"), spec);
+
+julia> issorted(spectral_axis(νspec); rev = true)
+true
+
+julia> flux_axis(νspec)[1] ≈ uconvert(u"Jy", 1.0u"W/m^2/μm" * (1.0u"μm")^2 / Unitful.c0)
+true
+```
+"""
+Unitful.uconvert(spectral_unit::Unitful.Units, spec::Spectrum{<:Unitful.Quantity}) = _convert_spectrum(spec, spectral_unit, nothing)
+Unitful.uconvert(units::Tuple{Unitful.Units, Unitful.Units}, spec::Spectrum{<:Unitful.Quantity}) = _convert_spectrum(spec, units[1], units[2])
+
+const _UNITLESS_AXIS_MSG = "`uconvert` requires a spectral axis with Unitful units. This spectrum's axis is unitless. Attach units first, e.g., `spectrum(axis * u\"angstrom\", flux)`"
+Unitful.uconvert(::Unitful.Units, ::AbstractSpectrum) = throw(ArgumentError(_UNITLESS_AXIS_MSG))
+Unitful.uconvert(::Tuple{Unitful.Units, Unitful.Units}, ::AbstractSpectrum) = throw(ArgumentError(_UNITLESS_AXIS_MSG))
diff --git a/src/transforms/transforms.jl b/src/transforms/transforms.jl
index fcd974a..fc6c02b 100644
--- a/src/transforms/transforms.jl
+++ b/src/transforms/transforms.jl
@@ -1,5 +1,6 @@
 include("redden.jl")
 include("redshift.jl")
+include("spectral_axis.jl")
 
 # ## Broadening ops
 # include("kernels.jl")
diff --git a/test/runtests.jl b/test/runtests.jl
index 3958875..b10d7fe 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,7 +4,7 @@ import SpectrumBase
 const init_code = quote
     using SpectrumBase: SpectrumBase, AbstractSpectrum, spectrum, spectral_axis, flux_axis
     using Measurements: Measurements, ±
-    using Unitful: @u_str, unit, ustrip
+    using Unitful: @u_str, uconvert, unit, ustrip
     using UnitfulAstro
     import Random
 
diff --git a/test/spectrum.jl b/test/spectrum.jl
index c757bd3..303b24f 100644
--- a/test/spectrum.jl
+++ b/test/spectrum.jl
@@ -1,4 +1,4 @@
-using SpectrumBase: SingleSpectrum, EchelleSpectrum, IFUSpectrum
+using SpectrumBase: SingleSpectrum, BinnedSpectrum, EchelleSpectrum, IFUSpectrum
 
 const RNG = Random.seed!(8675309)
 
@@ -20,6 +20,9 @@ const RNG = Random.seed!(8675309)
     wave[3] = 10
     @test_throws ArgumentError spectrum(wave, flux)
     @test_throws ArgumentError spectrum(wave[begin:end-1], flux)
+
+    wave, flux = [1 2; 2 3; 3 4; 4 5; 5 6], 1:5
+    @test spectrum(wave, flux) isa BinnedSpectrum
 end
 
 @testset "Spectrum - Single" begin
@@ -62,6 +65,51 @@ end
     @test spec.name == "test spectrum"
 end
 
+@testset "Spectrum - Binned" begin
+    wave = [1.0 2.0; 2.0 3.0; 3.0 4.0; 4.0 5.0]
+    flux = [10.0, 20.0, 30.0, 40.0]
+
+    spec = spectrum(wave, flux, name = "test binned spectrum")
+    spec_i = spec[2]
+    spec_I = spec[2:3]
+
+    @test spec isa BinnedSpectrum
+    @test spec_i isa BinnedSpectrum
+    @test spec_I isa BinnedSpectrum
+    @test spectral_axis(spec_i) == [2.0 3.0]
+    @test flux_axis(spec_i) == [20.0]
+    @test spectral_axis(spec_I) == [2.0 3.0; 3.0 4.0]
+    @test flux_axis(spec_I) == [20.0, 30.0]
+    @test firstindex(spec) == 1
+    @test lastindex(spec) == 4
+    @test spec[firstindex(spec)] == spec[1]
+    @test spec[lastindex(spec)] == spec[4]
+
+    # Any integer index selects a single bin
+    @test spec[Int32(2)] == spec[2] == spec[2:2]
+
+    # m × 2 shape and row-order invariants
+    @test_throws ArgumentError spectrum([1.0 2.0 3.0; 2.0 3.0 4.0], [1.0, 2.0]) # not two columns
+    @test_throws ArgumentError spectrum(reshape([1.0, 2.0], 2, 1), [1.0, 2.0]) # not two columns
+    @test_throws ArgumentError spectrum([2.0 1.0; 3.0 2.0], [1.0, 2.0]) # edges swapped within rows
+    @test_throws ArgumentError spectrum([1.0 2.0; 3.0 2.5], [1.0, 2.0]) # inconsistent row order
+    @test_throws ArgumentError spectrum([1.0 1.0; 2.0 2.0], [1.0, 2.0]) # zero-width bins
+    @test spectrum([4.0 3.0; 3.0 2.0], [1.0, 2.0]) isa BinnedSpectrum # descending edges are valid
+
+    # Empty selections construct and display
+    empty_spec = spec[falses(4)]
+    @test length(empty_spec) == 0
+    @test occursin("empty", sprint(show, empty_spec))
+
+    expected = """
+    BinnedSpectrum(Float64, Float64)
+      spectral axis (4, 2): 1.0 .. 5.0
+      flux axis (4,): 10.0 .. 40.0
+      meta: Dict{Symbol, Any}(:name => "test binned spectrum")"""
+    @test sprint(show, spec) == expected
+    @test spec.name == "test binned spectrum"
+end
+
 @testset "Spectrum - Echelle" begin
     n_wavs = 1000
     n_orders = 3
@@ -321,3 +369,18 @@ end
         @test flux_axis(s) ≈ flux_axis(spec) ./ 10
     end
 end
+
+@testset "Unitful Spectrum - Binned" begin
+    edges = [1.0 2.0; 2.0 3.0; 3.0 4.0]u"keV"
+    flux = [1.0, 2.0, 3.0]u"Jy"
+
+    spec = spectrum(edges, flux, name = "test xray spectrum")
+
+    @test spec isa BinnedSpectrum
+    @test unit(spec) == (u"keV", u"Jy")
+    @test spectral_axis(spec) == edges
+    @test flux_axis(spec) == flux
+    @test spec.name == "test xray spectrum"
+
+    @test_throws AssertionError spectrum([1.0 2.0; 2.0 3.0]u"kg", flux)
+end
diff --git a/test/transforms/redden.jl b/test/transforms/redden.jl
index 88cd372..71fc2b6 100644
--- a/test/transforms/redden.jl
+++ b/test/transforms/redden.jl
@@ -67,3 +67,11 @@ end
         @test flux_axis(dereddened) ≈ flux_axis(spec)
     end
 end
+
+@testset "Binned spectra are rejected" begin
+    binned = spectrum([4000.0 4010.0; 4010.0 4020.0], [1.0, 2.0])
+    @test_throws ArgumentError redden(binned, 0.3)
+    @test_throws ArgumentError redden!(binned, 0.3)
+    @test_throws ArgumentError deredden(binned, 0.3)
+    @test_throws ArgumentError deredden!(binned, 0.3)
+end
diff --git a/test/transforms/redshift.jl b/test/transforms/redshift.jl
index f1c7e05..d4bedf8 100644
--- a/test/transforms/redshift.jl
+++ b/test/transforms/redshift.jl
@@ -1,4 +1,5 @@
 using SpectrumBase:
+    BinnedSpectrum,
     redshift,
     redshift!,
     doppler_shift,
@@ -81,6 +82,36 @@ end
         spec = SpectrumBase.Spectrum(axis, ones(10) * u"erg/s/cm^2/angstrom", Dict{Symbol, Any}())
         @test_throws ArgumentError redshift(spec, 0.5)
     end
+
+    @testset "Unphysical redshift rejected" begin
+        spec = mock_spectrum_redshift()
+        @test_throws DomainError redshift(spec, -1.0) # Exactly z = -1
+        @test_throws DomainError redshift(spec, -1.5)
+        @test_throws DomainError redshift!(spec, -1.0)
+        # Blueshifts above the bound remain valid
+        @test redshift(spec, -0.5) isa typeof(spec)
+    end
+
+    @testset "Binned spectra" begin
+        edges = [1.0 2.0; 2.0 3.0; 3.0 4.0]
+        flux = [10.0, 20.0, 30.0]
+        binned = spectrum(edges, flux)
+
+        shifted = @inferred redshift(binned, 0.5)
+        @test shifted isa BinnedSpectrum
+        @test spectral_axis(shifted) ≈ edges .* 1.5
+        @test flux_axis(shifted) == flux
+        @test spectral_axis(binned) == edges
+
+        @test redshift!(binned, 0.5) === binned
+        @test spectral_axis(binned) ≈ edges .* 1.5
+
+        # Energy axis compresses
+        kev = spectrum(edges * u"keV", flux * u"Jy")
+        @test spectral_axis(@inferred redshift(kev, 0.5)) ≈ edges * u"keV" ./ 1.5
+        @test spectral_axis(doppler_shift(kev, 1000.0)) ≈ edges * u"keV" ./ (1 + 1000 / C_KMPS)
+        @test flux_axis(redshift(kev, 0.5)) == flux_axis(kev)
+    end
 end
 
 @testset "doppler_shift" begin
diff --git a/test/transforms/spectral_axis.jl b/test/transforms/spectral_axis.jl
new file mode 100644
index 0000000..2d231b0
--- /dev/null
+++ b/test/transforms/spectral_axis.jl
@@ -0,0 +1,102 @@
+using SpectrumBase: redshift
+
+const H_JS = 6.62607015e-34u"J*s" # Planck constant, exact by SI definition
+
+@testset "Wavelength and frequency conversion" begin
+    wavelength = [1.0, 1.5, 2.0]u"μm"
+    flux = [1.0, 2.0, 3.0]u"Jy"
+    spec = spectrum(wavelength, flux, name = "example")
+    c = C_MPS * u"m/s"
+
+    frequency_spec = @inferred uconvert(u"THz", spec)
+    expected_frequency = uconvert.(u"THz", c ./ wavelength)
+
+    # The conversion is elementwise, so an ascending wavelength axis becomes descending
+    @test issorted(spectral_axis(frequency_spec); rev = true)
+    @test spectral_axis(frequency_spec) ≈ expected_frequency
+    @test flux_axis(frequency_spec) == flux
+    @test frequency_spec.name == "example"
+    @test spectral_axis(spec) == wavelength
+
+    # Round trip restores the original axis and ordering
+    wavelength_spec = @inferred uconvert(u"μm", frequency_spec)
+    @test spectral_axis(wavelength_spec) ≈ wavelength
+    @test flux_axis(wavelength_spec) == flux
+
+    shifted = redshift(frequency_spec, 1.0)
+    @test spectral_axis(shifted) ≈ expected_frequency ./ 2
+
+    # Same-dimension targets are plain unit conversions, so conversions are idempotent
+    @test spectral_axis(uconvert(u"GHz", frequency_spec)) ≈
+        uconvert.(u"GHz", spectral_axis(frequency_spec))
+    @test spectral_axis(uconvert(u"nm", spec)) ≈ uconvert.(u"nm", wavelength)
+
+    # A `unit(spec)` tuple converts one spectrum to another's units
+    matched = uconvert(unit(frequency_spec), spec)
+    @test spectral_axis(matched) == spectral_axis(frequency_spec)
+    @test flux_axis(matched) == flux_axis(frequency_spec)
+
+    @test_throws MethodError uconvert(1u"Hz", spec) # a quantity, not a unit
+    @test_throws ArgumentError uconvert(u"kg", spec)
+    @test_throws ArgumentError uconvert(u"THz", spectrum([1.0, 2.0], [3.0, 4.0]))
+end
+
+@testset "Energy axis conversion" begin
+    energies = [1.0, 2.0, 5.0]u"keV"
+    flux = [1.0, 2.0, 3.0]u"Jy"
+    xray = spectrum(energies, flux)
+    c = C_MPS * u"m/s"
+
+    λspec = @inferred uconvert(u"angstrom", xray)
+    @test spectral_axis(λspec) ≈ uconvert.(u"angstrom", H_JS * c ./ energies)
+    @test flux_axis(λspec) == flux
+
+    # E → ν is linear, so the axis direction is preserved
+    νspec = @inferred uconvert(u"Hz", xray)
+    @test spectral_axis(νspec) ≈ uconvert.(u"Hz", energies ./ H_JS)
+    @test flux_axis(νspec) == flux
+
+    back = @inferred uconvert(u"keV", λspec)
+    @test spectral_axis(back) ≈ energies
+    @test flux_axis(back) == flux
+end
+
+@testset "Flux density conversion" begin
+    wavelength = [1.0, 1.5, 2.0]u"μm"
+    F_λ = [1.0, 2.0, 3.0]u"W/m^2/μm"
+    spec = spectrum(wavelength, F_λ)
+    c = C_MPS * u"m/s"
+
+    converted = @inferred uconvert((u"THz", u"Jy"), spec)
+    expected_F_ν = uconvert.(u"Jy", F_λ .* wavelength .^ 2 ./ c)
+    @test spectral_axis(converted) ≈ uconvert.(u"THz", c ./ wavelength)
+    @test flux_axis(converted) ≈ expected_F_ν
+
+    # Round trip back to F_λ
+    back = uconvert((u"μm", u"W/m^2/μm"), converted)
+    @test spectral_axis(back) ≈ wavelength
+    @test flux_axis(back) ≈ F_λ
+
+    # Per-energy density, the X-ray convention
+    F_E_spec = @inferred uconvert((u"keV", u"W/m^2/keV"), spec)
+    @test flux_axis(F_E_spec) ≈ uconvert.(u"W/m^2/keV", expected_F_ν ./ H_JS)
+
+    # A same-dimension flux target is a plain unit conversion
+    same = uconvert((u"THz", u"W/m^2/μm"), spec)
+    @test flux_axis(same) == F_λ
+end
+
+@testset "Binned spectrum axis conversion" begin
+    edges = [1.0 2.0; 2.0 3.0; 3.0 4.0]u"μm"
+    flux = [10.0, 20.0, 30.0]u"Jy"
+    binned = spectrum(edges, flux)
+    c = C_MPS * u"m/s"
+
+    converted = uconvert(u"THz", binned)
+    @test spectral_axis(converted) ≈ uconvert.(u"THz", c ./ edges)
+    @test flux_axis(converted) == flux
+    @test all(col -> issorted(col; rev = true), eachcol(spectral_axis(converted)))
+
+    # Bin values are integrated quantities; flux density conversion is undefined for them
+    @test_throws ArgumentError uconvert((u"THz", u"W/m^2/μm"), binned)
+end