We construct a new analytic model for the description of the gravitational wave (GW) emission in the post-merger phase of binary neutron star (BNS) mergers. The model consists of exponentially decaying sinusoidal functions, attributed to the various oscillation modes, quasi-linear combination tones, or non-linear features that appear in the post-merger phase. We consider a time dependence of the main post-merger frequency peak which is described by a two-segment linear expression. The effectiveness of the model is assessed, in terms of the fitting factor, along a sequence of equal-mass simulations of varying total binary mass. We identify new spectral features, appearing in high-mass configurations, originating from a non-linear coupling between the quasi-radial oscillation and antipodal tidal deformations. The model achieves high fitting factors, and so can be used for the parameter estimation of detections in upcoming searches with aLIGO+ and aVirgo+, or with future detectors such as Einstein Telescope, Cosmic Explorer, or with high-frequency detectors.