Selection of the amplitude of magnetic field modulation for continuous wave electron paramagnetic resonance (EPR) often is a trade-off between sensitivity and resolution. Increasing the modulation amplitude improves the signal-to-noise ratio, S/N, at the expense of broadening the signal. Combining information from multiple harmonics of the field-modulated signal is proposed as a method to obtain the first derivative spectrum with minimal broadening and improved signal-to-noise. The harmonics are obtained by digital phase-sensitive detection of the signal at the modulation frequency and its integer multiples. Reconstruction of the first derivative EPR line is done in the Fourier conjugate domain where each harmonic can be represented as the product of the Fourier transform of the 1st derivative signal with an analytical function. The analytical function for each harmonic can be viewed as a filter. The Fourier transform of the 1st derivative spectrum can be calculated from all available harmonics by solving an optimization problem with the goal of maximizing the S/N. Inverse Fourier transformation of the result produces the 1st derivative EPR line in the magnetic field domain. The use of modulation amplitude greater than linewidth improves the S/N, but does not broaden the reconstructed spectrum. The method works for an arbitrary EPR line shape, but is limited to the case when magnetization instantaneously follows the modulation field, which is known as the adiabatic approximation.