Harmonic Amplitude Summation for Frequency-tagging Analysis

It has long been known that a stimulus presented at a periodic rate elicits a response from an observer's brain at exactly that rate. For example, a light flickering on and off at a periodic rate, 14 times a second, elicits a measurable response in the EEG of a human observer 14 times a second (Adrian & Mathews, 1934). In the time domain, a response is evident as periodic changes in the brain's response amplitude across time. Following Fourier transformation into a frequency domain representation (Danielson & Lanczos, 1942; Fourier, 1822), the response is evident as a high amplitude “peak” at exactly the fundamental stimulus presentation rate (frequency = F) and/or its higher harmonics, that is, at frequencies that are integer multiples of F (2F, 3F, etc.; Regan, 1966, 1989).

The approach of presenting stimuli and analyzing neural responses at the frequency of stimulation is referred to by many names: “Frequency tagging” (Srinivasan, Russell, Edelman, & Tononi, 1999; Tononi, Srinivasan, Russell, & Edelman, 1998) is the one that will be used here. Other names for this approach differ mainly on their point of reference: to the responses that appear consistently periodic to stimuli presented at high rates, that is, “steady-state” responses, for example, “steady-state visual-evoked potentials” (Norcia, Appelbaum, Ales, Cottereau, & Rossion, 2015; Heinrich, 2010; Di Russo, Teder-Salajarvi, & Hillyard, 2002; Regan, 1966, 1989) and “auditory steady state potentials/responses” (Watkin, 2008; Geisler, 1960) or “traveling wave” responses (Engel, Glover, & Wandell, 1997); to the “stimulation” mode itself (“fast periodic visual stimulation”; Rossion, 2014; Rossion, Retter, & Liu-Shuang, 2020); or the “analysis” occurring in the frequency domain (“Fourier analysis/synthesis”; Zhou, Melloni, Poeppel, & Ding, 2016; Bach & Meigen, 1999; Movshon, Thompson, & Tolhurst, 1978; or simply “frequency (domain) analysis,” for example, as in McKeefry, Russell, Murray, & Kulikowski, 1996). Despite the varying terminologies, the principles of the approach are the same. In a similar vein, various types of stimulation modalities (visual, auditory, somatosensory, cross-modal) and recording methods (EEG, electroretinogram, fMRI, single-cell recordings, etc.) may be applied with various participant groups (human adults, children, infants, nonhuman primates, cats, rodents, frogs, insects, etc.), resulting in some practical differences, but the same fundamentals, of the approach.

In (cognitive) neuroscience research, the frequency-tagging approach is associated with undeniable advantages. As noted early on, this approach is well suited for specifically relating brain processes to external events: “This gives a method of tracing the visual messages in the brain, for by means of the flicker rhythm they can be made easy to recognize” (Adrian, 1944, p. 361). More recently, its objectivity and sensitivity (i.e., high signal-to-noise ratio) have been highlighted, and the use of the paradigm is undoubtedly on the rise, having been extended from the study of basic sensory processes and their modulation by spatial/selective attention to the direct measurement of higher levels of cognition in recent years (see Norcia et al., 2015, for a review). However, frequency tagging is still fundamentally limited by outstanding conceptual and methodological ambiguities in dealing with responses occurring across harmonics.

Another way of describing frequency tagging is the following: Given a periodic stimulus, responses of the brain periodic to that stimulus are investigated. In this formulation, it is evident that the brain's responses may occur at the rate of stimulation, F, but also at the other rates periodic to the stimulation: the higher harmonics (2F, 3F, etc.). For example, a stimulus modulated 8 times a second, at 8 Hz, may generate responses that are evident as amplitude peaks in the frequency domain representation of the brain recording at 8 Hz (F, the first harmonic, corresponding to the fundamental frequency¹), but also at 16 Hz (2F, the second harmonic) and 24 Hz (3F, the third harmonic). Because only responses at higher harmonics are periodic to the fundamental frequency, it is uniquely at the higher harmonics, rather than at a diffuse band, that higher frequency constituents of frequency-tagged brain responses are present.

Although responses of the brain are not always generated at the higher harmonics, they often do occur (Rossion et al., 2020; Zhou et al., 2016; Norcia et al., 2015; Rossion, 2014; Heinrich, 2010; Vialatte, Maurice, Dauwels, & Cichocki, 2010; Bach & Meigen, 1999; Regan, 1966). Note that responses are not always generated at F either; for a classic example, in the case of alternating symmetrical stimulus inputs (e.g., pattern-reversing checkerboards), the brain responds only at 2F and higher even harmonics (Cobb, Morton, & Ettlinger, 1967; reviewed in Norcia et al., 2015; for different examples, Zhou et al., 2016; Heinrich, 2010; Movshon et al., 1978). Furthermore, note that, throughout this article, only harmonics that are specific to their fundamental frequency are addressed, which is always the case when a single stimulus presentation frequency is used (but for an extension to other cases, with more complex stimulation paradigms, please see the Which Harmonics to Consider? section).

At present, although higher harmonic responses are an integral part of brain responses, they are not systematically addressed in frequency-tagging research. In many studies, higher harmonic responses are not even reported (e.g., Bekhtereva, Pritschmann, Keil, & Müller, 2018; Min, Dähne, Ahn, Noh, & Müller, 2016; Paulk, Kirszenblat, Zhou, & van Swinderen, 2015; Coia, Jones, Duncan, & Crognale, 2014; Kuś et al., 2013; Cottereau et al., 2011; Wattam-Bell et al., 2010; Di Russo et al., 2007; Müller et al., 2006; Braddick, Birtles, Wattam-Bell, & Atkinson, 2005; Chen, Seth, Gally, & Edelman, 2003; Heinrich & Bach, 2001; Regan & Regan, 1988a; Tononi et al., 1998; Peterzell & Norcia, 1997; Müller, Teder, & Hillyard, 1997; Morgan, Hansen, & Hillyard, 1996, or are extirpated by narrow [band-pass, Gabor, etc.] filtering, e.g., Davidson, Mithen, Hogendoorn, van Boxtel, & Tsuchiya, 2020; Miskociv & Keil, 2013; Anderson & Müller, 2010; Regan, 1975). Does it matter? Yes, at least when there is considerable amplitude at the higher harmonics, relative to F. In this case, higher harmonics do contribute significantly to the response measurement (for a dramatic illustration to this effect, see Figure 1).

How often is there considerable amplitude at the higher harmonics? Can the cost of omitting higher harmonics in published studies be evaluated? Unfortunately, most studies do not report whether or not there were responses at higher harmonics, as mentioned above. Moreover, when higher harmonic responses were reported to be present, they were often not described, for example, “Peaks were also present at the harmonics of the stimulus frequency but were not analyzed in this study” (Srinivasan et al., 1999, p. 5438); “Higher harmonics may play a role, especially at lower temporal frequencies (see, for example, the double peaks in the 12-Hz data in Figure 1), but these are not considered here” (Kremers, Rodrigues, Silveira, & da Silva Filho, 2010, p. 579); “Note however that 2 Hz is a harmonic of 1 Hz and may actually be a relevant spectral region to consider (albeit outside the scope of this report)” (Kosem, Gramfort, & Van Wassenhove, 2014, Figure S2; evident in figures but not discussed or included in analyses: Schettino, Porcu, Gundlach, Keitel, & Müller, 2020; Eidelman-Rothman et al., 2019; Chadnova et al., 2018; Winawer et al., 2013; Hönegger et al., 2011; Katzner et al., 2009; Pastor, Artieda, Arbizu, Valencia, & Masdeu, 2003; Pastor et al., 2002; Kaspar, Hassler, Martens, Trujillo-Barreto, & Gruber, 2001; Müller et al., 1998), preventing a wide-scale review.

From some studies that have described the presence or absence of higher harmonics, there appear to be a couple, specific cases in which their amplitude was not considerable and may be neglected with little cost. The most well-documented case is that in which high stimulus presentation rates are used (as will be demonstrated in the following section; with electrophysiology: Retter & Rossion, 2016a; Heinrich, Groten, & Bach, 2015; Alonso-Prieto, Van Belle, Liu-Shuang, Norcia, & Rossion, 2013; Capilla, Pazo-Alvarez, Darriba, Campo, & Gross, 2011; Tlumak, Durrant, Delgado, & Robert Boston, 2011; Vialatte, Maurice, Dauwels, & Cichocki, 2009; Luck, 2005; Kremers & Scholl, 2001; Ross, Borgmann, Draganova, Roberts, & Pantev, 2000; Regan, 1989; Van der Tweel & Verduyn Lunel, 1965; with neuroimaging: Retter, Webster, & Jiang, 2019; Gao, Gentile, & Rossion, 2018; Puce, Allison, Gore, & McCarthy, 1995). However, the case of high stimulus presentation rates cannot be readily identified across studies. This is because “high” is dependent on the relationship of the stimulus presentation rate to the duration of the brain responses being measured, by the recording technique, in the assessed population (Retter, Jiang, Webster, & Rossion, 2020; Keysers & Perrett, 2002; see also Heinrich, 2010).

A second case is that in which low-amplitude first harmonic responses were reported, such as responses elicited with subtle stimuli (e.g., with no higher harmonic amplitude above noise: Lochy, Schiltz, & Rossion, 2020; Park, 2018; Retter & Rossion, 2017; Brazier, 1964; with very low higher harmonic amplitudes: Moungou, Thonnard, & Mouraux, 2016; McFadden et al., 2014; Ales, Farzin, Rossion, & Norcia, 2012). However, this case also does not allow general inferences: low first harmonic amplitudes may still coincide with large higher harmonic amplitudes (see the Frequency-tagged Responses in the Frequency Domain section; also, e.g., Gaume, Vialatte, & Dreyfus, 2014; Alonso-Prieto et al., 2013; Capilla et al., 2011; Vialatte et al., 2009), such that the amplitude of the first harmonic itself is not diagnostic.

On the other hand, considerable amplitude at higher harmonics has been reported in a wide array of studies. For example, higher harmonics often exceed the fundamental in studies on the brain's responses for an extensive range of processes (with F usually below 8 Hz, recorded to visual stimuli with EEG/MEG): from luminance (patterns: Gaume et al., 2014; Capilla et al., 2011; Vialatte et al., 2009; Davila, Srebro, & Ghaleb, 1998; electrocorticography: Winawer et al., 2013; electroretinogram: Baker & Hess, 1984) to color and motion (McKeefry et al., 1999; Tyler & Kaitz, 1977) to face perception (Liu-Shuang, Norcia, & Rossion, 2014; Alonso-Prieto et al., 2013; below about 40 Hz to auditory stimuli: Ross et al., 2000; see also Tlumak et al., 2011). Higher harmonics may also be present, each with a lower amplitude than the fundamental, but with their amplitude distributed across a large range of harmonic frequencies (e.g., Cunningham, Baker, & Pierce, 2017; Retter & Rossion, 2016a; Gaume et al., 2014; Painter, Dux, Travis, & Mattingley, 2014; Alonso-Prieto et al., 2013; Capilla et al., 2011; Vialatte et al., 2009). Considerable higher harmonics have also been demonstrated with low-temporal resolution techniques, including fMRI, given appropriately slow stimulus presentation frequencies (e.g., with F well below 0.1 Hz; motor activity: Bandettini, Jesmanowicz, Wong, & Hyde, 1993; luminance patterns: Engel et al., 1997).

In some cases, the conclusions of studies considering and not considering higher harmonics can be compared. For example, in frequency-tuning studies, without considering higher harmonics, maximal visual responses were reported to stimuli modulated at about 10–15 Hz with EEG (e.g., Ding, Sperling, & Srinivasan, 2006; Pastor et al., 2003; Regan, 1966; see also Vialatte et al., 2009). However, when higher harmonics were considered, the lowest stimulation frequency tested (3 Hz) yielded the maximal visual EEG responses, being over 3 times higher than the responses to 12 Hz stimulation with natural images (Retter et al., 2020, Figure S3). Similarly, maximal auditory responses were reported to stimuli modulated at about 40 Hz with EEG (e.g., Pastor et al., 2002; Ross et al., 2000; Galambos, Makeig, & Talmachoff, 1981), but a consideration of higher harmonics produced the maximal auditory EEG responses at the lowest stimulation frequency tested (0.75 Hz; Tlumak et al., 2011, Figure 4). Overall, although it is thus impossible to ascertain what the impact of unreported or uncharacterized higher harmonic responses in most studies may have been, it is likely that it was often considerable.

At present, harmonics are surrounded by many questions: Why do they occur? What do they represent? Which, or how many, harmonics should be considered? Should they be taken into account for response identification and measurement, and if so, how? Indeed, the lack of understanding and standard practice regarding higher harmonics has limited the ease (i.e., objectivity) of frequency-tagged response identification and measurement. This is particularly significant because objectivity is given as a primary advantage of the frequency-tagging technique, contributing to its increasing application in (cognitive) neuroscience research and clinical applications (e.g., see Rossion et al., 2020; Norcia et al., 2015).

In previous studies that reported harmonic responses, most often these responses have been described individually (e.g., at F, 2F, 3F) and have not been taken into account for response measurement (e.g., Cunningham et al., 2017; Moungou et al., 2016; Painter et al., 2014; Alonso-Prieto et al., 2013; Ales et al., 2012; Capilla et al., 2011; Rossion & Boremanse, 2011; Vialatte et al., 2009; Herrmann, 2001; Ross et al., 2000; Srinivasan et al., 1999; Bandettini et al., 1993; Tyler & Kaitz, 1977). Although considering harmonic responses separately is considerably better than not at all, individual harmonic responses do not represent independent aspects of a time domain response (e.g., see Tang & Norcia, 1995; for dependent harmonic amplitude examples: Retter & Rossion, 2016a; Zhou et al., 2016; for qualitatively similar neighboring harmonic examples: Rossion et al., 2020; Zemon & Gordon, 2018; Jacques, Retter, & Rossion, 2016; Rossion, 2014; see the Interpreting Harmonics section).

In practice, considering higher harmonic responses improves response detection, measurement, and classification (e.g., Zemon & Gordon, 2018; Retter & Rossion, 2016a; Tlumak et al., 2011; Cebulla, Stürzebecher, & Elberling, 2006; Davila et al., 1998; for brain–computer/machine interfaces: Cetin, Ozekes, & Varol, 2020; Chen et al., 2015; Muller-Putz, Scherer, Brauneis, & Pfurtscheller, 2005). Combining harmonics is particularly useful for comparing response amplitudes across experimental conditions. Otherwise, if one input produces a response with larger amplitudes than another input at some harmonics but not others, how could these responses be evaluated overall? Or, how could the relative change (e.g., percent increase) of one response relative to another be calculated overall? In some previous studies, harmonic responses have been combined with various approaches, such as root-mean-square summation (that is, the square root of the summed squared amplitudes, also known as the “root sum square”, or “summation in quadrature”; e.g., Tlumak et al., 2011; Appelbaum, Ales, Cottereau, & Norcia, 2010; Appelbaum, Wade, Vildavski, Pettet, & Norcia, 2006; Hou, Pettet, Sampath, Candy, & Norcia, 2003), (weighted) summation (e.g., Zhang, Jamison, Engel, He, & He, 2011; Wang, Gao, Hong, Jia, & Gao, 2008; Cheng, Gao, Gao, & Xu, 2002), or averaging (e.g., Milton et al., 2020; Lochy, Van Belle, & Rossion, 2015; Liu-Shuang et al., 2014). However, these approaches have not been justified and have not been related to approaches analyzing the brain's response amplitudes in the time domain (or physiologically: see Heinrich, 2010).

In the following, a validated methodology for combining (baseline-corrected) harmonic amplitudes through simple summation will be provided. This approach derives from a theoretical basis of how signals over time are represented through mathematical transformations into the frequency domain (see the Frequency Domain Representations section), extended to experimental responses in practice (see the Frequency-tagged Responses in the Frequency Domain section). This approach was indicated empirically by Retter and Rossion (2016a), and it has since been applied in a number of studies, however, primarily by those authors or associated research groups (e.g., van de Walle de Ghelcke, Rossion, Schiltz, & Lochy, 2021; Damon, Leleu, Rekow, Foncet, & Baudouin, 2020; Fisher, Towler, Rossion, & Eimer, 2020; Gwinn & Jiang, 2019; Dwyer, Xu, & Tanaka, 2019; Van der Donck et al., 2019; Beck, Rossion, & Samson, 2018; Chemin, Huang, Mulders, & Mouraux, 2018; De Keyser, Mouraux, Quek, Torta, & Legrain, 2018; Guillaume, Mejias, Rossion, Dzhelyova, & Schiltz, 2018; Gwinn, Matera, O'Neil, & Webster, 2018; Leleu et al., 2018; Xu, Liu-Shuang, Rossion, & Tanaka, 2017). To be of further use to the scientific community, the approach requires deeper methodological evaluation and, especially, evaluation in a theoretical context, which is the goal throughout this article. From this, some practical guidelines are offered (see the Combining Harmonic Responses section) and implications are drawn for the interpretation of harmonic responses more generally (see the Interpreting Harmonics section).

When a signal is transformed into the frequency domain (by means of a Fourier transform), it becomes represented through a combination of sine waves, which are the fundamental units of the frequency domain. Although there are many texts on the mathematics of frequency transformations and representations (e.g., Gonzalez & Woods, 2018; Forinash & Christian, 2016; Patel, 2012; Strang, 2007; Smith, 1997; Press, Falnnery, & Teukolsky, 1993), a basic understanding of sine waves and their combination is a sufficient foundation for the interpretation of multiharmonic responses of the brain (Regan, 1989).

Briefly, sine waves are trigonometric functions that describe periodic signals in terms of frequency, amplitude, and phase (Figure 2). The frequency of a sine wave describes the number of cycles (of 360°; equivalent to 2π radians) per unit of time or space (time is typically given in units of cycles/sec = sec⁻¹ = Hertz = Hz). Note that the cycles of sine waves are periodic and could repeat their pattern infinitely, as a circle could be endlessly traced (the sine wave, as in Figure 2B, derives from the y-axis values of a unit circle, as shown in Figure 2A). The amplitude of a sine wave is defined along the y-axis (e.g., in Figure 2B, the sine wave, spanning from −1 to 1, has an amplitude of 1; its unit in EEG recordings is typically microvolts). The phase of a sine wave is a measure of its starting angle (indicated by theta in Figure 2A), with an arbitrary beginning at zero (in units of degrees or radians), as is shown in Figure 2B. Changes in frequency, amplitude, and phase, which help demonstrate these properties, are illustrated in Figure 2C. For those whom it helps to see it mathematically, the expression of a sine wave, as a function of x, is: y(x) = asin(2πfx + ϕ), where a = amplitude (scaling on the y-axis), f = frequency (by cycles), and ϕ = phase (x-axis shifts).

A frequency domain representation of a signal is essentially a lot of sine waves. That is, when a signal is transformed into the frequency domain, the resultant x-axis describes the frequency of its constituent sine waves. The other descriptors of sine waves, amplitude, and phase are described in the transformed, complex-valued y-axis at each frequency, which is typically plotted as separate amplitude and/or phase frequency spectra.² The resolution (x-axis sampling) of the frequency domain spectrum is the inverse of the signal recording length and the range spans from zero to half of the signal sampling rate (note that these properties have practical implications for frequency-tagging experimental design, e.g., as addressed in Bach & Meigen, 1999). The combination, through summation, of these sine waves described in the frequency domain reconstructs the original signal in the time domain. Here, the focus will be on periodic signals over time, but note that frequency domain analyses can be applied in many settings (e.g., signals over space or over two dimensions).

In the simplest case, a periodic signal that is a perfect sine wave is represented in the frequency domain by a single frequency, representing a single sine wave (Figure 3A). Another way to understand this is to observe that, in this case, the frequency, amplitude, and phase of a single sine wave in the frequency domain are sufficient to reconstruct the original signal in the time domain. In most cases, signals are more complex (i.e., nonsinusoidal), but this does not pose a problem: A combination of sine waves at different frequencies can sum to model any signal. In a classic example, a periodic squarewave signal is shown to be represented with a sum of sine waves specific to its periodicity (Figure 3B). Nonperiodic signals, for example, ERPs to temporally jittered stimuli, can also be represented in the frequency domain, but because they are not specific to limited frequencies, their interpretation does not correspond to that of frequency-tagged signals (Figure 3C). Although frequency-based analyses of nonperiodic signals may be applied (e.g., see Herrmann, Rach, Vosskuhl, & Strüber, 2014; Başar & Schürmann, 1994; Regan, 1989; see also Chemin et al., 2018), these are outside the present focus on frequency tagging.

It may be observed that a simple sine wave ranging from −1 to 1 in the time domain has an amplitude of 1 in the frequency domain, but that the relationship between the time domain and frequency domain amplitudes for multiharmonic signals is more complex (compare Figure 3A and B). However, there is a direct relationship between these dimensions: Because energy is conserved from the time to the frequency domain, the sum of the root-mean-square amplitudes of the time domain signal equals the sum of the squared root mean amplitudes of the frequency domain signal (Parseval's relation; Parseval des Chênes, 1806; e.g., see Smith, 1997).³ For example, the sum of the squared amplitudes per cycle of the time domain signal in Figure 3A is equal to 0.5, and the sum of the squared root-mean-square amplitudes of its discrete frequency domain signal is equal to 0.5. Multiharmonic signals also preserve this relationship, although their time domain amplitude range does not directly relate to their frequency domain amplitude (being affected by phase; see the What About Phase? section).

Frequency tagging is an approach in which stimuli are presented periodically to generate periodic responses of the brain that can thus be identified in the frequency domain at specific frequencies harmonic to the stimulation (i.e., the fundamental and the higher harmonic frequencies). In the following, examples will be taken only for harmonics that are specific to a single tagged frequency: Again, for determining specific harmonics in the context of multiple tagged frequencies, see the Which Harmonics to Consider? section).

According to the principles of frequency analyses, a simple, sinusoidal brain response would be represented only at the fundamental frequency F, whereas more complex brain responses would be represented with a combination of F and its higher harmonics, 2F, 3F, and so on. Generally, this is evidenced with experimental brain responses. In the event that the brain responses are nearly sinusoidal, the response is dominated by amplitude at F (e.g., at a high stimulus presentation rate: Figure 4A); in the event that the responses of the brain are complex, a combination of sine waves at different frequencies (i.e., the higher harmonics) can sum to model any signal (Figure 4C). In many studies, complex, nonsinusoidal responses of the brain evoked over time are represented in the frequency domain not only at F, but with considerable amplitude at its higher harmonics (as addressed in the Introduction, e.g., Rossion et al., 2020; Retter & Rossion, 2016a; Zhou et al., 2016; Norcia et al., 2015; Alonso-Prieto et al., 2013; Tlumak et al., 2011; Heinrich, 2010; Vialatte et al., 2009, 2010; Kremers & Scholl, 2001; Ross et al., 2000; Bach & Meigen, 1999; Sieving, Arnold, Jamison, Liepa, & Coats, 1998; Regan, 1966, 1989; Donker, 1975; Van der Tweel & Verduyn Lunel, 1965; Brazier, 1964).

Indeed, higher harmonic responses may be accounted for in relation to the complex (i.e., nonsinusoidal) responses of the brain, in accordance with the principles of frequency domain analysis of periodic signals (as in Rossion et al., 2020; Zhou et al., 2016; Norcia et al., 2015; Heinrich, 2010; Regan, 1989). This account explains that higher harmonics are present when complex brain responses are present, but does not implicate a specific source of complex brain responses (see Heinrich, 2010). However, it is important to note that complex brain responses are not a product of frequency tagging and may equivalently occur with nonperiodic (event-related) stimulus presentation modes.⁴

The harmonics do not represent new information, specific to the frequency domain: They are merely highlighted in an alternative, frequency domain representation of the original time domain signal (certain variations may be represented more or less clearly in each domain). These domains are interchangeable: As time domain data can be transformed into the frequency domain, frequency domain data can also be inversely transformed back into the time domain; indeed, time domain responses can be reconstructed from frequency domain harmonic amplitudes and phases (e.g., Ruhnau, Keitel, Lithari, Weisz, & Neuling, 2016; Sieving et al., 1998). The next section focuses on the combination of harmonic responses; however, individual harmonics will be further addressed in the Interpreting Harmonics section.

In theoretical examples of signal transformation into the frequency domain (as in the Frequency Domain Representations section), the signal is pure signal. In frequency tagging, as in all brain recordings, the signal (here, i.e., the responses of the brain at the tagged frequencies) also carries “noise,” a term that refers to both non-event-related brain activity and artifacts (e.g., Luck, 2005; Regan, 1989). Note that there are many methods for correcting for noise (i.e., isolating signal) in frequency-tagging research, although a discussion of these is more general than the scope of this article (see instead, e.g., Meigen & Bach, 2000; Appendix 2 of Norcia et al., 2015). In the examples given here, a simple correction for noise will be applied that subtracts a local baseline from the amplitude of the frequencies of interest (e.g., Retter & Rossion, 2016a).⁵ The baseline is defined as the mean amplitude of a symmetrical range of neighboring (i.e., continuously adjacent⁶) frequency bins (for theoretical justification, see Norcia et al., 2015; Regan, 1989; e.g., Boremanse, Norcia, & Rossion, 2013; Meigen & Bach, 2000; Peterzell & Norcia, 1997; with power rather than amplitude: Mouraux et al., 2011; Vialatte et al., 2009; Srinivasan et al., 1999). This method is used to provide a measure of signal amplitude in the frequency domain that is relatable (i.e., both in the same unit) to amplitude in the time domain and with a noise level at zero, while compensating for local variations of noise inherent to brain recordings across the frequency spectrum.

The combination of sine waves is simple: Sine waves sum linearly to reconstruct a signal. However, with the goal of identifying and measuring overall response amplitude in the frequency domain, because sine waves carry both amplitude and phase information, their sum is not intuitive to interpret in terms of amplitude only (or phase only). Note that there are alternative approaches for combining harmonics that incorporate both amplitude and phase; however, these approaches make use of phase as an indicator of reliability (coherence), typically across short stimulation durations (e.g., Delorme & Makeig, 2004; Strasburger, 1987; Jervis, Nichols, Johnson, Allen, & Hudson, 1983). Furthermore, although the amplitude of a single sine wave in the frequency domain directly relates to its time domain peak amplitude, that is, half its positive to negative peak range (e.g., Figure 3A; similarly, see Figure 4A), the frequency domain amplitude of complex time domain signals (summed sine waves) is not as easily visualized from the time domain. Perhaps for these reasons, various approaches have been taken for the combination of frequency-tagged multiharmonic brain response amplitude in the frequency domain (e.g., as mentioned previously, averaging or root-mean-square summing, respectively, Liu-Shuang et al., 2014; Hou et al., 2003).

The summation of harmonic amplitudes is recommended here for identifying and measuring the overall brain response (based on Retter & Rossion, 2016a; see also Heinrich, 2009). In the study of Retter and Rossion (2016a), this approach was validated empirically by qualitative comparison of time and frequency domain responses, in the situation where several equivalent time domain EEG responses were produced by several slow target stimulus presentation frequencies. There, it was observed that despite different distributions of harmonic amplitudes stemming from the different fundamental target stimulus presentation frequencies (1.1–2.5 Hz), the summation of baseline-subtracted harmonic amplitude across a common frequency range led to equivalent overall amplitudes that related to approximately equivalent response amplitude peaks in the time domain by visual inspection (see Figure 5A and B, Row 1, here, for examples of reprocessing of that data in combination with the underlying harmonic distributions). Moreover, a faster stimulus presentation rate (4.2 Hz), which produced visually lower amplitude deflections in the time domain (Figure 5C, Row 1), also produced a lower summed harmonic response amplitude.

Here, these data are revisited quantitatively, with typical time domain interpretations of response amplitudes: peak-to-peak amplitude of the largest deflections and the area under the curve of the response deflections (Figure 5, Row 4). Note that an exact comparison of specific deflections, as is more commonly done in relating brain responses, is possible for the conditions at 1.1 and 1.4 Hz, but that a different response pattern is observed for the condition at 4.2 Hz, preventing such a direct comparison.

The summed harmonic response amplitude is shown to be congruent with these measures (Figure 5D). Critically, other approaches for harmonic combination would not have led to these conclusions when comparing conditions. Here, in the frequency domain, a large fundamental harmonic amplitude relates to a smaller number of harmonic responses (with an amplitude above 0.1 μV). Therefore, for example, averaging the harmonic responses would have generated lower amplitude responses the slower the stimulus presentation frequency (1.1 Hz < 1.4 Hz < 4.2 Hz; Figure 5D). For another example, the root-mean-square harmonic amplitude would have generated the highest amplitude for the highest stimulus presentation rate (4.2 Hz; Figure 5D). For a last example, using non-baseline-corrected amplitudes would have produced a larger response at 1.1 Hz than 1.4 Hz, because “noise” would have been included at more and lower frequency (noisier) harmonics at 1.1 Hz.

Thus, summing baseline-subtracted harmonic amplitudes is advantageous for a correspondence with interpretations of time domain brain responses. This approach has been used to quantify and compare overall response amplitude in a number of studies following Retter and Rossion (2016a), as mentioned previously (e.g., including time domain correspondences: De Keyser et al., 2018; Leleu et al., 2018; frequency domain analyses only: Damon et al., 2020; Dwyer et al., 2019; Gwinn & Jiang, 2019; Beck et al., 2018; Chemin et al., 2018; Guillaume et al., 2018; Gwinn et al., 2018; Xu et al., 2017).⁷

If measures other than amplitude are desired (e.g., signal-to-noise ratio, z scores, or another statistic), the harmonic amplitudes can be extracted with an inclusion of a baseline frequency range (i.e., as a “chunk” of X Hz, centered around each frequency of interest), and then summed before these baseline-relative computations (Retter & Rossion, 2016a; see also Appendix 2 of Norcia et al., 2015; Box 2 of Rossion et al., 2020). In this way, a single statistical measure can be applied to the combined harmonic amplitude relative to its combined baseline amplitude (i.e., “noise”). Note that different approaches for combining harmonics may serve different ends, in that they describe different aspects of the signal, for example, the root-mean-square amplitude relates to the equivalent power of a flat (nonsinusoidal) signal; however, these aspects must be justified in relation to their physiological meaning.

As addressed previously, there is a direct relationship between signal amplitude in the time and frequency domains. As a reminder, this relationship is given by Parseval's relation, which states that energy is conserved across the time domain (where energy equals the sum of the squared amplitudes) and frequency domain (where energy equals the sum of the squared root-mean-square amplitudes). In light of this, the amplitude across the harmonics relates to the overall amplitude of the signal in the time domain, regardless of phase.

However, to fully relate signals across the time and frequency domains, both the amplitude and phase of the representative frequency domain sine waves need to be taken into account. Without phase information, the fluctuation of amplitude across time (e.g., affecting local amplitude peaks) cannot be determined. Therefore, there is a cost toward relating time and frequency domain signals when excluding phase information. However, this cost is reasonably minor, for example, as relative phase changes across harmonics, it is possible that the latency of signal peaks varies, but that their amplitude does not (Figure 6A and B). When relative phase does affect peak amplitudes, this influence is limited (e.g., compare Figure 6B and C). Moreover, despite relative phase changes, the area under the curve of the time domain signal may remain approximately constant (Figure 6A–C; see also Heinrich, 2010).⁸

Moreover, in frequency tagging, it is worth remembering that the phase is not arbitrary: The phase of each relevant harmonic is determined relative to the time domain signal. In other words, the aligning positive and negative peaks of the sine waves across harmonic frequencies correspond to the time of the positive and negative peaks of the signal in the time domain (see again Figure 3B, for an example of deconstructive [when the signal is 0] and constructive [when the signal is 1], phase-locked harmonic sine wave superpositioning). This leads to phase differences across harmonics that are similar to describe time domain signals with similar temporal dynamics, despite the use of different stimulation frequencies (Figure 6D; see also Strasburger, 1987). Thus, the influence of phase on combined harmonics is largely invariant of the stimulus presentation frequency, given consistent temporal dynamics of the response (as hinted at empirically, e.g., Retter & Rossion, 2016a; Appelbaum et al., 2006). Finally, it is worth noting a couple of helpful restrictions in the context of frequency tagging for combining harmonic responses: Only one sine wave is represented at each frequency bin, and nonharmonic frequencies are not considered, such that the response is fully periodic at the cycle duration of the fundamental frequency, F.⁹

Before combining harmonics, a decision of which harmonics to consider is required. To this extent, harmonics of interest (similarly to a ROI) must be defined. The first criterion for determining harmonics of interest is whether a higher harmonic is specific to its fundamental frequency. As mentioned previously, in frequency-tagging paradigms using a single-stimulation frequency, the higher harmonics are always specific to the fundamental frequency. However, in paradigms deploying multiple stimulation frequencies, unspecific harmonics may occur, which are often excluded from the analyses (for a paradigm-focused review, see Norcia et al., 2015). For example, two stimuli may be simultaneously presented at different spatial locations, one at 8 Hz (F₁) and the other at 6 Hz (F₂). If a response occurred at 24 Hz, it would not be specific to either stimulus, being the third harmonic of 8 Hz (3F₁) and the fourth harmonic of 6 Hz (4F₂), and would therefore be excluded from the analyses of responses to each stimulus (for further examples, see Table 1).

A second consideration for determining harmonics of interest aims to exclude extreme harmonic frequencies (e.g., the 30th harmonic of 8 Hz, at 240 Hz) at which no signal is expected or found. A limited selection of harmonics has been made based on various types of criteria: (1) amplitude, power, signal-to-noise ratio, or significance thresholds (e.g., Rossion, Torfs, Jacques, & Liu-Shuang, 2015; Hou, Pettet, & Norcia, 2008; Donker, 1975); (2) frequency range (e.g., Leleu et al., 2018; Zemon & Gordon, 2018; Jacques et al., 2016; Sieving et al., 1998); (3) harmonic series number (e.g., Wittevrongel et al., 2018; Appelbaum et al., 2006; Donker, 1975); (4) in relation to other stimulation frequencies (e.g., Milton et al., 2020; Heinrich, Mell, & Bach, 2009); and (5) correlation with the time domain response (e.g., Engel et al., 1997; Bandettini et al., 1993).

The use of a limiting frequency range is recommended here, either as determined a priori or from an assessment of the highest harmonic meeting a threshold (in terms of amplitude, signal-to-noise ratio, or significance). This is recommended because the upper frequency limit of harmonic responses, although affected by the overall strength of the signal, generally relates to the highest frequency that is strongly represented in the signal (see the Interpreting Harmonics section). The upper frequency limit of harmonic responses is thus often conserved across fundamental stimulation frequencies (see Figure 2 of Retter & Rossion, 2016a).

Note that the highest harmonic of interest can be determined either at the group level across conditions or as presented by any participant for any condition, but that typically a common range of frequencies of interest should be used across participants and conditions (e.g., Jacques et al., 2016). In this approach, there may be harmonic frequencies included for consideration at which there is no signal (e.g., in some participants, conditions, or regions of interest); however, including a small number of such frequencies is likely less detrimental (given that an appropriate baseline noise correction is applied, e.g., so that approximately zero amplitude values are added) than missing some frequencies containing a weak signal. Similarly, although responses are typically expected to occur consecutively across harmonic frequencies, in the event that a small number of within-range harmonic frequencies do not contain signal (above threshold), including them is typically tolerable (e.g., Rossion et al., 2020; Liu-Shuang et al., 2014).

Finally, in some cases, harmonic responses appear to be qualitatively different from one another. This may occasionally be related to physiological sources: for example, different harmonic response patterns are generated from the recordings of frequency-tagged responses from single- versus double-opponent cortical cells (Movshon et al., 1978). However, more often, physiological sources may only be tentatively inferred, for example, when different EEG scalp topographies are observed at different harmonic frequency ranges (e.g., Rossion, 2014; see the What Do Higher Harmonic Responses Represent? section). In this case, is it appropriate to select subranges of qualitatively homogeneous harmonics to consider and/or combine? Perhaps, although it should be remembered that harmonic responses are not independent of one another (e.g., Retter & Rossion, 2016a; Zhou et al., 2016) and therefore should also be described individually and/or summed all together (see the Should Higher Harmonic Responses Be Combined, and If So, How? (Reprise) section).

It is not advised to select or subgroup harmonics a priori in accordance to only their number, unless this is explicitly derived from the stimulation paradigm (see Table 1). For example, there is a persistent history of considering the first versus second harmonic response (Saupe, Schröger, Andersen, & Müller, 2009; Pastor, Valencia, Artieda, Alegre, & Masdeu, 2007; Kremers & Scholl, 2001; Falsini et al., 1999; Burns, Elsner, & Kreitz, 1992; Baker & Hess, 1984; the first vs. second harmonic, rather than odd vs. even harmonics: Kim, Grabowecky, Paller, & Suzuki, 2011; Kim, Grabowecky, Paller, Muthu, & Suzuki, 2007). This relates to early interpretations of the first harmonic reflecting asymmetries in responses following on and off stimulation cycles (e.g., Clynes, Kohn, & Lifshitz, 1964) and the second harmonic being typically dominant with pattern reversal stimulation (see the following section). However, the presence of third, fourth, and further higher harmonics—and their interdependence—is indicative of the limits of such an oversimplification in most cases.

Higher harmonic responses represent complex neural responses in the time domain. At a fundamental level, these harmonic responses are like any other frequency domain representations: They are sine waves described by frequency, amplitude, and phase (Figure 2). Although only one sine wave is required to describe a sinusoidal signal in the time domain, a combination of (a lot of) sine waves is required to describe complex signals in the time domain (Figure 3). Frequency-tagged brain responses are periodic in the time domain, and thus, only sine waves periodic to their fundamental frequency (i.e., the harmonics) are mathematically available to describe them. Simple brain responses require few harmonics, whereas complex responses require more harmonics (Figure 4). For example, lower stimulation frequency responses often have more harmonics, because there are relatively more harmonic frequencies available within a relevant frequency range ceiling (Figure 5).

Higher harmonic responses have often been interpreted as being caused by nonlinearities in the stimulus presentation and/or brain responses (e.g., Gordon, Hohwy, Davidson, van Boxtel, & Tsuchiya, 2019; Norcia et al., 2015; Burns et al., 1992; Regan, 1966, 1989; Troelstra, 1971; Spekreijse, 1969; Van der Tweel & Verduyn Lunel, 1964, 1965; see also Shapley, 2009). This relates to early attempts to present stimuli perfectly sinusoidally (e.g., a uniform field being modulated sinusoidally in luminance; since van der Tweel et al., 1958). The rationale was that if the brain's response to a sinusoidal stimulus was linear at the level of recording, it too would be perfectly sinusoidal in following this stimulus and would be represented in the frequency domain by a response only at the fundamental frequency, that is, without higher harmonics (see the Frequency Domain Representations section). Contrarily, higher harmonic brain responses were produced in most cases and were attributed to nonlinearities in the brain's responses themselves (e.g., nonlinear action potential firing, neural population response dynamics; Shapley, 2009; Skottun et al., 1991; Movshon et al., 1978).¹⁰

However, more recent studies suggest that the amount of nonlinearity, or complex temporal frequency content, in stimulus presentation may not correspond with the amount of nonlinearity in the brain's response at the population level: There was little to no difference in the higher harmonic amplitude distributions in response to (imperfect) sinusoidal versus squarewave (i.e., abrupt on/off) complex stimulus presentation (Dzhelyova, Jacques, & Rossion, 2017; Retter, 2016; Fawcett, Barnes, Hillebrand, & Singh, 2004; Burns et al., 1992). Moreover, although the inherent nonlinearity of the brain's responses could account for higher harmonics, in practice, the amplitude of the higher harmonics is not always above noise level (as addressed in the Introduction) or is very low, suggesting only a modest contribution of this factor (see the Frequency-tagged Responses In The Frequency Domain section). That being said, one source of complexity in the brain's responses (even those underlying the first harmonic) may be these nonlinearities.

Higher harmonic responses represent the relevant frequency characteristics of the response in the time domain (e.g., Rossion et al., 2020; Zemon & Gordon, 2018; Retter & Rossion, 2016a; Zhou et al., 2016; Gaume et al., 2014; Heinrich, 2010; Norcia, Sato, Shinn, & Mertus, 1986; Baker & Hess, 1984; Galambos et al., 1981). That is, dynamics of the time domain response best represented at different frequency ranges will produce more amplitude in those frequency ranges in the frequency domain (Figure 7; compare with the frequency domain representation in Figure 3C; for another example of frequency representations over time, see Makeig et al., 2002). To visualize the impact of (a range of) individual harmonics, partial harmonic time domain reconstructions have been plotted (Sieving et al., 1998; Baker & Hess, 1984). Note that the amplitude distribution of harmonics across frequencies is, however, affected by the fundamental stimulation frequency and the overall amplitude of the signal. Moreover, individual harmonics do not represent independent or temporally separated aspects of a time domain response (e.g., Retter & Rossion, 2016a; Zhou et al., 2016).

In line with the above interpretation, harmonic responses may be (gradually) influenced, quantitatively and/or qualitatively, by the frequency at which they fall. For example, Retter and Rossion (2016a, Figure 2C) described harmonic EEG responses that were gradually characterized by frequency, in terms of amplitude and scalp lateralization, commonly across conditions with different, low stimulation frequencies. Harmonic responses may also appear to group into somewhat distinct frequency ranges, for example, visually evoked EEG responses above about 10 Hz having a more medial-occipital (“low-level”) scalp topography (Jacques et al., 2016; Rossion, 2014; see also Zemon & Gordon, 2018). In general, higher harmonics may be more associated with earlier response dynamics, as the onset slopes of (ERP) responses are typically steeper than the offsets (Norcia et al., 1986; as in Figure 7). In some cases, different experimental effects may be pronounced at some harmonic frequencies and not others (e.g., an increased response because of attention: Saupe et al., 2009; Pei, Pettet, & Norcia, 2002; opposing effects following current stimulation: Ruhnau et al., 2016). Indeed, limited frequency ranges may be appropriate for measuring certain properties (e.g., a lower frequency range for chromatic than luminance signals; see Parry et al., 2012; Burns et al., 1992). Thus, in addition to combining harmonics, investigating individual (ranges of) harmonics may provide insight into the functional dynamics of the neural processes that occur at their respective frequencies.

Differing harmonic response distributions may stem from the same fundamental stimulation frequency, such as are observed across conditions (along with time domain response dynamic differences: e.g., Dzhelyova et al., 2017; Jacques et al., 2016), populations (e.g., clinical: Van der Donck et al., 2019; Falsini et al., 1999; age/grade levels: van de Walle de Ghelcke et al., 2021; Hou et al., 2003), and individuals (as observed in Heinrich et al., 2015; Tlumak et al., 2011; Heinrich & Bach, 2001; used in classifying individuals across testing sessions: Dzhelyova et al., 2019). In past studies, individual harmonic amplitudes have been displayed with frequency domain spectra or tables, sometimes with additional descriptors, such as functional images or topographies (e.g., Jacques et al., 2016; Liu-Shuang et al., 2014; Bandettini et al., 1993). In more global visualizations, colored/shaded matrices of harmonic responses by ROIs or participant have been plotted (Liu-Shuang et al., 2014, Figure 3; Rossion et al., 2015, Figure 3), as well as close-shaped “fingerprint” plots, using connected vectors of which the angle is determined by harmonic frequency and length by amplitude (van de Walle de Ghelcke et al., 2021, Figure 6).

It is not recommended to display harmonics independently by their sequential number (e.g., across different fundamental stimulation frequencies: Herrmann, 2001; Ross et al., 2000; Troelstra, 1971). Indeed, because harmonic responses are characterized by their frequency, they are generally not well characterized by sequential number, irrespective of frequency, unless differentially tagged in the stimulation paradigm (see the Which Harmonics to Consider? section). However, one case in which harmonic frequency does not characterize responses well is at high stimulus presentation frequencies, for which the first harmonic is dominant. At high frequencies, overlapping responses interfere with one another, decreasing the measured response's complexity and amplitude and therefore harmonic contents (see the Frequency-tagged Responses in the Frequency Domain section). Thus, a dominant exponential decrease of harmonic amplitude may be observed across sequential harmonics (compare Figure 4B and C; Retter et al., 2020, Figure S2). However, this must not be taken as evidence of distinct brain responses being evidenced at the first versus second (vs. third, etc.) harmonic (e.g., see the discussion of Saupe et al., 2009).

Finally, it may be noted that consecutive response interference is not only dependent upon the duration of the brain response being measured but is also influenced by the spatiotemporal dynamics of the interfering response components. In models of consecutive responses, destructive superpositioning of overlapping (ERP) response components may reduce or even eliminate (selective harmonic) responses (constructive superpositioning is also possible; see Heinrich, 2010, Figures 1 and 2; see also Footnote ⁴). In practice, overlapping brain responses do not occur independently of one another (see Retter et al., 2020; Keysers & Perrett, 2002). However, particularly when distinct response sources are implicated (e.g., light- and dark-preferring neural populations or across sensory modalities), differential interactions of stimulation frequency and selective harmonic amplitudes could also provide insights into the underlying dynamics of the brain responses.

Limiting the analysis of frequency-tagged responses to a (nonpredominant) fundamental harmonic is not recommended: Taking harmonic responses into account leads to substantially improved response quantification and classification (Figure 8; Cetin et al., 2020; Zemon & Gordon, 2018; Retter & Rossion, 2016a; Chen et al., 2015; Tlumak et al., 2011; Cebulla et al., 2006; Muller-Putz et al., 2005; Davila et al., 1998). The combination of harmonic amplitudes through summation is justified through the principles of frequency-based analyses, and it leads to a combined response measurement that relates to typical time domain amplitude measurements. It is useful for comparing brain response amplitudes overall, especially those with different temporal dynamics or following different stimulus presentation rates. Note that it is does not preclude, but rather is complementary, to the description of harmonic responses individually.

Stimuli that are presented periodically generate periodic responses of the brain that are often complex (i.e., nonsinusoidal). To capture and describe these complex brain responses overall, (baseline-corrected) frequency-tagged harmonic response amplitude can be combined through simple summation.

The authors thank those who gave feedback on earlier drafts of this article: Carrie Georges, Aliette Lochy, Milena Dzhelyova, Fang Jiang, Michael Webster, Anthony Norcia, three anonymous reviewers, and journal editors. Synthetic signals were generated in MATLAB, R2019b (The Mathworks).

Reprint requests should be sent to Talia L. Retter, University of Luxembourg, Maison des Science Humaines, 11 Porte des Sciences, 4366 Esch-sur-Alzette, Luxembourg, or via e-mail: [email protected].

Talia L. Retter: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Project administration; Software; Visualization; Writing—original draft; Writing—review & editing. Bruno Rossion: Conceptualization; Funding acquisition; Methodology; Writing—review & editing. Christine Schiltz: Conceptualization; Funding acquisition; Resources; Supervision; Writing—review & editing.

This work was supported by the Face Perception INTER project (INTER/FNRS/15/11015111 to C. S.), funded by the Luxembourgish Fund for Scientific Research (https://dx.doi.org/10.13039/501100001866) and the Belgian National Foundation for Scientific Research (https://dx.doi.org/10.13039/501100002661), grant number PDR T.0207.16 FNRS to B. R.

A retrospective analysis of the citations in every article published in this journal from 2010 to 2020 has revealed a persistent pattern of gender imbalance: Although the proportions of authorship teams (categorized by estimated gender identification of first author/last author) publishing in the Journal of Cognitive Neuroscience (JoCN) during this period were M(an)/M = .408, W(oman)/M = .335, M/W = .108, and W/W = .149, the comparable proportions for the articles that these authorship teams cited were M/M = .579, W/M = .243, M/W = .102, and W/W = .076 (Fulvio et al., JoCN, 33:1, pp. 3–7). Consequently, JoCN encourages all authors to consider gender balance explicitly when selecting which articles to cite and gives them the opportunity to report their article's gender citation balance.