The topic of functional connectivity in neuroimaging is expanding rapidly and many studies now focus on coupling between spatially separate brain regions. These studies show that a relatively small number of large scale networks exist within the brain, and that healthy function of these networks is disrupted in many clinical populations. To date, the vast majority of studies probing connectivity employ techniques that compute time averaged correlation over several minutes, and between specific pre-defined brain locations. However, increasing evidence suggests that functional connectivity is non-stationary in time. Further, electrophysiological measurements show that connectivity is dependent on the frequency band of neural oscillations. It is also conceivable that networks exhibit a degree of spatial inhomogeneity, i.e. the large scale networks that we observe may result from the time average of multiple transiently synchronised sub-networks, each with their own spatial signature. This means that the next generation of neuroimaging tools to compute functional connectivity must account for spatial inhomogeneity, spectral non-uniformity and temporal non-stationarity. Here, we present a means to achieve this via application of windowed canonical correlation analysis (CCA) to source space projected MEG data. We describe the generation of time–frequency connectivity plots, showing the temporal and spectral distribution of coupling between brain regions. Moreover, CCA over voxels provides a means to assess spatial non-uniformity within short time–frequency windows. The feasibility of this technique is demonstrated in simulation and in a resting state MEG experiment where we elucidate multiple distinct spatio-temporal-spectral modes of covariation between the left and right sensorimotor areas.