Numerical methods for time-harmonic wave propagation phenomena are often computationally intensive, especially in mid- and high-frequency regimes, making a direct frequency-response analysis prohibitively expensive. Starting from few expensive solves of the problem, model order reduction (MOR) methods can provide a reliable and cheap-to-evaluate approximation of the frequency response of the system. In this talk, we describe an MOR approach for parametric frequency-response problems, where the high-fidelity problem models not only the impact of the frequency z on the system response, but also that of additional (e.g., design or uncertain) parameters p. Our proposed method relies on minimal rational interpolation for the surrogate modeling of the frequency dependence, for few fixed values of the parameters p. In this step, only the dependence on z is taken into account, and the resulting surrogates are a collection of rational functions of frequency only. Then, the different z-surrogates are combined to obtain a global approximation with respect to both frequency and parameters. Our approach is non-intrusive, i.e., we do not require access to the matrices/operators defining the underlying problem. Moreover, we do not assume the high-fidelity problem to have a specific structure. This is of practical interest, for instance, when the (finite element) solver that computes the frequency response is a black box, which happens when closed-source proprietary code is used. Similarly, non-intrusive methods may be attractive when the data comes from lab experiments rather than from numerical simulations. Moreover, we note that, by its non-intrusive nature, our method is easily applied even when nonlinearities are present. We also describe how, in the interest of making our approach more accurate and efficient, one can select the sample points adaptively, still within a non-intrusive framework. This is achieved by using locally adaptive sparse grids over p-space, which, among other benefits, weaken the curse of dimension that is incurred if many parameters are present. The resulting method is referred to as "double-greedy parametric minimal rational interpolation", since adaptivity acts both on frequency and on the parameters. We present numerical examples in the modeling of electrical circuits and of elasto-dynamic systems, providing evidence of the approximation quality and computational efficiency of the surrogate models obtained with the proposed technique.