icweaver · Sep 17, 2026 · Aug 7, 2026 · Aug 11, 2026 · Aug 13, 2026 · Aug 17, 2026
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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,14 +59,15 @@ 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.
 
 # 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
@@ -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
@@ -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
@@ -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