The horizontal-to-vertical spectral ratio (HVSR) of ambient noise is commonly used to infer a site’s resonance frequency (f0,site). HVSR calculations are performed most commonly using the Fourier amplitude spectrum obtained from a single merged horizontal component (e.g. the geometric mean component) from a three-component sensor. However, the use of a single merged horizontal component implicitly relies on the assumptions of azimuthally isotropic seismic noise and 1D surface and subsurface conditions. These assumptions may not be justified at many sites, leading to azimuthal variability in HVSR measurements that cannot be accounted for using a single merged component. This paper proposes a new statistical method to account for azimuthal variability in the peak frequency of HVSR curves (f0,HVSR). The method uses rotated horizontal components at evenly distributed azimuthal intervals to investigate and quantify azimuthal variability. To ensure unbiased statistics for f0,HVSR are obtained, a frequency-domain window-rejection algorithm is applied at each azimuth to automatically remove contaminated time windows in which the f0,HVSR values are statistical outliers relative to those obtained from the majority of windows at that azimuth. Then, a weighting scheme is used to account for different numbers of accepted time windows at each azimuth. The new method is applied to a dataset of 114 HVSR measurements with significant azimuthal variability in f0,HVSR, and is shown to reliably account for this variability. The methodology is also extended to the estimation of a complete lognormal-median HVSR curve that accounts for azimuthal variability. To encourage the adoption of this statistical approach to accounting for azimuthal variability in single-station HVSR measurements, the methods presented in this paper have been incorporated into hvsrpy, an open-source Python package for HVSR processing.