diff --git a/docs/Project.toml b/docs/Project.toml
index 2a6dfef..737407d 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,8 @@
 [deps]
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
diff --git a/docs/make.jl b/docs/make.jl
index 459e71e..95c09cc 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -33,4 +33,4 @@ deploydocs(;
     devbranch = "main",
     push_preview = true,
     versions = ["stable" => "v^", "v#.#"], # Restrict to minor releases
-)
+)
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
index 176afb9..a4bd36d 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -25,9 +25,9 @@ Here is a quick demo of some of our features.
 ```jldoctest guide
 julia> using Spectra, FITSIO, Unitful, UnitfulAstro, Plots
 
-julia> fitsurl = "https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12";
+julia> # fitsurl = "https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12";
 
-julia> # download(fitsurl, "sdss.fits");
+julia> # f = FITS(HTTP.get(fitsurl).body)
 
 julia> f = FITS("sdss.fits")
 File: sdss.fits
@@ -65,4 +65,4 @@ TODO
 
 ## Contributing
 
-Please see [Contributing](@ref contrib) for information on contributing and extending Spectra.jl.
+Please see [Contributing](@ref contrib) for information on contributing and extending Spectra.jl.
\ No newline at end of file
diff --git a/docs/src/transforms.md b/docs/src/transforms.md
index 57c5c31..571792a 100644
--- a/docs/src/transforms.md
+++ b/docs/src/transforms.md
@@ -2,7 +2,7 @@
 
 ## Extinction
 
-By levaraging [DustExtinction.jl](https://github.com/juliaastro/dustextinction.jl) we can apply common reddening laws to our spectra.
+By leveraging [DustExtinction.jl](https://github.com/juliaastro/dustextinction.jl) we can apply common reddening laws to our spectra.
 
 ```jldoctest
 julia> using Spectra, Unitful, Measurements, Random
@@ -59,3 +59,11 @@ redden!
 deredden
 deredden!
 ```
+
+## 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:
+
+```@docs
+SpectrumResampler
+```
\ No newline at end of file
diff --git a/src/Spectra.jl b/src/Spectra.jl
index e1ad5f9..4c61508 100644
--- a/src/Spectra.jl
+++ b/src/Spectra.jl
@@ -2,14 +2,18 @@ module Spectra
 
 # Uniform API
 export AbstractSpectrum, Spectrum, spectrum, spectral_axis, flux_axis
-# spectra_single.jl, spectra_ifu.jl, spectra_echelle.jl, spectra_binned.jl
-export SingleSpectrum, IFUSpectrum, EchelleSpectrum, spectral_axis, flux_axis
-# utils.jl
+
+# AbstractSpectrum types
+export SingleSpectrum, IFUSpectrum, EchelleSpectrum
+
+# Transforms
+export SpectrumResampler, redden, redden!, deredden, deredden!
+
+# Utilities
 export blackbody #, line_flux, equivalent_width
-# fitting/fitting.jl
+
+# Fitting
 #export continuum, continuum!
-# transforms/redden.jl
-export redden, redden!, deredden, deredden!
 
 using RecipesBase: @recipe
 using Measurements: Measurements, Measurement
@@ -82,7 +86,7 @@ Return the spectral axis of `spec`.
 spectral_axis(spec::AbstractSpectrum) = spec.spectral_axis
 
 """
-    flux(spec::AbstractSpectrum)
+    flux_axis(spec::AbstractSpectrum)
 
 Return the flux axis of `spec`.
 """
@@ -294,8 +298,8 @@ function spectrum(spectral_axis::AbstractMatrix{<:Quantity}, flux_axis::Abstract
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
-# tools
 include("utils.jl")
+include("transforms/resampler.jl")
 include("transforms/transforms.jl")
 include("plotting.jl")
 #include("fitting/fitting.jl")
diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl
new file mode 100644
index 0000000..1de68ad
--- /dev/null
+++ b/src/transforms/resampler.jl
@@ -0,0 +1,88 @@
+"""
+    SpectrumResampler(spec::Spectrum, interp)
+
+Type representing the spectrum `spec` with interpolator `interp`.
+
+Interpolation methods from many packages can be used without issue. Below we show example usage of [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) and [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl).
+
+First, we set up an arbitrary spectrum and a linear interpolator from DataInterpolations.jl:
+
+```jldoctest resampling
+julia> using Spectra: SpectrumResampler, spectrum, spectral_axis, flux_axis
+
+julia> using DataInterpolations: LinearInterpolation, ExtrapolationType
+
+julia> spec = spectrum([20, 40, 120, 160, 200], [1, 3, 7, 6, 20]);
+
+julia> interp = LinearInterpolation(flux_axis(spec), spectral_axis(spec);
+           extrapolation = ExtrapolationType.Constant);
+```
+
+Now, we construct the `SpectrumResampler` and define the new wavelength grid that we want to resample the original spectrum to:
+
+```jldoctest resampling
+julia> resampler = SpectrumResampler(spec, interp);
+
+julia> wave_sampled = [10, 50, 90, 130, 140, 170, 210, 220, 230];
+```
+
+To perform the resampling, you call the resampler with the desired wavelength grid.
+
+```jldoctest resampling
+julia> result = resampler(wave_sampled);
+```
+
+The resampled wavelength and flux can be obtained with the `wave` and `flux` methods.
+
+```jldoctest resampling
+julia> result isa Spectrum
+true
+
+julia> spectral_axis(result) == wave_sampled
+true
+
+julia> flux_axis(result) == interp(wave_sampled)
+true
+```
+
+Use of [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) follows the same general procedure, but using a different `interp`:
+
+```jldoctest resampling
+julia> using Interpolations: linear_interpolation, Flat
+
+julia> interp = linear_interpolation(spectral_axis(spec), flux_axis(spec); extrapolation_bc = Flat());
+
+julia> resampler = SpectrumResampler(spec, interp);
+
+julia> result = resampler(wave_sampled);
+
+julia> result isa Spectrum
+true
+
+julia> spectral_axis(result) == wave_sampled
+true
+
+julia> flux_axis(result) == interp(wave_sampled)
+true
+```
+"""
+struct SpectrumResampler{A <: Spectrum, B}
+    spectrum::A
+    interp::B
+end
+
+spectrum(s::SpectrumResampler) = s.spectrum
+spectral_axis(s::SpectrumResampler) = (spectral_axis ∘ spectrum)(s)
+flux_axis(s::SpectrumResampler) = (flux_axis ∘ spectrum)(s)
+meta(s::SpectrumResampler) = (meta ∘ spectrum)(s)
+
+function (s::SpectrumResampler)(wave_sampled)
+    flux_resampled = (s.interp)(wave_sampled)
+    return Spectrum(wave_sampled, flux_resampled, meta(s))
+end
+
+function Base.show(io::IO, s::SpectrumResampler)
+    println(io, "SpectrumResampler(", (eltype ∘ spectral_axis)(s), ", ", (eltype ∘ flux_axis)(s), ")")
+    println(io, "  spec: ", (typeof ∘ spectrum)(s))
+    print(io, "  interpolator: ", typeof(s.interp))
+end
diff --git a/test/Project.toml b/test/Project.toml
index 99cdd30..512392c 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DustExtinction = "fb44c06c-c62f-5397-83f5-69249e0a3c8e"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
diff --git a/test/transforms/resample.jl b/test/transforms/resample.jl
index b5c219f..1a599f1 100644
--- a/test/transforms/resample.jl
+++ b/test/transforms/resample.jl
@@ -1,72 +1,99 @@
-#function mock_spectrum(n::Int = Int(1e3); use_units::Bool = false)
-#    wave = range(1e4, 3e4, length = n)
-#    sigma = 0.1 .* sin.(wave)
-#    T = 6700
-#    flux = @. 1e14 / (wave^5 * (exp(1 / (wave * T)) - 1)) ± sigma
-#    if use_units
-#        wave *= u"angstrom"
-#        flux *= u"erg/s/cm^2/angstrom"
-#    end
-#    spectrum(wave, flux, name="Test Spectrum")
-#end
-#
-#@testset "resample" begin
-#    @testset "Resampling" begin
-#        spec = mock_spectrum()
-#        new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) ÷ 2.4))
-#        res_spec = resample(spec, new_wave)
-#
-#        @test res_spec.wave == new_wave
-#        @test length(res_spec.flux) == length(new_wave)
-#
-#        resample!(spec, new_wave)
-#        @test res_spec.wave == spec.wave
-#        @test res_spec.flux == spec.flux
-#
-#        # Unitful
-#        spec = mock_spectrum(use_units=true)
-#        new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) ÷ 2.4))
-#        @assert unit(eltype(new_wave)) == unit(eltype(spec.wave))
-#        res_spec = resample(spec, new_wave)
-#
-#        @test res_spec.wave == new_wave
-#        @test length(res_spec.flux) == length(new_wave)
-#
-#        resample!(spec, new_wave)
-#        @test res_spec.wave == spec.wave
-#        @test res_spec.flux == spec.flux
-#
-#        # Test resampling to another Spectrum
-#        spec1 = mock_spectrum(Integer(1e4))
-#        spec2 = mock_spectrum(Integer(1e3))
-#        res_spec = resample(spec1, spec2)
-#        @test res_spec.wave == spec2.wave
-#
-#        resample!(spec1, spec2)
-#        @test spec1.wave == spec2.wave
-#        @test spec1.flux == res_spec.flux
-#
-#        # Unitful
-#        spec1 = mock_spectrum(Integer(1e4), use_units=true)
-#        spec2 = mock_spectrum(Integer(1e3), use_units=true)
-#        res_spec = resample(spec1, spec2)
-#
-#        @test res_spec.wave == spec2.wave
-#        @test unit(eltype(spec1.flux)) == unit(eltype(spec2.flux)) == unit(eltype(res_spec.flux))
-#
-#        # Test when spectra2 has different units
-#        spec1 = mock_spectrum(Integer(1e4), use_units=true)
-#        spec2 = mock_spectrum(Integer(1e3), use_units=true)
-#        spec1.wave = uconvert.(u"cm", spec1.wave)
-#        resample!(spec1, spec2)
-#        @test ustrip.(spec1.wave) ≈ ustrip.(unit(eltype(spec1.wave)), spec2.wave)
-#
-#        # address bug when doing the same thing but affecting spec2
-#        spec1 = mock_spectrum(Integer(1e4), use_units=true)
-#        spec2 = mock_spectrum(Integer(1e3), use_units=true)
-#        spec2.wave = uconvert.(u"cm", spec2.wave)
-#        resample!(spec1, spec2)
-#        @test ustrip.(spec1.wave) ≈ ustrip.(unit(eltype(spec1.wave)), spec2.wave)
-#
-#    end
-#end
+using Spectra: Spectra, AbstractSpectrum, SpectrumResampler, spectrum, spectral_axis, flux_axis
+using DataInterpolations: LinearInterpolation, ExtrapolationType
+using Unitful: @u_str, uconvert
+using UnitfulAstro
+using Measurements: ±
+
+# TODO: See if it makes sense to have an exported API for resample/resample! even though we are using external interpolators.
+
+function resample(spec, new_wave)
+    interp = LinearInterpolation(flux_axis(spec), spectral_axis(spec); extrapolation = ExtrapolationType.Constant)
+    resampler = SpectrumResampler(spec, interp)
+    return resampler(new_wave)
+end
+
+function resample(spec1::AbstractSpectrum, spec2::AbstractSpectrum)
+    new_wave = spectral_axis(spec2)
+    return resample(spec1, new_wave)
+end
+
+function mock_spectrum(n::Int = Int(1e3); use_units::Bool = false)
+    wave = range(1e4, 3e4, length = n)
+    sigma = 0.1 .* sin.(wave)
+    T = 6700
+    flux = @. 1e14 / (wave^5 * (exp(1 / (wave * T)) - 1)) ± sigma
+    if use_units
+        wave *= u"angstrom"
+        flux *= u"erg/s/cm^2/angstrom"
+    end
+    spectrum(wave, flux; name = "Test Spectrum")
+end
+
+@testset "Resampler" begin
+    spec = mock_spectrum()
+    s, f = spectral_axis(spec), flux_axis(spec)
+    interp = LinearInterpolation(f, s; extrapolation = ExtrapolationType.Constant)
+    resampler = SpectrumResampler(spec, interp)
+    expected = """
+    SpectrumResampler(Float64, Measurements.Measurement{Float64})
+      spec: Spectra.SingleSpectrum{Float64, Measurements.Measurement{Float64}}
+      interpolator: DataInterpolations.LinearInterpolation{Vector{Measurements.Measurement{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Measurements.Measurement{Float64}}, Vector{Measurements.Measurement{Float64}}, Measurements.Measurement{Float64}}"""
+
+    @test sprint(show, resampler) == expected
+    @test spectral_axis(resampler) == s
+    @test flux_axis(resampler) == f
+end
+
+@testset "Resampling" begin
+    spec = mock_spectrum()
+    s, f = spectral_axis(spec), flux_axis(spec)
+    new_wave = range(minimum(s), maximum(s); length = Integer(length(s) ÷ 2.4))
+    res_spec = resample(spec, new_wave)
+    expected = """
+    SingleSpectrum(Float64, Measurements.Measurement{Float64})
+      spectral axis (416,): 10000.0 .. 30000.0
+      flux axis (416,): 67.0 ± 0.031 .. 0.827 ± 0.08
+      meta: Dict{Symbol, Any}(:name => "Test Spectrum")"""
+
+    @test sprint(show, res_spec) == expected
+    @test spectral_axis(res_spec) == new_wave
+    @test length(flux_axis(res_spec)) == length(new_wave)
+
+    # Unitful
+    spec = mock_spectrum(; use_units = true)
+    s, f = spectral_axis(spec), flux_axis(spec)
+    new_wave = range(minimum(s), maximum(s); length = Integer(length(s) ÷ 2.4))
+    @assert unit(eltype(new_wave)) == unit(eltype(s))
+    res_spec = resample(spec, new_wave)
+
+    @test spectral_axis(res_spec) == new_wave
+    @test length(flux_axis(res_spec)) == length(new_wave)
+
+    # Test resampling to another Spectrum
+    spec1 = mock_spectrum(Integer(1e4))
+    spec2 = mock_spectrum(Integer(1e3))
+    res_spec = resample(spec1, spec2)
+    @test spectral_axis(res_spec) == spectral_axis(spec2)
+
+    # Unitful
+    spec1 = mock_spectrum(Integer(1e4); use_units = true)
+    spec2 = mock_spectrum(Integer(1e3); use_units = true)
+    res_spec = resample(spec1, spec2)
+
+    @test spectral_axis(res_spec) == spectral_axis(spec2)
+    @test unit(eltype(flux_axis(spec1))) == unit(eltype(flux_axis(spec2))) == unit(eltype(flux_axis(res_spec)))
+
+    # Test when spectra2 has different units
+    #spec1 = mock_spectrum(Integer(1e4), use_units=true)
+    #spec2 = mock_spectrum(Integer(1e3), use_units=true)
+    #spec1.wave = uconvert.(u"cm", spec1.wave)
+    #resample!(spec1, spec2)
+    #@test ustrip.(spectral_axis(spec1)) ≈ ustrip.(unit(eltype(spectral_axis(spec1))), spectral_axis(spec2))
+
+    # address bug when doing the same thing but affecting spec2
+    #spec1 = mock_spectrum(Integer(1e4), use_units=true)
+    #spec2 = mock_spectrum(Integer(1e3), use_units=true)
+    #spec2.wave = uconvert.(u"cm", spec2.wave)
+    #resample!(spec1, spec2)
+    #@test ustrip.(spectral_axis(spec1)) ≈ ustrip.(unit(eltype(spectral_axis(spec1))), spectral_axis(spec2))
+end