diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9c59465b1..1363bc0c4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## Unreleased
 
 - Change `SingleSiteOperator` with the more general `MultiSiteOperator`. ([#324])
+- Make `spectrum` and `correlation` functions align with `Python QuTiP`, introduce spectrum solver `PseudoInverse`, remove spectrum solver `FFTCorrelation`, and introduce `spectrum_correlation_fft`. ([#330])
 
 ## [v0.22.0] (2024-11-20)
 
@@ -36,3 +37,4 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 [#315]: https://github.com/qutip/QuantumToolbox.jl/issues/315
 [#318]: https://github.com/qutip/QuantumToolbox.jl/issues/318
 [#324]: https://github.com/qutip/QuantumToolbox.jl/issues/324
+[#330]: https://github.com/qutip/QuantumToolbox.jl/issues/330
diff --git a/benchmarks/correlations_and_spectrum.jl b/benchmarks/correlations_and_spectrum.jl
index a5eab65a8..04f63997f 100644
--- a/benchmarks/correlations_and_spectrum.jl
+++ b/benchmarks/correlations_and_spectrum.jl
@@ -1,3 +1,9 @@
+function _calculate_fft_spectrum(H, tlist, c_ops, A, B)
+    corr = correlation_2op_1t(H, nothing, tlist, c_ops, A, B; progress_bar = Val(false))
+    ωlist, spec = spectrum_correlation_fft(tlist, corr)
+    return nothing
+end
+
 function benchmark_correlations_and_spectrum!(SUITE)
     N = 15
     ω = 1
@@ -9,11 +15,18 @@ function benchmark_correlations_and_spectrum!(SUITE)
     c_ops = [sqrt(γ * (nth + 1)) * a, sqrt(γ * nth) * a']
 
     ω_l = range(0, 3, length = 1000)
+    t_l = range(0, 333 * π, length = 1000)
+
+    PI_solver = PseudoInverse()
 
     SUITE["Correlations and Spectrum"]["FFT Correlation"] =
-        @benchmarkable spectrum($H, $ω_l, $(a'), $a, $c_ops, solver = FFTCorrelation(), progress_bar = false)
+        @benchmarkable _calculate_fft_spectrum($H, $t_l, $c_ops, $(a'), $a)
+
+    SUITE["Correlations and Spectrum"]["Spectrum"]["Exponential Series"] =
+        @benchmarkable spectrum($H, $ω_l, $c_ops, $(a'), $a)
 
-    SUITE["Correlations and Spectrum"]["Exponential Series"] = @benchmarkable spectrum($H, $ω_l, $(a'), $a, $c_ops)
+    SUITE["Correlations and Spectrum"]["Spectrum"]["Pseudo Inverse"] =
+        @benchmarkable spectrum($H, $ω_l, $c_ops, $(a'), $a, solver = $PI_solver)
 
     return nothing
 end
diff --git a/docs/make.jl b/docs/make.jl
index 3aca7289f..5512e8621 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -50,8 +50,7 @@ const PAGES = [
             "Solving Problems with Time-dependent Hamiltonians" => "users_guide/time_evolution/time_dependent.md",
         ],
         "Solving for Steady-State Solutions" => "users_guide/steadystate.md",
-        "Symmetries" => [],
-        "Two-time correlation functions" => [],
+        "Two-time correlation functions" => "users_guide/two_time_corr_func.md",
         "Extensions" => [
             "users_guide/extensions/cuda.md",
         ],
diff --git a/docs/src/resources/api.md b/docs/src/resources/api.md
index d19c4b86c..d20fc6805 100644
--- a/docs/src/resources/api.md
+++ b/docs/src/resources/api.md
@@ -230,9 +230,13 @@ lr_mesolve
 
 ```@docs
 correlation_3op_2t
+correlation_3op_1t
 correlation_2op_2t
 correlation_2op_1t
+spectrum_correlation_fft
 spectrum
+ExponentialSeries
+PseudoInverse
 ```
 
 ## [Metrics](@id doc-API:Metrics)
diff --git a/docs/src/resources/bibliography.bib b/docs/src/resources/bibliography.bib
index d83932eed..3cddf83b5 100644
--- a/docs/src/resources/bibliography.bib
+++ b/docs/src/resources/bibliography.bib
@@ -1,3 +1,13 @@
+@book{Gardiner-Zoller2004,
+  title = {Quantum Noise},
+  ISBN = {9783540223016},
+  url = {https://link.springer.com/book/9783540223016},
+  publisher = {Springer Berlin, Heidelberg},
+  author = {Gardiner, Crispin and Zoller, Peter},
+  year = {2004},
+  month = aug 
+}
+
 @book{Nielsen-Chuang2011,
   title = {Quantum Computation and Quantum Information: 10th Anniversary Edition},
   ISBN = {9780511976667},
diff --git a/docs/src/users_guide/two_time_corr_func.md b/docs/src/users_guide/two_time_corr_func.md
new file mode 100644
index 000000000..b3ad59730
--- /dev/null
+++ b/docs/src/users_guide/two_time_corr_func.md
@@ -0,0 +1,266 @@
+# [Two-time Correlation Functions](@id doc:Two-time-Correlation-Functions)
+
+## Introduction
+
+With the `QuantumToolbox.jl` time-evolution function [`mesolve`](@ref), a state vector ([`Ket`](@ref)) or density matrix ([`Operator`](@ref)) can be evolved from an initial state at ``t_0`` to an arbitrary time ``t``, namely
+
+```math
+\hat{\rho}(t) = \mathcal{G}(t, t_0)\{\hat{\rho}(t_0)\},
+```
+where ``\mathcal{G}(t, t_0)\{\cdot\}`` is the propagator defined by the equation of motion. The resulting density matrix can then be used to evaluate the expectation values of arbitrary combinations of same-time operators.
+
+To calculate two-time correlation functions on the form ``\left\langle \hat{A}(t+\tau) \hat{B}(t) \right\rangle``, we can use the quantum regression theorem [see, e.g., [Gardiner-Zoller2004](@cite)] to write
+
+```math
+\left\langle \hat{A}(t+\tau) \hat{B}(t) \right\rangle = \textrm{Tr} \left[\hat{A} \mathcal{G}(t+\tau, t)\{\hat{B}\hat{\rho}(t)\} \right] = \textrm{Tr} \left[\hat{A} \mathcal{G}(t+\tau, t)\{\hat{B} \mathcal{G}(t, 0)\{\hat{\rho}(0)\}\} \right],
+```
+
+We therefore first calculate ``\hat{\rho}(t) = \mathcal{G}(t, 0)\{\hat{\rho}(0)\}`` using [`mesolve`](@ref) with ``\hat{\rho}(0)`` as initial state, and then again use [`mesolve`](@ref) to calculate ``\mathcal{G}(t+\tau, t)\{\hat{B}\hat{\rho}(t)\}`` using ``\hat{B}\hat{\rho}(t)`` as initial state.
+
+Note that if the initial state is the steady state, then ``\hat{\rho}(t) = \mathcal{G}(t, 0)\{\hat{\rho}_{\textrm{ss}}\} = \hat{\rho}_{\textrm{ss}}`` and
+
+```math
+\left\langle \hat{A}(t+\tau) \hat{B}(t) \right\rangle = \textrm{Tr} \left[\hat{A} \mathcal{G}(t+\tau, t)\{\hat{B}\hat{\rho}_{\textrm{ss}}\} \right] = \textrm{Tr} \left[\hat{A} \mathcal{G}(\tau, 0)\{\hat{B} \hat{\rho}_{\textrm{ss}}\} \right] = \left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle,
+```
+which is independent of ``t``, so that we only have one time coordinate ``\tau``.
+
+`QuantumToolbox.jl` provides a family of functions that assists in the process of calculating two-time correlation functions. The available functions and their usage is shown in the table below.
+
+| **Function call** | **Correlation function** |
+|:------------------|:-------------------------|
+| [`correlation_2op_2t`](@ref) | ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle`` or ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \right\rangle`` |
+| [`correlation_2op_1t`](@ref) | ``\left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle`` or ``\left\langle \hat{A}(0) \hat{B}(\tau) \right\rangle`` |
+| [`correlation_3op_1t`](@ref) | ``\left\langle \hat{A}(0) \hat{B}(\tau) \hat{C}(0) \right\rangle`` |
+| [`correlation_3op_2t`](@ref) | ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle`` |
+
+The most common used case is to calculate the two time correlation function ``\left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle``, which can be done by [`correlation_2op_1t`](@ref).
+
+```@setup correlation_and_spectrum
+using QuantumToolbox
+
+using CairoMakie
+CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
+```
+
+## Steadystate correlation function
+
+The following code demonstrates how to calculate the ``\langle \hat{x}(t) \hat{x}(0)\rangle`` correlation for a leaky cavity with three different relaxation rates ``\gamma``.
+
+```@example correlation_and_spectrum
+tlist = LinRange(0, 10, 200)
+a = destroy(10)
+x = a' + a
+H = a' * a
+
+# if the initial state is specified as `nothing`, the steady state will be calculated and used as the initial state.
+corr1 = correlation_2op_1t(H, nothing, tlist, [sqrt(0.5) * a], x, x)
+corr2 = correlation_2op_1t(H, nothing, tlist, [sqrt(1.0) * a], x, x)
+corr3 = correlation_2op_1t(H, nothing, tlist, [sqrt(2.0) * a], x, x)
+
+# plot by CairoMakie.jl
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], xlabel = L"Time $t$", ylabel = L"\langle \hat{x}(t) \hat{x}(0) \rangle")
+lines!(ax, tlist, real(corr1), label = L"\gamma = 0.5", linestyle = :solid)
+lines!(ax, tlist, real(corr2), label = L"\gamma = 1.0", linestyle = :dash)
+lines!(ax, tlist, real(corr3), label = L"\gamma = 2.0", linestyle = :dashdot)
+
+axislegend(ax, position = :rt)
+
+fig
+```
+
+## Emission spectrum
+
+Given a correlation function ``\left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle``, we can define the corresponding power spectrum as
+
+```math
+S(\omega) = \int_{-\infty}^\infty \left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle e^{-i \omega \tau} d \tau
+```
+
+In `QuantumToolbox.jl`, we can calculate ``S(\omega)`` using either [`spectrum`](@ref), which provides several solvers to perform the Fourier transform semi-analytically, or we can use the function [`spectrum_correlation_fft`](@ref) to numerically calculate the fast Fourier transform (FFT) of a given correlation data.
+
+The following example demonstrates how these methods can be used to obtain the emission (``\hat{A} = \hat{a}^\dagger`` and ``\hat{B} = \hat{a}``) power spectrum.
+
+```@example correlation_and_spectrum
+N = 4             # number of cavity fock states
+ωc = 1.0 * 2 * π  # cavity frequency
+ωa = 1.0 * 2 * π  # atom frequency
+g  = 0.1 * 2 * π  # coupling strength
+κ  = 0.75         # cavity dissipation rate
+γ  = 0.25         # atom dissipation rate
+
+# Jaynes-Cummings Hamiltonian
+a  = tensor(destroy(N), qeye(2))
+sm = tensor(qeye(N), destroy(2))
+H = ωc * a' * a + ωa * sm' * sm + g * (a' * sm + a * sm')
+
+# collapse operators
+n_th = 0.25
+c_ops = [
+    sqrt(κ * (1 + n_th)) * a,
+    sqrt(κ *      n_th)  * a',
+    sqrt(γ)              * sm,
+];
+
+# calculate the correlation function using mesolve, and then FFT to obtain the spectrum. 
+# Here we need to make sure to evaluate the correlation function for a sufficient long time and 
+# sufficiently high sampling rate so that FFT captures all the features in the resulting spectrum.
+tlist = LinRange(0, 100, 5000)
+corr = correlation_2op_1t(H, nothing, tlist, c_ops, a', a; progress_bar = Val(false))
+ωlist1, spec1 = spectrum_correlation_fft(tlist, corr)
+
+# calculate the power spectrum using spectrum
+# using Exponential Series (default) method
+ωlist2 = LinRange(0.25, 1.75, 200) * 2 * π
+spec2 = spectrum(H, ωlist2, c_ops, a', a; solver = ExponentialSeries())
+
+# calculate the power spectrum using spectrum
+# using Pseudo-Inverse method
+spec3 = spectrum(H, ωlist2, c_ops, a', a; solver = PseudoInverse())
+
+# plot by CairoMakie.jl
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], title = "Vacuum Rabi splitting", xlabel = "Frequency", ylabel = "Emission power spectrum")
+lines!(ax, ωlist1 / (2 * π), spec1, label = "mesolve + FFT", linestyle = :solid)
+lines!(ax, ωlist2 / (2 * π), spec2, label = "Exponential Series", linestyle = :dash)
+lines!(ax, ωlist2 / (2 * π), spec3, label = "Pseudo-Inverse", linestyle = :dashdot)
+
+xlims!(ax, ωlist2[1] / (2 * π), ωlist2[end] / (2 * π))
+axislegend(ax, position = :rt)
+
+fig
+```
+
+## Non-steadystate correlation function
+
+More generally, we can also calculate correlation functions of the kind ``\left\langle \hat{A}(t_1 + t_2) \hat{B}(t_1) \right\rangle``, i.e., the correlation function of a system that is not in its steady state. In `QuantumToolbox.jl`, we can evaluate such correlation functions using the function [`correlation_2op_2t`](@ref). The default behavior of this function is to return a matrix with the correlations as a function of the two time coordinates (``t_1`` and ``t_2``).
+
+```@example correlation_and_spectrum
+t1_list = LinRange(0, 10.0, 200)
+t2_list = LinRange(0, 10.0, 200)
+
+N = 10
+a = destroy(N)
+x = a' + a
+H = a' * a
+
+c_ops = [sqrt(0.25) * a]
+
+α = 2.5
+ρ0 = coherent_dm(N, α)
+
+corr = correlation_2op_2t(H, ρ0, t1_list, t2_list, c_ops, x, x; progress_bar = Val(false))
+
+# plot by CairoMakie.jl
+fig = Figure(size = (500, 400))
+
+ax = Axis(fig[1, 1], title = L"\langle \hat{x}(t_1 + t_2) \hat{x}(t_1) \rangle", xlabel = L"Time $t_1$", ylabel = L"Time $t_2$")
+
+heatmap!(ax, t1_list, t2_list, real(corr))
+
+fig
+```
+
+### Example: first-order optical coherence function
+
+This example demonstrates how to calculate a correlation function on the form ``\left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle`` for a non-steady initial state. Consider an oscillator that is interacting with a thermal environment. If the oscillator initially is in a coherent state, it will gradually decay to a thermal (incoherent) state. The amount of coherence can be quantified using the first-order optical coherence function
+
+```math
+g^{(1)}(\tau) = \frac{\left\langle \hat{a}^\dagger(\tau) \hat{a}(0) \right\rangle}{\sqrt{\left\langle \hat{a}^\dagger(\tau) \hat{a}(\tau) \right\rangle \left\langle \hat{a}^\dagger(0) \hat{a}(0)\right\rangle}}.
+```
+For a coherent state ``\vert g^{(1)}(\tau) \vert = 1``, and for a completely incoherent (thermal) state ``g^{(1)}(\tau) = 0``. The following code calculates and plots ``g^{(1)}(\tau)`` as a function of ``\tau``:
+
+```@example correlation_and_spectrum
+τlist = LinRange(0, 10, 200)
+
+# Hamiltonian
+N = 15
+a = destroy(N)
+H = 2 * π * a' * a
+
+# collapse operator
+G1 = 0.75
+n_th = 2.00  # bath temperature in terms of excitation number
+c_ops = [
+    sqrt(G1 * (1 + n_th)) * a,
+    sqrt(G1 *      n_th)  * a'
+]
+
+# start with a coherent state of α = 2.0
+ρ0 = coherent_dm(N, 2.0)
+
+# first calculate the occupation number as a function of time
+n = mesolve(H, ρ0, τlist, c_ops, e_ops = [a' * a], progress_bar = Val(false)).expect[1,:]
+n0 = n[1] # occupation number at τ = 0
+
+# calculate the correlation function G1 and normalize with n to obtain g1
+g1 = correlation_2op_1t(H, ρ0, τlist, c_ops, a', a, progress_bar = Val(false))
+g1 = g1 ./ sqrt.(n .* n0)
+
+# plot by CairoMakie.jl
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], title = "Decay of a coherent state to an incoherent (thermal) state", xlabel = L"Time $\tau$")
+lines!(ax, τlist, real(g1), label = L"g^{(1)}(\tau)", linestyle = :solid)
+lines!(ax, τlist, real(n), label = L"n(\tau)", linestyle = :dash)
+
+axislegend(ax, position = :rt)
+
+fig
+```
+
+### Example: second-order optical coherence function
+
+The second-order optical coherence function, with time-delay ``\tau``, is defined as
+
+```math
+g^{(2)}(\tau) = \frac{\left\langle \hat{a}^\dagger(0) \hat{a}^\dagger(\tau) \hat{a}(\tau) \hat{a}(0) \right\rangle}{\left\langle \hat{a}^\dagger(0) \hat{a}(0) \right\rangle^2}.
+```
+
+For a coherent state ``g^{(2)}(\tau) = 1``, for a thermal state ``g^{(2)}(\tau = 0) = 2`` and it decreases as a function of time (bunched photons, they tend to appear together), and for a Fock state with ``n``-photons ``g^{(2)}(\tau = 0) = n(n-1)/n^2 < 1`` and it increases with time (anti-bunched photons, more likely to arrive separated in time).
+
+To calculate this type of correlation function with `QuantumToolbox.jl`, we can use [`correlation_3op_1t`](@ref), which computes a correlation function on the form ``\left\langle \hat{A}(0) \hat{B}(\tau) \hat{C}(0) \right\rangle`` (three operators and one delay-time vector). We first have to combine the central two operators into one single one as they are evaluated at the same time, e.g. here we do ``\hat{B}(\tau) = \hat{a}^\dagger(\tau) \hat{a}(\tau) = (\hat{a}^\dagger\hat{a})(\tau)``.
+
+The following code calculates and plots ``g^{(2)}(\tau)`` as a function of ``\tau`` for a coherent, thermal and Fock state:
+
+```@example correlation_and_spectrum
+τlist = LinRange(0, 25, 200)
+
+# Hamiltonian
+N = 25
+a = destroy(N)
+H = 2 * π * a' * a
+
+κ = 0.25
+n_th = 2.0  # bath temperature in terms of excitation number
+c_ops = [
+    sqrt(κ * (1 + n_th)) * a,
+    sqrt(κ *      n_th)  * a'
+]
+
+cases = [
+    Dict("state" => coherent_dm(N, sqrt(2)), "label" => "coherent state", "lstyle" => :solid),
+    Dict("state" => thermal_dm(N, 2), "label" => "thermal state", "lstyle" => :dash),
+    Dict("state" => fock_dm(N, 2), "label" => "Fock state", "lstyle" => :dashdot),
+]
+
+# plot by CairoMakie.jl
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], xlabel = L"Time $\tau$", ylabel = L"g^{(2)}(\tau)")
+
+for case in cases
+    ρ0 = case["state"]
+
+    # calculate the occupation number at τ = 0
+    n0 = expect(a' * a, ρ0)
+
+    # calculate the correlation function g2
+    g2 = correlation_3op_1t(H, ρ0, τlist, c_ops, a', a' * a, a, progress_bar = Val(false))
+    g2 = g2 ./ n0^2
+
+    lines!(ax, τlist, real(g2), label = case["label"], linestyle = case["lstyle"])
+end
+
+axislegend(ax, position = :rt)
+
+fig
+```
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index 676c61e32..6b711223d 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -60,7 +60,7 @@ import DiffEqNoiseProcess: RealWienerProcess!
 # other dependencies (in alphabetical order)
 import ArrayInterface: allowed_getindex, allowed_setindex!
 import Distributed: RemoteChannel
-import FFTW: fft, fftshift
+import FFTW: fft, ifft, fftfreq, fftshift
 import Graphs: connected_components, DiGraph
 import IncompleteLU: ilu
 import Pkg
@@ -107,6 +107,7 @@ include("time_evolution/time_evolution_dynamical.jl")
 # Others
 include("permutation.jl")
 include("correlations.jl")
+include("spectrum.jl")
 include("wigner.jl")
 include("spin_lattice.jl")
 include("arnoldi.jl")
diff --git a/src/correlations.jl b/src/correlations.jl
index ebd3c4666..5b7068d9e 100644
--- a/src/correlations.jl
+++ b/src/correlations.jl
@@ -1,89 +1,95 @@
-export SpectrumSolver, FFTCorrelation, ExponentialSeries
-export correlation_3op_2t, correlation_2op_2t, correlation_2op_1t, spectrum
-
-abstract type SpectrumSolver end
-
-struct FFTCorrelation <: SpectrumSolver end
-
-struct ExponentialSeries{T<:Real,CALC_SS} <: SpectrumSolver
-    tol::T
-    ExponentialSeries(tol::T, calc_steadystate::Bool = false) where {T} = new{T,calc_steadystate}(tol)
+export correlation_3op_2t, correlation_3op_1t, correlation_2op_2t, correlation_2op_1t
+
+function _check_correlation_time_list(tlist::AbstractVector)
+    any(t -> t == 0, tlist) ||
+        throw(ArgumentError("The time list for calculating correlation function must contain the element `0`"))
+    all(>=(0), tlist) ||
+        throw(ArgumentError("All the elements in the time list for calculating correlation function must be positive."))
+    return nothing
 end
 
-ExponentialSeries(; tol = 1e-14, calc_steadystate = false) = ExponentialSeries(tol, calc_steadystate)
-
 @doc raw"""
-    correlation_3op_2t(H::QuantumObject,
-        ψ0::QuantumObject,
-        t_l::AbstractVector,
-        τ_l::AbstractVector,
+    correlation_3op_2t(H::AbstractQuantumObject,
+        ψ0::Union{Nothing,QuantumObject},
+        tlist::AbstractVector,
+        τlist::AbstractVector,
+        c_ops::Union{Nothing,AbstractVector,Tuple},
         A::QuantumObject,
         B::QuantumObject,
-        C::QuantumObject,
-        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
+        C::QuantumObject;
         kwargs...)
 
-Returns the two-times correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle``
+Returns the two-times correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
 
-for a given initial state ``\ket{\psi_0}``.
+If the initial state `ψ0` is given as `nothing`, then the [`steadystate`](@ref) will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
 """
 function correlation_3op_2t(
-    H::QuantumObject{<:AbstractArray{T1},HOpType},
-    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
-    t_l::AbstractVector,
-    τ_l::AbstractVector,
-    A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
-    C::QuantumObject{<:AbstractArray{T5},OperatorQuantumObject},
-    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+    H::AbstractQuantumObject{DataType,HOpType},
+    ψ0::Union{Nothing,QuantumObject{<:AbstractArray{T1},StateOpType}},
+    tlist::AbstractVector,
+    τlist::AbstractVector,
+    c_ops::Union{Nothing,AbstractVector,Tuple},
+    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
+    C::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject};
     kwargs...,
 ) where {
+    DataType,
     T1,
     T2,
     T3,
     T4,
-    T5,
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
 }
-    (H.dims == ψ0.dims && H.dims == A.dims && H.dims == B.dims && H.dims == C.dims) ||
+    # check tlist and τlist
+    _check_correlation_time_list(tlist)
+    _check_correlation_time_list(τlist)
+
+    L = liouvillian(H, c_ops)
+    if ψ0 isa Nothing
+        ψ0 = steadystate(L)
+    end
+
+    allequal((L.dims, ψ0.dims, A.dims, B.dims, C.dims)) ||
         throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
 
-    kwargs2 = merge((saveat = collect(t_l),), (; kwargs...))
-    ρt = mesolve(H, ψ0, t_l, c_ops; kwargs2...).states
+    kwargs2 = merge((saveat = collect(tlist),), (; kwargs...))
+    ρt_list = mesolve(L, ψ0, tlist; kwargs2...).states
 
-    corr = map((t, ρ) -> mesolve(H, C * ρ * A, τ_l .+ t, c_ops, e_ops = [B]; kwargs...).expect[1, :], t_l, ρt)
+    corr = map((t, ρt) -> mesolve(L, C * ρt * A, τlist .+ t, e_ops = [B]; kwargs...).expect[1, :], tlist, ρt_list)
 
-    return corr
+    # make the output correlation Matrix align with QuTiP
+    # 1st dimension corresponds to tlist
+    # 2nd dimension corresponds to τlist
+    return reduce(vcat, transpose.(corr))
 end
 
 @doc raw"""
-    correlation_2op_2t(H::QuantumObject,
-        ψ0::QuantumObject,
-        t_l::AbstractVector,
-        τ_l::AbstractVector,
+    correlation_3op_1t(H::AbstractQuantumObject,
+        ψ0::Union{Nothing,QuantumObject},
+        τlist::AbstractVector,
+        c_ops::Union{Nothing,AbstractVector,Tuple},
         A::QuantumObject,
         B::QuantumObject,
-        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
-        reverse::Bool=false,
+        C::QuantumObject;
         kwargs...)
 
-Returns the two-times correlation function of two operators ``\hat{A}`` and ``\hat{B}`` 
-at different times: ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle``.
+Returns the one-time correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(0) \hat{B}(\tau) \hat{C}(0) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
 
-When `reverse=true`, the correlation function is calculated as ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \right\rangle``.
+If the initial state `ψ0` is given as `nothing`, then the [`steadystate`](@ref) will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
 """
-function correlation_2op_2t(
-    H::QuantumObject{<:AbstractArray{T1},HOpType},
-    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
-    t_l::AbstractVector,
-    τ_l::AbstractVector,
-    A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
-    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    reverse::Bool = false,
+function correlation_3op_1t(
+    H::AbstractQuantumObject{DataType,HOpType},
+    ψ0::Union{Nothing,QuantumObject{<:AbstractArray{T1},StateOpType}},
+    τlist::AbstractVector,
+    c_ops::Union{Nothing,AbstractVector,Tuple},
+    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
+    C::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject};
     kwargs...,
 ) where {
+    DataType,
     T1,
     T2,
     T3,
@@ -91,156 +97,90 @@ function correlation_2op_2t(
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
 }
-    C = eye(prod(H.dims), dims = H.dims)
-    if reverse
-        corr = correlation_3op_2t(H, ψ0, t_l, τ_l, A, B, C, c_ops; kwargs...)
-    else
-        corr = correlation_3op_2t(H, ψ0, t_l, τ_l, C, A, B, c_ops; kwargs...)
-    end
+    corr = correlation_3op_2t(H, ψ0, [0], τlist, c_ops, A, B, C; kwargs...)
 
-    return reduce(hcat, corr)
+    return corr[1, :] # 1 means tlist[1] = 0
 end
 
 @doc raw"""
-    correlation_2op_1t(H::QuantumObject,
-        ψ0::QuantumObject,
-        τ_l::AbstractVector,
+    correlation_2op_2t(H::AbstractQuantumObject,
+        ψ0::Union{Nothing,QuantumObject},
+        tlist::AbstractVector,
+        τlist::AbstractVector,
+        c_ops::Union{Nothing,AbstractVector,Tuple},
         A::QuantumObject,
-        B::QuantumObject,
-        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
+        B::QuantumObject;
         reverse::Bool=false,
         kwargs...)
 
-Returns the one-time correlation function of two operators ``\hat{A}`` and ``\hat{B}`` at different times ``\left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle``.
+Returns the two-times correlation function of two operators ``\hat{A}`` and ``\hat{B}`` : ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
 
-When `reverse=true`, the correlation function is calculated as ``\left\langle \hat{A}(0) \hat{B}(\tau) \right\rangle``.
+If the initial state `ψ0` is given as `nothing`, then the [`steadystate`](@ref) will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
+
+When `reverse=true`, the correlation function is calculated as ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \right\rangle``.
 """
-function correlation_2op_1t(
-    H::QuantumObject{<:AbstractArray{T1},HOpType},
-    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
-    τ_l::AbstractVector,
-    A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
-    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+function correlation_2op_2t(
+    H::AbstractQuantumObject{DataType,HOpType},
+    ψ0::Union{Nothing,QuantumObject{<:AbstractArray{T1},StateOpType}},
+    tlist::AbstractVector,
+    τlist::AbstractVector,
+    c_ops::Union{Nothing,AbstractVector,Tuple},
+    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject};
     reverse::Bool = false,
     kwargs...,
 ) where {
+    DataType,
     T1,
     T2,
     T3,
-    T4,
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
 }
-    corr = correlation_2op_2t(H, ψ0, [0], τ_l, A, B, c_ops; reverse = reverse, kwargs...)
+    C = eye(prod(H.dims), dims = H.dims)
+    if reverse
+        corr = correlation_3op_2t(H, ψ0, tlist, τlist, c_ops, A, B, C; kwargs...)
+    else
+        corr = correlation_3op_2t(H, ψ0, tlist, τlist, c_ops, C, A, B; kwargs...)
+    end
 
-    return corr[:, 1]
+    return corr
 end
 
 @doc raw"""
-    spectrum(H::QuantumObject,
-        ω_list::AbstractVector,
-        A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-        B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
-        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
-        solver::MySolver=ExponentialSeries(),
+    correlation_2op_1t(H::AbstractQuantumObject,
+        ψ0::Union{Nothing,QuantumObject},
+        τlist::AbstractVector,
+        c_ops::Union{Nothing,AbstractVector,Tuple},
+        A::QuantumObject,
+        B::QuantumObject;
+        reverse::Bool=false,
         kwargs...)
 
-Returns the emission spectrum 
+Returns the one-time correlation function of two operators ``\hat{A}`` and ``\hat{B}`` : ``\left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
 
-```math
-S(\omega) = \int_{-\infty}^\infty \left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle e^{-i \omega \tau} d \tau
-```
+If the initial state `ψ0` is given as `nothing`, then the [`steadystate`](@ref) will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
+
+When `reverse=true`, the correlation function is calculated as ``\left\langle \hat{A}(0) \hat{B}(\tau) \right\rangle``.
 """
-function spectrum(
-    H::QuantumObject{MT1,HOpType},
-    ω_list::AbstractVector,
+function correlation_2op_1t(
+    H::AbstractQuantumObject{DataType,HOpType},
+    ψ0::Union{Nothing,QuantumObject{<:AbstractArray{T1},StateOpType}},
+    τlist::AbstractVector,
+    c_ops::Union{Nothing,AbstractVector,Tuple},
     A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
-    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    solver::MySolver = ExponentialSeries(),
+    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject};
+    reverse::Bool = false,
     kwargs...,
 ) where {
-    MT1<:AbstractMatrix,
+    DataType,
+    T1,
     T2,
     T3,
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
-    MySolver<:SpectrumSolver,
+    StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
 }
-    return _spectrum(H, ω_list, A, B, c_ops, solver; kwargs...)
-end
-
-function _spectrum(
-    H::QuantumObject{<:AbstractArray{T1},HOpType},
-    ω_list::AbstractVector,
-    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
-    c_ops,
-    solver::FFTCorrelation;
-    kwargs...,
-) where {T1,T2,T3,HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
-    Nsamples = length(ω_list)
-    ω_max = abs(maximum(ω_list))
-    dω = 2 * ω_max / (Nsamples - 1)
-    ω_l = -ω_max:dω:ω_max
-
-    T = 2π / (ω_l[2] - ω_l[1])
-    τ_l = range(0, T, length = length(ω_l))
-
-    ρss = steadystate(H, c_ops)
-    corr = correlation_2op_1t(H, ρss, τ_l, A, B, c_ops; kwargs...)
-
-    S = fftshift(fft(corr)) / length(τ_l)
-
-    return ω_l, 2 .* real.(S)
-end
-
-function _spectrum_get_rates_vecs_ss(L, solver::ExponentialSeries{T,true}) where {T}
-    result = eigen(L)
-    rates, vecs = result.values, result.vectors
-
-    return rates, vecs, steadystate(L).data
-end
-
-function _spectrum_get_rates_vecs_ss(L, solver::ExponentialSeries{T,false}) where {T}
-    result = eigen(L)
-    rates, vecs = result.values, result.vectors
-
-    ss_idx = findmin(abs2, rates)[2]
-    ρss = vec2mat(@view(vecs[:, ss_idx]))
-    ρss = (ρss + ρss') / 2
-    ρss ./= tr(ρss)
-
-    return rates, vecs, ρss
-end
-
-function _spectrum(
-    H::QuantumObject{<:AbstractArray{T1},HOpType},
-    ω_list::AbstractVector,
-    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
-    c_ops,
-    solver::ExponentialSeries;
-    kwargs...,
-) where {T1,T2,T3,HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
-    (H.dims == A.dims == B.dims) || throw(DimensionMismatch("The dimensions of H, A and B must be the same"))
-
-    L = liouvillian(H, c_ops)
-
-    ω_l = ω_list
-
-    rates, vecs, ρss = _spectrum_get_rates_vecs_ss(L, solver)
-
-    ρ0 = B.data * ρss
-    v = vecs \ mat2vec(ρ0)
-
-    amps = map(i -> v[i] * tr(A.data * vec2mat(@view(vecs[:, i]))), eachindex(rates))
-    idxs = findall(x -> abs(x) > solver.tol, amps)
-    amps, rates = amps[idxs], rates[idxs]
-
-    # spec = map(ω -> 2 * real(sum(@. amps * (1 / (1im * ω - rates)))), ω_l)
-    amps_rates = zip(amps, rates)
-    spec = map(ω -> 2 * real(sum(x -> x[1] / (1im * ω - x[2]), amps_rates)), ω_l)
+    corr = correlation_2op_2t(H, ψ0, [0], τlist, c_ops, A, B; reverse = reverse, kwargs...)
 
-    return ω_l, spec
+    return corr[1, :] # 1 means tlist[1] = 0
 end
diff --git a/src/deprecated.jl b/src/deprecated.jl
index 8cc4db3e5..98d92193d 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -13,3 +13,93 @@ function deprecated_foo(args...; kwargs...)
     error("`deprecated_foo` has been deprecated and will be removed in next major release, please use `new_foo` instead.")
 end
 =#
+
+export FFTCorrelation
+
+FFTCorrelation() = error(
+    "`FFTCorrelation` has been deprecated and will be removed in next major release, please use `spectrum_correlation_fft` to calculate the spectrum with FFT method instead.",
+)
+
+correlation_3op_2t(
+    H::QuantumObject{<:AbstractArray{T1},HOpType},
+    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
+    t_l::AbstractVector,
+    τ_l::AbstractVector,
+    A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
+    C::QuantumObject{<:AbstractArray{T5},OperatorQuantumObject},
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+    kwargs...,
+) where {
+    T1,
+    T2,
+    T3,
+    T4,
+    T5,
+    HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
+    StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
+} = error(
+    "The parameter order of `correlation_3op_2t` has been changed, please use `?correlation_3op_2t` to check the updated docstring.",
+)
+
+correlation_3op_1t(
+    H::QuantumObject{<:AbstractArray{T1},HOpType},
+    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
+    τ_l::AbstractVector,
+    A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
+    C::QuantumObject{<:AbstractArray{T5},OperatorQuantumObject},
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+    kwargs...,
+) where {
+    T1,
+    T2,
+    T3,
+    T4,
+    T5,
+    HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
+    StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
+} = error(
+    "The parameter order of `correlation_3op_1t` has been changed, please use `?correlation_3op_1t` to check the updated docstring.",
+)
+
+correlation_2op_2t(
+    H::QuantumObject{<:AbstractArray{T1},HOpType},
+    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
+    t_l::AbstractVector,
+    τ_l::AbstractVector,
+    A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+    reverse::Bool = false,
+    kwargs...,
+) where {
+    T1,
+    T2,
+    T3,
+    T4,
+    HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
+    StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
+} = error(
+    "The parameter order of `correlation_2op_2t` has been changed, please use `?correlation_2op_2t` to check the updated docstring.",
+)
+
+correlation_2op_1t(
+    H::QuantumObject{<:AbstractArray{T1},HOpType},
+    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
+    τ_l::AbstractVector,
+    A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+    reverse::Bool = false,
+    kwargs...,
+) where {
+    T1,
+    T2,
+    T3,
+    T4,
+    HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
+    StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
+} = error(
+    "The parameter order of `correlation_2op_1t` has been changed, please use `?correlation_2op_1t` to check the updated docstring.",
+)
diff --git a/src/spectrum.jl b/src/spectrum.jl
new file mode 100644
index 000000000..a62a69ca6
--- /dev/null
+++ b/src/spectrum.jl
@@ -0,0 +1,177 @@
+export spectrum, spectrum_correlation_fft
+export SpectrumSolver, ExponentialSeries, PseudoInverse
+
+abstract type SpectrumSolver end
+
+@doc raw"""
+    ExponentialSeries(; tol = 1e-14, calc_steadystate = false)
+
+A solver which solves [`spectrum`](@ref) by finding the eigen decomposition of the Liouvillian [`SuperOperator`](@ref) and calculate the exponential series.
+"""
+struct ExponentialSeries{T<:Real,CALC_SS} <: SpectrumSolver
+    tol::T
+    ExponentialSeries(tol::T, calc_steadystate::Bool = false) where {T} = new{T,calc_steadystate}(tol)
+end
+
+ExponentialSeries(; tol = 1e-14, calc_steadystate = false) = ExponentialSeries(tol, calc_steadystate)
+
+@doc raw"""
+    PseudoInverse(; alg::SciMLLinearSolveAlgorithm = KrylovJL_GMRES())
+
+A solver which solves [`spectrum`](@ref) by finding the inverse of Liouvillian [`SuperOperator`](@ref) using the `alg`orithms given in [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/).
+"""
+struct PseudoInverse{MT<:SciMLLinearSolveAlgorithm} <: SpectrumSolver
+    alg::MT
+end
+
+PseudoInverse(; alg::SciMLLinearSolveAlgorithm = KrylovJL_GMRES()) = PseudoInverse(alg)
+
+@doc raw"""
+    spectrum(H::QuantumObject,
+        ωlist::AbstractVector,
+        c_ops::Union{Nothing,AbstractVector,Tuple},
+        A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
+        B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject};
+        solver::SpectrumSolver=ExponentialSeries(),
+        kwargs...)
+
+Calculate the spectrum of the correlation function
+
+```math
+S(\omega) = \int_{-\infty}^\infty \lim_{t \rightarrow \infty} \left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle e^{-i \omega \tau} d \tau
+```
+
+See also the following list for `SpectrumSolver` docstrings:
+- [`ExponentialSeries`](@ref)
+- [`PseudoInverse`](@ref)
+"""
+function spectrum(
+    H::QuantumObject{MT1,HOpType},
+    ωlist::AbstractVector,
+    c_ops::Union{Nothing,AbstractVector,Tuple},
+    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject};
+    solver::SpectrumSolver = ExponentialSeries(),
+    kwargs...,
+) where {MT1<:AbstractMatrix,T2,T3,HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
+    return _spectrum(liouvillian(H, c_ops), ωlist, A, B, solver; kwargs...)
+end
+
+function _spectrum_get_rates_vecs_ss(L, solver::ExponentialSeries{T,true}) where {T}
+    result = eigen(L)
+    rates, vecs = result.values, result.vectors
+
+    return rates, vecs, steadystate(L).data
+end
+
+function _spectrum_get_rates_vecs_ss(L, solver::ExponentialSeries{T,false}) where {T}
+    result = eigen(L)
+    rates, vecs = result.values, result.vectors
+
+    ss_idx = findmin(abs2, rates)[2]
+    ρss = vec2mat(@view(vecs[:, ss_idx]))
+    ρss = (ρss + ρss') / 2
+    ρss ./= tr(ρss)
+
+    return rates, vecs, ρss
+end
+
+function _spectrum(
+    L::QuantumObject{<:AbstractArray{T1},SuperOperatorQuantumObject},
+    ωlist::AbstractVector,
+    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
+    solver::ExponentialSeries;
+    kwargs...,
+) where {T1,T2,T3}
+    allequal((L.dims, A.dims, B.dims)) ||
+        throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+
+    rates, vecs, ρss = _spectrum_get_rates_vecs_ss(L, solver)
+
+    ρ0 = B.data * ρss
+    v = vecs \ mat2vec(ρ0)
+
+    amps = map(i -> v[i] * tr(A.data * vec2mat(@view(vecs[:, i]))), eachindex(rates))
+    idxs = findall(x -> abs(x) > solver.tol, amps)
+    amps, rates = amps[idxs], rates[idxs]
+
+    # spec = map(ω -> 2 * real(sum(@. amps * (1 / (1im * ω - rates)))), ωlist)
+    amps_rates = zip(amps, rates)
+    spec = map(ω -> 2 * real(sum(x -> x[1] / (1im * ω - x[2]), amps_rates)), ωlist)
+
+    return spec
+end
+
+function _spectrum(
+    L::QuantumObject{<:AbstractArray{T1},SuperOperatorQuantumObject},
+    ωlist::AbstractVector,
+    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
+    solver::PseudoInverse;
+    kwargs...,
+) where {T1,T2,T3}
+    allequal((L.dims, A.dims, B.dims)) ||
+        throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+
+    ωList = convert(Vector{_FType(L)}, ωlist) # Convert it to support GPUs and avoid type instabilities
+    Length = length(ωList)
+    spec = Vector{_FType(L)}(undef, Length)
+
+    # calculate vectorized steadystate, multiply by operator B on the left (spre)
+    ρss = mat2vec(steadystate(L))
+    b = (spre(B) * ρss).data
+
+    # multiply by operator A on the left (spre) and then perform trace operation
+    D = prod(L.dims)
+    _tr = SparseVector(D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(_CType(L), D)) # same as vec(system_identity_matrix)
+    _tr_A = transpose(_tr) * spre(A).data
+
+    cache = nothing
+    I_cache = I(D^2)
+    for (idx, ω) in enumerate(ωList)
+        if idx == 1
+            cache = init(LinearProblem(L.data - 1im * ω * I_cache, b), solver.alg, kwargs...)
+            sol = solve!(cache)
+        else
+            cache.A = L.data - 1im * ω * I_cache
+            sol = solve!(cache)
+        end
+
+        # trace over the Hilbert space of system (expectation value)
+        spec[idx] = -2 * real(dot(_tr_A, sol.u))
+    end
+
+    return spec
+end
+
+@doc raw"""
+    spectrum_correlation_fft(tlist, corr; inverse=false)
+
+Calculate the power spectrum corresponding to a two-time correlation function using fast Fourier transform (FFT).
+
+# Parameters
+- `tlist::AbstractVector`: List of times at which the two-time correlation function is given.
+- `corr::AbstractVector`: List of two-time correlations corresponding to the given time point in `tlist`.
+- `inverse::Bool`: Whether to use the inverse Fourier transform or not. Default to `false`.
+
+# Returns
+- `ωlist`: the list of angular frequencies ``\omega``.
+- `Slist`: the list of the power spectrum corresponding to the angular frequencies in `ωlist`.
+"""
+function spectrum_correlation_fft(tlist::AbstractVector, corr::AbstractVector; inverse::Bool = false)
+    N = length(tlist)
+    dt_list = diff(tlist)
+    dt = dt_list[1]
+
+    all(≈(dt), dt_list) || throw(ArgumentError("tlist must be equally spaced for FFT."))
+
+    # power spectrum list
+    F = inverse ? N * ifft(corr) : fft(corr)
+    Slist = 2 * dt * real(fftshift(F))
+
+    # angular frequency list
+    ωlist = 2 * π * fftshift(fftfreq(N, 1 / dt))
+
+    return ωlist, Slist
+end
diff --git a/test/core-test/correlations_and_spectrum.jl b/test/core-test/correlations_and_spectrum.jl
index 5c5079096..381cd9e08 100644
--- a/test/core-test/correlations_and_spectrum.jl
+++ b/test/core-test/correlations_and_spectrum.jl
@@ -1,13 +1,22 @@
 @testset "Correlations and Spectrum" begin
-    a = destroy(10)
+    N = 10
+    Id = qeye(N)
+    a = destroy(N)
     H = a' * a
     c_ops = [sqrt(0.1 * (0.01 + 1)) * a, sqrt(0.1 * (0.01)) * a']
 
-    ω_l = range(0, 3, length = 1000)
-    ω_l1, spec1 = spectrum(H, ω_l, a', a, c_ops, solver = FFTCorrelation(), progress_bar = Val(false))
-    ω_l2, spec2 = spectrum(H, ω_l, a', a, c_ops)
+    t_l = range(0, 333 * π, length = 1000)
+    corr1 = correlation_2op_1t(H, nothing, t_l, c_ops, a', a; progress_bar = Val(false))
+    corr2 = correlation_3op_1t(H, nothing, t_l, c_ops, Id, a', a; progress_bar = Val(false))
+    ω_l1, spec1 = spectrum_correlation_fft(t_l, corr1)
+
+    ω_l2 = range(0, 3, length = 1000)
+    spec2 = spectrum(H, ω_l2, c_ops, a', a)
+    spec3 = spectrum(H, ω_l2, c_ops, a', a; solver = PseudoInverse())
+
     spec1 = spec1 ./ maximum(spec1)
     spec2 = spec2 ./ maximum(spec2)
+    spec3 = spec3 ./ maximum(spec3)
 
     test_func1 = maximum(real.(spec1)) * (0.1 / 2)^2 ./ ((ω_l1 .- 1) .^ 2 .+ (0.1 / 2)^2)
     test_func2 = maximum(real.(spec2)) * (0.1 / 2)^2 ./ ((ω_l2 .- 1) .^ 2 .+ (0.1 / 2)^2)
@@ -15,9 +24,36 @@
     idxs2 = test_func2 .> 0.05
     @test sum(abs2.(spec1[idxs1] .- test_func1[idxs1])) / sum(abs2.(test_func1[idxs1])) < 0.01
     @test sum(abs2.(spec2[idxs2] .- test_func2[idxs2])) / sum(abs2.(test_func2[idxs2])) < 0.01
+    @test all(corr1 .≈ corr2)
+    @test all(spec2 .≈ spec3)
 
     @testset "Type Inference spectrum" begin
-        @inferred spectrum(H, ω_l, a', a, c_ops, solver = FFTCorrelation(), progress_bar = Val(false))
-        @inferred spectrum(H, ω_l, a', a, c_ops)
+        @inferred correlation_2op_1t(H, nothing, t_l, c_ops, a', a; progress_bar = Val(false))
+        @inferred spectrum_correlation_fft(t_l, corr1)
+        @inferred spectrum(H, ω_l2, c_ops, a', a)
+        @inferred spectrum(H, ω_l2, c_ops, a', a; solver = PseudoInverse())
+    end
+
+    # tlist and τlist checks
+    t_fft_wrong = [0, 1, 10]
+    t_wrong1 = [1, 2, 3]
+    t_wrong2 = [-1, 0, 1]
+    @test_throws ArgumentError spectrum_correlation_fft(t_fft_wrong, corr1)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_l, t_wrong1, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_l, t_wrong2, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong1, t_l, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong2, t_l, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong1, t_wrong1, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong1, t_wrong2, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong2, t_wrong1, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong2, t_wrong2, c_ops, Id, a', a)
+
+    @testset "Deprecated Errors" begin
+        ρ0 = rand_dm(N)
+        @test_throws ErrorException FFTCorrelation()
+        @test_throws ErrorException correlation_3op_2t(H, ρ0, t_l, t_l, a, a', a, c_ops)
+        @test_throws ErrorException correlation_3op_1t(H, ρ0, t_l, a, a', a, c_ops)
+        @test_throws ErrorException correlation_2op_2t(H, ρ0, t_l, t_l, a', a, c_ops)
+        @test_throws ErrorException correlation_2op_1t(H, ρ0, t_l, a', a, c_ops)
     end
 end