Oscillations of magma in volcanic conduits are thought to be the source of certain seismic and infrasonic signals observed near active volcanoes. However, the multiphase and stratified nature of magma within the conduit complicates the calculation of resonant modes that is required to interpret observations. Here we present a linearized mathematical framework to describe small-amplitude oscillations and waves in a stably stratified column of two-phase magma (liquid melt and gas bubbles) with a traction-free upper surface (a lava lake). We explore the role of time-dependent mass exchange between the phases, depth-varying fluid properties and gravity on the modes of oscillation of inviscid magma within an axisymmetric, vertical conduit. Non-equilibrium phase exchange, which we refer to as bubble growth and resorption (BGR), is parameterized by introduction of a kinetic time scale quantifying mass exchange between the liquid and gas phases that evolves the mixture towards a state of thermodynamic equilibrium. Using a provably stable finite difference method, we solve the eigenvalue problem for the resonance frequencies, decay rates, and spatial structure of the conduit eigenmodes. The numerical method is then extended to time-domain simulations of waves excited by internal volumetric sources in the conduit or forces applied to the surface of the lava lake. We connect time-dependent wave propagation simulations to the modal analysis by identifying the primary modes that are excited by representative excitation processes. Waves propagating through bubbly magma are dispersive, and their behaviour is determined by three dimensionless parameters. One quantifies the importance of buoyancy and gravitational restoring forces relative to compressibility, the second quantifies differences between fluid properties (e.g. mixture compressibility) under equilibrium and non-equilibrium conditions, and the third compares the wave period to the BGR time scale. Pronounced depth variations in background fluid properties, such as the transition from liquid melt with dissolved volatiles at the high pressures at depth to bubbly magma above the gas exsolution depth, segment the conduit into distinct regions. The longest-period modes, which are expressed with the largest amplitudes for typical excitation processes, are most sensitive to the length of the bubbly region and properties of the bubbly magma within it. While the boundary condition at the bottom of the conduit determines whether the fundamental mode is affected by the total conduit length, modes localized above the exsolution depth are remarkably insensitive to the overall conduit length. Our analysis suggests that parameters affecting eruption style, such as total volatile content and kinetic time scales of BGR, along with excitation source characteristics, are imprinted on long-period seismic and infrasonic signals at active volcanoes.