Resolvent analysis identifies the most responsive forcings and most receptive states of a dynamical system, in an input--output sense, based on its governing equations. Interest in the method has continued to grow during the past decade due to its potential to reveal structures in turbulent flows, to guide sensor/actuator placement, and for flow control applications. However, resolvent analysis requires access to high-fidelity numerical solvers to produce the linearized dynamics operator. In this work, we develop a purely data-driven algorithm to perform resolvent analysis to obtain the leading forcing and response modes, without recourse to the governing equations, but instead based on snapshots of the transient evolution of linearly stable flows. The formulation of our method follows from two established facts: $1)$ dynamic mode decomposition can approximate eigenvalues and eigenvectors of the underlying operator governing the evolution of a system from measurement data, and $2)$ a projection of the resolvent operator onto an invariant subspace can be built from this learned eigendecomposition. We demonstrate the method on numerical data of the linearized complex Ginzburg--Landau equation and of three-dimensional transitional channel flow, and discuss data requirements. The ability to perform resolvent analysis in a completely equation-free and adjoint-free manner will play a significant role in lowering the barrier of entry to resolvent research and applications.