We analyzed over 3000 Fourier spectra from 370 earthquakes of energy magnitude ( M E) 1.1–6.0 recorded by the Southern Ontario Seismic Network (SOSN)/POLARIS networks during the period 1991–2010 in the area of southern Ontario and western Quebec. We employed a range of velocity stacking methods to significantly reduce the problem of variability due to wave scattering. This enabled us to determine underlying nonrandom spectral features, including source effects, site effects, and anelastic attenuation effects on spectral shape. The analysis technique is that we stack the velocity spectra of the whole observed data set into one or two bins and then compare that sum (the observed stack) with the theoretical expectation for corresponding stacks of simulated signals (the theoretical stack) for a given set of input parameters. A grid‐search technique is used to find the input‐parameter combination that optimizes the agreement between the observed and theoretical stacks. By stacking the spectra in different ways, different underlying spectral features are explored. We find the method works surprisingly well, allowing us to determine the apparent anelastic attenuation effects on the spectral shape, the average effect of site response, and some basic features of the source spectra. The key results of our paper: (1) there is no unique pair of values of the coefficients Q and n of the frequency‐dependent Quality factor relationship Q = Q f n , but there exist pairs of Q and n along a curve in Q – n space that are equivalent in terms of their effect on spectral shape; (2) the relationship between log corner frequency and energy magnitude ( M E) is linear, with a slope close to (−0.22) that is consistent with constant‐width faulting for the studied small‐to‐moderate events; (3) the relation between moment magnitude M and energy magnitude M E was found to be M = 8/9 M E.