Integrated trimodal SSEP experimental setup for visual, auditory and tactile stimulation

Objective. Steady-state evoked potentials (SSEPs), the brain responses to repetitive stimulation, are commonly used in both clinical practice and scientific research. Particular brain mechanisms underlying SSEPs in different modalities (i.e. visual, auditory and tactile) are very complex and still not completely understood. Each response has distinct resonant frequencies and exhibits a particular brain topography. Moreover, the topography can be frequency-dependent, as in case of auditory potentials. However, to study each modality separately and also to investigate multisensory interactions through multimodal experiments, a proper experimental setup appears to be of critical importance. The aim of this study was to design and evaluate a novel SSEP experimental setup providing a repetitive stimulation in three different modalities (visual, tactile and auditory) with a precise control of stimuli parameters. Results from a pilot study with a stimulation in a particular modality and in two modalities simultaneously prove the feasibility of the device to study SSEP phenomenon. Approach. We developed a setup of three separate stimulators that allows for a precise generation of repetitive stimuli. Besides sequential stimulation in a particular modality, parallel stimulation in up to three different modalities can be delivered. Stimulus in each modality is characterized by a stimulation frequency and a waveform (sine or square wave). We also present a novel methodology for the analysis of SSEPs. Main results. Apart from constructing the experimental setup, we conducted a pilot study with both sequential and simultaneous stimulation paradigms. EEG signals recorded during this study were analyzed with advanced methodology based on spatial filtering and adaptive approximation, followed by statistical evaluation. Significance. We developed a novel experimental setup for performing SSEP experiments. In this sense our study continues the ongoing research in this field. On the other hand, the described setup along with the presented methodology is a considerable improvement and an extension of methods constituting the state-of-the-art in the related field. Device flexibility both with developed analysis methodology can lead to further development of diagnostic methods and provide deeper insight into information processing in the human brain.

Objective. Steady-state evoked potentials (SSEPs), the brain responses to repetitive stimulation, are commonly used in both clinical practice and scientific research. Particular brain mechanisms underlying SSEPs in different modalities (i.e. visual, auditory and tactile) are very complex and still not completely understood. Each response has distinct resonant frequencies and exhibits a particular brain topography. Moreover, the topography can be frequency-dependent, as in case of auditory potentials. However, to study each modality separately and also to investigate multisensory interactions through multimodal experiments, a proper experimental setup appears to be of critical importance. The aim of this study was to design and evaluate a novel SSEP experimental setup providing a repetitive stimulation in three different modalities (visual, tactile and auditory) with a precise control of stimuli parameters. Results from a pilot study with a stimulation in a particular modality and in two modalities simultaneously prove the feasibility of the device to study SSEP phenomenon. Approach. We developed a setup of three separate stimulators that allows for a precise generation of repetitive stimuli. Besides sequential stimulation in a particular modality, parallel stimulation in up to three different modalities can be delivered. Stimulus in each modality is characterized by a stimulation frequency and a waveform (sine or square wave). We also present a novel methodology for the analysis of SSEPs. Main results. Apart from constructing the experimental setup, we conducted a pilot study with both sequential and simultaneous stimulation paradigms. EEG signals recorded during this study were analyzed with advanced methodology based on spatial filtering and adaptive approximation, followed by statistical evaluation. Significance. We developed a novel experimental setup for performing SSEP experiments. In this sense our study continues the ongoing research in this field. On the other hand, the described setup along with the presented methodology is a considerable improvement and an extension of methods constituting the stateof-the-art in the related field. Device flexibility both with developed analysis methodology can lead to further development of diagnostic methods and provide deeper insight into information processing in the human brain.
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Evoked potential (EP) is a neural response to a stimulus observed in electroencephalographical (EEG) and magnetoencephalographical (MEG) recordings. It is widely used in clinical practice and research applications to assess brain activity related to the occurrence of a stimulus. Steady-state evoked potential (SSEP) is a brain response evoked by repetitive sensory stimulation with a constant frequency. As an oscillatory waveform, SSEP response has the same fundamental frequency as the stimulus [14, 16, 47-49, 55, 60, 62]. Three types of SSEP can be distinguished as stimulation modalities: (1) steady-state visually evoked potentials (SSVEPs), which occur after looking at a visual stimulus flickering with constant frequency [48], (2) auditory steady-state responses (ASSRs), induced by repetitive auditory stimuli [16], and (3) somatosensory steady-state evoked potentials (SSSEPs), observed as a response to a repetitive tactile stimulation [14].
The strongest amplitude response of SSVEPs has been observed in stimulus frequency ranges of 5-10 Hz, 13-25 Hz and 40-60 Hz by [49] and in the range of 5-25 Hz with a peak at around 15 Hz by [64] and [42]. Herrmann et al identified a number of resonant frequencies in SSVEP signals, including 10, 20, 40 and 80 Hz [19]. Auditory periodic stimulation evokes maximal responses for stimulus rate around 40 Hz [16], while 26 Hz was suggested an optimal frequency to evoke SSSEP [55]. EEG/MEG source modeling, fMRI and PET studies identified sources of SSVEP and SSSEP in cortical areas [2,10,12,20,36,42,58,59]. Generation of auditory steadystate responses might additionally involve subcortical structures [18]. Furthermore, a common feature among steady-state responses in all three modalities is a presence of harmonic components in their spectra. These harmonic frequencies can be triggered even by sine stimulation, what suggests their nonlinear origin. On the other hand the properties of harmonic spectra in the three modalities exhibit subtle differences. In SSVEP and ASSR magnitude of harmonic peaks increases at resonant frequencies and monotonically decreases away from them [7,19,24,42]. On the contrary, in somatosensory system the second harmonic component may be sometimes stronger than the fundamental response [6,25,61]. The strong frequency-doubling phenomenon could possibly arise due to rectification of 'on' and 'off' responses in rapidly adapting sensory receptors [43].
Due to their unique properties, SSEPs are widely used in basic research and clinical applications. ASSRs can be used as a diagnostic method in audiology [47], while SSVEPs have found applications in neurology [63]. Compared to EPs, SSEPs are characterized by a higher signal-to-noise ratio (SNR) and are less affected by low frequency artifacts, like eye blinks. SSEPs are commonly used in neuropsychological research to study cognitive processes and monitor subject's attention [63]. In addition to EPs and event-related desynchronization/synchronization, they are one of the basic paradigms in brain-computer interface (BCI) applications [1,8,15,32,33,46,65]. SSEPs can also help with understanding the neural mechanisms responsible for their occurrence, which can provide insight into the functioning of the brain.
For SSVEPs, several aspects of stimulation are crucial in different types of applications. In clinical practice, where photoparoxysmal response is evoked, a flash lamp (i.e. Grass, PS33-PLUS, www.grasstechnologies.com) is commonly used [21]. This technology allows for the generation of stimuli in a broad frequency range (∼1-60 Hz) with a constant intensity, but the shape of the stimulus cannot be adjusted. For BCIs, the complexity of the system affects stimulator selection. Wu et al compared three different types of displays: light-emitting diode (LED), liquid crystal display (LCD) and cathode ray tube (CRT). Their results showed that CRT and LCD induce lower fundamental frequency amplitude in SSVEPs, which is consistent with the differences among the stimulators themselves [66].
The core of investigation on ASSR has been audiological research. Sound stimuli are provided by commercial closed systems (i.e. Neuroaudio, www.neurosoft.com; R330, www. metrohealth.com.sg; Sentiero, www.pathme.de). Repetitive auditory stimuli are also used in EEG-based BCI systems. In this case, stimuli generated by a computer and delivered via headphones or speakers are sufficient for getting a response. Systems used in clinical practice are adjusted only to standard diagnostic procedures and do not allow for flexible control of stimulation parameters. However, ASSR can be evoked with many different types of stimuli, such as: beats, chirps, clicks, pulses, sine waves, sinusoidally modulated noise, narrowband noise, frequency modulated tones, amplitude modulated tones or tones with mixed modulation (e.g. 100% amplitude and 25% frequency modulation) with or without phase difference between modulation types [44]. The magnitude of the ASSR strongly depends on the stimulus but commercial systems do not allow for a precise control of generated stimuli types and their parameters. Moreover, such systems usually record EEG signals during stimulation from limited number of electrodes (mainly Fz and Cz).
Although research into SSEPs has mainly focused on the practical application of the results in such areas as medical diagnostics and BCIs, the three modalities of SSEPs have most commonly been examined separately. By presenting sensory stimulations simultaneously in different modalities one can investigate whether the responses are enhanced, reduced or unchanged with respect to unisensory stimulation. In this way interactions between sensory processing across the senses can be studied. Although multisensory interactions were investigated in a number of studies (e.g. [31,51]), they were rarely examined using periodic stimulation. Combined SSEP paradigm in two modalities has been used in only few studies [6,17,38,50] with limited set of stimuli parameters.
Our integrated multimodal setup is well suited to investigate processes of multisensory integration. It allows for simultaneous delivery of repetitive stimulation in three different sensory modalities: visual, tactile and auditory with a precise control of stimuli parameters. Additionally, our setup may be used to investigate whether unisensory evoked steady state responses in various modalities share common characteristics in individual subjects. Such possibility was suggested by [16] more than three decades ago, but no studies to examine this idea have been set out, possibly due to a lack of proper experimental setup. The scheme of the ASSR stimulator is presented in figure 1(A). It includes two main paths: acoustic (marked on figure 1(A) with a red line) and synchronization (marked on figure 1(A) with a blue line). The acoustic signal path consists of a D/A converter integrated with a headphone amplifier (Asus Xonar Essence One), a 50 dB passive attenuator, and audiometric headphones (Sennheiser HDA280). The attenuator serves two functions. It lowers the signal level and it improves the impedance matching between the headphone amplifier and the headphones. The acoustic path uses one of the channels of the Asus Xonar Essence One converter. The second path provides synchronization between stimulus presentation and registration of the EEG signal. It is connected to the second channel of the Asus Xonar Essence One converter. In order to provide electrical isolation between the sound card and EEG amplifier, the synchronization path is supplied with a linear optocoupler. The sound stimulus is generated digitally using a PC-compatible computer with a 44.1 kHz sampling rate and is led to the D/A converter through a USB 2.0 interface. The whole system is powered by an isolated power supply unit.

Methods
Before the experiments, the ASSR stimulator was calibrated as follows. Headphone measurements were performed in an acoustically isolated room with the use of an artificial ear (B&K, type 4153 with a DB 0843 adaptor for circumaural headphones) coupled with a half-inch microphone (B&K, type 4134), a microphone preamplifier (B&K, type 2619) and a spectrum analyzer (B&K, type 2144). The spectrum analyzer was equipped with software for FFT and CPB analysis (7651 and 7667 software options). The complete stimulator setup was calibrated using a pistonphone (B&K, type 4220). The spectrum analyzer and measuring amplifier (B&K, type 2606) were used for measuring signals from the headphone amplifier output and attenuator output. The non-linearity of the audio path (which mainly refers to the headphones) manifested itself in the presence of the second and third harmonics with overall total harmonic distortion (THD) below 0.16% in the frequency range of 200-1250 Hz.
The comparison of the intended and generated waveform is presented on the figure 2(A). The signal in the form of a full-scale sine wave at a frequency of 1 kHz directed to the D/A converter generated an output signal of 22.1 mV rms (with the headphone amplifier's volume level knob in the maximum position). The sound pressure levels in the headphones were 92.5 and 93.3 dB for the left and right channel respectively. The frequency response curve irregularities in the range of 200-1250 Hz did not exceed 1.7 dB (see figure 2(B)). An example of an audio signal used to generate auditory stimuli and its spectral content is presented in figures 2(C) and (D). The stimulus was a 480 Hz tone modulated according to the following formula: where A 0 is the amplitude of the carrier wave, t is the time, f c = 480 Hz is the carried frequency and x(t) is the modulating signal. The modulating signal could have been a sine or square wave obtained as follows: • sine modulation wave: • square modulation wave: where f s is the frequency of the modulating signal and k is the highest harmonic number used for the generation of a square wave. As can be seen, we obtained a square wave as a product of the finite Fourier time series. This is due to the fact that the spectrum of an amplitude modulated signal contains harmonic components with frequencies f c ± i · f m , where f m is the spectral content of the modulating wave and i = 1, 3, 5, ...∞. To avoid aliasing, which could occur when In order to minimize the impact of transient states on the quality of the delivered stimulus, the audio signal was made smooth during the first 100 ms by means of a Tukey window.

Visual stimulator.
The scheme of the stimulator dedicated to evoking SSVEP is presented in figure 1(B). It consists of a custom-made stroboscope lamp, an arbitrary waveform generator (AWG) and a linear optocoupler. The display (10 cm by 10 cm) of the stroboscope lamp is backlit by a single LED, emitting warm white light. To obtain a uniform illumination of the display, two diffusers are positioned in front of the LED. The specially designed LED control electronics provide the stroboscope lamp with a light output linearly proportional to the input voltage (THD less than 5%) in the frequency range of 1-240 Hz. The maximal luminance of the stimulus emitted by the display and measured at its surface is 250 lx. All stimuli are normalized to have the same total luminescence per their period time. The digital stimulus signals are generated by a PC-compatible computer and led to the AWG through a USB 2.0 interface. This circuit is marked on figure 1(B) with a red line. Such a construction of a SSVEP stimulator gives the unique opportunity to generate stable light stimuli of various time courses, frequencies, intensities and modulation depths. The synchronization of the stimuli and the EEG recording is provided by the circuit marked in figure 1(B) with a blue line. The circuit includes a photodiode (not shown in the figure) placed inside the stroboscope lamp. The photodiode measures the light emitted by the LED and sends the corresponding electrical signal by the linear optocoupler directly to the EEG amplifier. In this way, the EEG signal is recorded with respect to the exact time course of the light stimuli. Examplary visual stimuli and their spectral content are presented in figure 3. Stimuli were generated according to the following formula: for the sine wave and for the square wave. A 0 is the light intensity.  human hand. The tactile stimulator is designed to stimulate the finger pad of the right or left hand's forefinger. The stimulation is delivered by a flat-tipped plastic cylinder, 6 mm in diameter, which comes into contact with the distal finger pad. Preliminary static skin tension on the finger pad is achieved through the use of two measures. Firstly, the plastic cylinder is surrounded by a flat ring, with an internal diameter of 7 mm. The gap of 1 mm between the edge of the ring and the side of the cylinder causes a desirable skin deflection. Secondly, the index finger is pressed from above by static force, provided by a suitable system of weights. The system allows the static force to be adjusted from 0 to 2.3 N. The digital stimuli signals that control the mechanical vibrations are generated by the PC and delivered to the stimulator by means of the AWG (see the blue line in figure 1(C)). The movements of the plastic cylinder are monitored by the accelerometer (Type 8636C5; Kistler Instruments) and then visualized on an oscilloscope. In order to synchronize the displacement of the plastic cylinder with the registration of the EEG signal, an acceleration signal is transmitted through the linear optocoupler to the EEG amplifier. The amplitude of vibration as well as initial displacement amplitude of the plastic cylinder is analytically determined by double integration of the signal from the accelerometer. The maximal static displacement of the flat-tipped cylinder is ±1500 μm. The maximal vibration amplitude of the cylinder is 1000 μm peak-to-peak. The non-linearity of stimulator manifested with overall THD below 5% in the range of 1-240 Hz. Due to interindividual variability of the skin and shape of the forefinger, the parameters of each type of stimulus are fitted to each participant individually, using custom-made calibration software. Examples of tactile stimuli and their spectral content are presented in figure 4. Stimuli were generated according to the formula: (4) and (5), provided A 0 is the amplitude of the tactile stimulus. Once the periodic stimulation begins, the plastic cylinder needs to accelerate slowly in order to minimize instability of the delivered tactile stimulus. To achieve this, the digital signal controlling the plastic cylinder starts to move 150 ms prior to the stimulation onset. Its time course over these 150 ms is a cosine lobe.

Electrophysiological experiment
Apart from standard testing procedures of the setup itself, we also conducted an electrophysiological experiment. In the experiment, three types of stimulation modality (visual, auditory and tactile), and two types of waveforms (sinusoidal and square) with three different frequencies (15,24,40 Hz) were used, which yielded 18 stimulation types. Five right handed adults of both sexes (4 women and 1 man, mean age = 26, standard deviation = 6) participated in the experiment. All were screened for photogenic epilepsy. They declared a lack of neurological and psychiatric disorders and that they were not on any medication that affected brain electrical activity. Written consent to take part in the experiment and the information about the procedure was signed by all of them. The EEG signal was recorded by means of Ag/Ag-Cl electrodes placed on the human scalp in 47 positions selected from the international 10-10 system (see figure 5). We selected the following electrodes: AFz (as a reference), F3, Fz, F4, FC5, FC3, FC1, FCz, FC2, FC4, FC6, C5, C3, C1, Cz, C2, C4, C6, CP5, CP3, CP1, CPz, CP2, CP4, CP6, P5, P3, P1, Pz, P2, P4, P6, PO3, POz, PO4, O1 and O2, which are localized above brain areas responsible for generating particular SSEP responses. The impedance of the electrodes was below 5 kΩ. The signal was acquired by a two coupled Porti7 (TMSI) amplifiers with a sampling frequency of 2048 Hz. The testing protocol consisted of two separate stimulation paradigms: sequential and simultaneous. In the former, sequential, three different types of stimulus were used sequentially: visual, auditory and tactile. In the latter paradigm, visual and tactile or visual and auditory stimuli were exposed simultaneously. In both experimental paradigms, the participant sat on a comfortable chair in a darkened and soundproofed room. Before the start of the experiment, a stroboscopic lamp was positioned 50 cm in front of the participant, the headphones of the auditory stimulator were put on and the participant's index finger was placed on the tactile stimulator. At the beginning of the procedure, the participant was instructed to concentrate on the cued stimulus and was asked to avoid making any movements or blinking. In this experiment, we used the following physical properties of stimuli. The intensity of the light emitted by the visual stimulator was 250 lx. Auditory stimuli were presented diotically at a sound pressure level of 50 ± 0.5 dB SPL. Tactile stimuli had an amplitude of 500 μm peak-to-peak.

Sequential experimental paradigm.
In the sequential experimental paradigm a single trial was made up of three phases (see figure 6(A)): sound information, a 1 s-long break, the stimulation period and an inter-trial interval. Delivered to the left headphone, a 500 ms long voice cue indicated the modality of the stimulus (visual, auditory or tactile) to which the participant should pay attention. Next, after 1 s-long break, the visual, auditory or tactile stimulus was presented to the participant. Despite modality, the delivered stimuli were differed in frequency (15 Hz, 24 Hz or 40 Hz), modulation wave (sine or square) and they were presented in a random order. Stimuli were generated according to following formulas: (2)-(5). The frequency of carrier wave of ASSR stimuli was 480 Hz.
To avoid an adaptation process, the duration of the intertrial interval was randomly selected from the range of 2.5-3.5 s. The stimulus with the same combination of parameters was repeated nine times during one session. The experiment was divided into five about 18 min -long sessions. Each session was followed by a short break of a duration determined by the participant. An example sequence of stimuli in the sequential paradigm is shown in figure 7.

Simultaneous experimental paradigm.
Two participants took part in the second experiment with a simultaneous stimulation paradigm. They were asked to focus on two simultaneously provided stimuli. Due to the fact that ASSR and SSSEP responses might be observed on the same electrodes, depending on body part being stimulated, an auditoryvisual and visual-tactile stimuli combination was used in this paradigm. The trial scheme was similar to the sequential paradigm, except that there was no voice cue at the beginning of each trial (see figure 6(B)).
The first participant was stimulated with a sine modulated wave in visual and auditory modalities, while the second one with a square modulated wave in visual and tactile modalities. To limit the duration of the experiment, for each participant three pairs of random frequencies from the testing set were chosen. The following pairs were used for the first partici-  Preprocessing the data consisted of the following steps. EEG signals were re-referenced to the common average montage and downsampled to 512 Hz. We applied a third order Butterworth filter in the range of 1-95 Hz. Signals were additionally notch-filtered by means of a fifth order Butterworth filter to remove the 50 Hz mains disturbance.
We extracted 6 s-long segments, containing a 1 s-long period preceding the stimulus onset, a 4 s-long stimulation epoch and a 1 s-long period following the end of the stimulation (see figure 6). The extracted segments were denoted as X c ∈ R k×N×M , where c ∈ {'auditory 15 Hz square modulation', 'auditory 15 Hz sinusoidal modulation', …, 'tactile 40 Hz sinusoidal modulation'} represents the 18 different types of simulations, R k×N×M is a k × N × M real number matrix, k is the number of channels, N is the number of samples, M is the  number of experimental repetitions, and c is the given stimulation type. Segments with amplitudes exceeding ±125 μV were marked as artifacts and removed from further analysis.

Detection of SSEP response.
In this study, we applied synchronized averaging and spatial filtering to detect SSEP signals. Firstly, we averaged the responses X c over exper imental trials. In the next step, we used spatial filtering based on the formula proposed in [33], originally intended to improve phase detection in SSVEP-BCI. We modified the formula published in [33] and applied it to the analysis of each type of SSEP response separately. Mathematical details can be found in appendix A. The spatial filtering technique allows for decomposition of EEG signal into components related to different brain electrcial activites. Each component corresponds to its spatial distribution in the EEG sensor space, which is called a spatial pattern. In this work we denoted component related to SSEP response and its corresponding spatial paterns as s(t) and w respectively. Spatial filtering, like other blind source separation (BSS) techniques, does not guarantee that the selected signal s(t) will contain any physiologically meaningful electrical activity of the brain. The spatial pattern w corresponding to the component s(t) can be visualized as a plot showing its spatial distribution on the scalp (called a topography). To validate obtained results, we analyzed both the topography and dynamics of the chosen component s(t).

Estimation of SSEP dynamics. Among many different methods of time-frequency analysis, matching pursuit (MP), introduced by Mallat [29]
, supplies the highest possible timefrequency resolution [4]. In addition, the method allows for accurate parametrization of structures contained in the analyzed signal. By parametrization, we mean the possibility to find parameters of structures contained in the analyzed signal, i.e. its position in time and in frequency, its amplitude and the width of its time span, defined typically as full width at half maximum. The MP algorithm has been repeatedly used in brain signal analysis [3,70,71]. The EPs analyzed in our experiment are indeed symmetrical, but they are characterized by steep slopes, so we expected similar difficulties with describing the signals using Gabor functions. We therefore decided to use the MP algorithm implementation proposed in [57]. Details can be found in appendix B.

Spatial-time-frequency representation.
For the purpose of the present study, we developed a method incorporating the best properties of spatial filtering and MP, reducing some of the drawbacks that each of the methods has when employed separately. BSS methodology allows for the separation of a multi-channel signal into components. The component s(t) with the highest SNR for SSEP response and desired topography w was chosen for further analysis. The accurate parametrization of structures presents in the signal s(t) was achieved by means of the MP algorithm. Such procedure results in a very precise spatial-time-frequency filter.
It is expected that the combination of both algorithms should result in a method that gives a deeper insight into the properties of the analyzed steady-state evoked responses.

Statistical analysis.
In this study we analyzed the amplitude of the SSEP response in the time-frequency domain. We tested the null hypothesis that the periodic stimulation does not elicit SSEP response against the hypothesis that it does. In terms of time-frequency representation it means that according to the null hypothesis there is no change in the EEG amplitude in the frequency of stimulation and its harmonics, and the alternative hypothesis states that there is increase of the power in the frequency of stimulation and its harmonics. We performed validation of the null hypothesis by means of bootstrap procedure with extreme value statistics. Details can be found in appendix C.

Results
We conducted two physiological experiments to observe the performance of the experimental setup. In the first one, a random sequence of tactile, visual and auditory stimuli with different parameters was applied to each participant. Figure 8 shows an example final outcome of an analysis of the experimental data. Starting from the top of the picture we can see the chosen spatial filtering component s(t) and its time-frequency representation obtained by means of MP procedure. The SSEP response lasts throughout the entire stimulus interval and occurs as quasi-periodic signal. The most characteristic feature of the SSEP is the increase of signal amplitude in the frequency of stimulation and its harmonics, marked respectively on the time-frequency map by structure no. 1 and structure no. 2-4. Apart from the widely used time-frequency analysis, we show that SSEP response consists of well identified and parametrized structures. The parameters of atoms matched to those structures are presented in the charts at the bottom of figure 8. For example, we can see that single atoms fitted to the spatial filtering components s(t) explain the SSEP responses in the frequency of stimulation. Their amplitudes are 19.44 µV, 3.27 µV and 2.80 µV respectively for SSVEP, ASSR and SSSEP. After projecting the atom on the channel where SSEP signals are typically measured in a given modality, we obtain an amplitude of 3.80 µV on channel O2 for SSVEP response, 0.27 µV on channel FCz for ASSR response, and 0.17 µV on channel CP3 for SSSEP response. The amplitudes of SSEP in the harmonics of stimulation frequencies are presented in successive rows of the charts. As can be seen, we detected up to the fourth, second and third harmonic respectively for SSVEP, ASSR and SSSEP responses. These results are consistent with earlier studies demonstrating that the amplitude of the SSEP depends on the stimulated modality and is a few µV in the case of SSVEP [42], while the amplitude of ASSR and SSSEP is several hundreds of nanovolts [62].
In addition to precisely parametrized amplitude and frequency of the response, we also obtained the length of SSEP potentials, which is presented in the third column (titled width) of the charts. We show that during a 4 s-long stimulation, the SSEP response ranges from 3.70 -3.88 s. This parameter is not possible to achieve using classical Fourier analysis commonly applied in the processing of SSEP responses. The obtained time-frequency patterns are in line with state-of-theart knowledge on SSEP morphology [13,22,23,42,47,55,62,63]. SSEP lasts even after the offset of the stimulation. It has been previously reported that with periodic visual and auditory stimulation, SSEP can persist for several cycles after the end of the stimulation [56]. Presence of even harmonics in the time-frequency maps, when using square modulated stimulation, cannot be attributed to the harmonics of the stimulus, as a perfect square wave does not contain even harmonics. Accordingly, their origin must be related to the nonlinearity of the stimulated system. Nonlinear transformations of the input signal may induce multiple responses for each cycle of the stimulus or simply produce evoked response waveform, which cannot be represented by a sine wave. Both phenomena will be manifested by harmonic frequencies in time-frequency representation of the output signal.
Our experimental setup makes it possible to differentiate the shape of modulating waveform. In figures 9 and 10, the time-frequency maps of a selected spatial filtering component s(t) averaged across all participants both for square and sine modulating waveforms are shown. A statistically significant response in fundamental frequency was evoked in all cases. In some cases, an increase in harmonic frequencies was also significant. However, different shapes of the stimulus elicited distinct patterns of activity. Additionally, visualizations of the spatial patterns w corresponding to the component s(t) are presented as topographic maps in the left corners in figures 9 and 10. The topographies were averaged across all participants. As shown in figures 9 and 10, the strongest activity of SSVEP response was observed in the primary visual cortex (V1) and in the motion sensitive (MT/V5) areas where it was possibly generated [10,63]. In the case of ASSR, its dominant magnitude is observed over central derivations, which is in line with other studies [30,50,52,68]. For the SSSEP responses, the strongest magnitude is observed in the left parietal lobe, where the sensory areas of the right hand are located [9,36,55,62]. One can notice that the spatial patterns have a bipolar character, with their centers localized in the area of the C3 electrode. This topography resembles the distribution of the potential derived from a horizontal dipole arranged between two poles [40]. A possible candidate for the dipole position in the brain is central sulcus [37], which is localized in the vicinity of the C3 electrode.
During the second experiment, participants were exposed to multimodal visual-auditory or visual-tactile stimuli. As shown in figure 11 during visual-auditory stimulation, magnitude

Discussion
In this study, we present a novel experimental setup dedicated to measure SSEPs in visual, auditory and somatosensory modalities. Our setup provides a precise shape of generated stimuli and manifests the THD below 0.16%, 5% and 5% respectively for auditory, visual and tactile modalities. Unlike other available devices, the flexibility of stimulus generation with particular parameters allows responses in different modalities to be directly compared. Apart from conducting experiments in each modality separately, the setup makes it possible to design multimodal studies.
To the best of our knowledge, only few studies focused on bimodal stimulation. Saupe et al investigated parallel appearance of SSVEP and ASSR [50], while similar experiments On the timefrequency maps, horizontal axis-time, vertical axis-frequency, the intensity of grey represents the amplitude of the SSEP signal, as shown on bars next to each map. The stimulation onset is at 1st second and the stimulation offset at 5th second. Areas that were statistically not significant were whiten out. On the topography map, color represents the amplitude of the response on the human scalp: red-high amplitude, green-low amplitude, blue-high amplitude of SSEP response with negative polarization.   on SSSEP and ASSR has been reported in [6,61]. Our setup additionally enables to investigate simultaneous SSVEP and SSSEP recordings as well.
For visual modality, our setup allows for the generation of various types of stimuli, including sine, square and triangle waveforms and their superposition. Additional parameters, like phase of stimulation or duty cycle of waveform are adjustable. Stimuli can also be provided in jittered sequences. For auditory tones, beats, chirps, amplitude-, frequency-and mixed-modulated sine waves can be generated. The tactile stimulator makes it possible to generate sinusoidally-and square-modulated waveforms.
In this study, we present also an experiment where SSEP stimuli in three modalities were provided in homogenous conditions, which means that the stimuli had the same properties in terms of stimulation frequency and stimuli shape. The results obtained both for each modality and for the simultaneous stimulation of two modalities are consistent with previous studies. The spectral contents of the obtained SSEP responses and their topographical distributions on the scalp are typical for each modality and reveal the neurophysiological processing of sensory inputs.
To improve detection of the SSEP response, we applied a variant of a common spatial filtering method to the averaged data. The selected component was decomposed using the MP algorithm to reveal its time-frequency properties. Apart from the detection and parametrization of the SSEP signal, the main goal of the applied methodology was to analyze the  Example of MP signal decomposition with atoms from standard and enlarged Gabor dictionary. Panel (A)-the decompositions of a simulated function: on the left-representation with atoms from the standard Gabor dictionary, on the right-representation with atoms from the enlarged dictionary. Panel (B)-actual amplitude and time-frequency representations. Panel (C)-resulting functions. In the Gabor dictionary, three functions are necessary to account for 95% of the energy of the simulated signal. Only one function is necessary to account for 99% of the energy in the second case. SSEP response by means of uniform tools. It has been shown that differences between response to stimulus and additional EEG activity can be made more evident if the multichannel data are subjected to a suitable transformation-spatial filtering. Spatial filtering has been widely applied in BCIs based on SSVEPs [8,15,26,69], but to the best of our knowledge there are only few reported applications of this method to the analysis of vibrotactile [35] or auditory [67] responses. MP is characterized by a high time-frequency resolution and has found various applications in EEG signal parametrization [3,4,11,27,28,53,54,70,71]. However, it has not yet been combined with spatial filtering and applied to analysis of SSEP response.
Standard methods of SSEP analysis in the time-frequency domain are based on Fourier methodology. This methodology implies the need to set some arbitrary parameters, i.e. window length, window shape and window shift. Such parameters do not change and are not matched to the structures contained in the analyzed signal. However, MP estimates the time-frequency representation that adopts decomposition to the local signal structures. This results in a representation which is basically free from any arbitrary settings, apart from the dictionary construction.
In conclusion, we propose a multimodal experimental setup with advanced methodology for SSEP data analysis. The results indicate the validity of using such a precise (in terms of the characteristics of the delivered stimulus and control of stimulation parameters) setup for stimuli application. Both the experimental setup and the methodology of data analysis seem to be effective tools for conducting SSEP research.
The experimental setup does not rely on any specific equipment. It could be easily reconstructed using commonly available parts and devices. Components used are not of critical importance either and can be replaced by others which are characterized by similar parameters. Calibration of the device also requires only a standard electro-acoustic or electromechanical measurement equipment. Readers interested in constructing such setup can contact the corresponding author to receive its full technical documentation.

Acknowledgments
This study was supported by the National Science Centre, Poland, grant number: 2012/07/D/NZ4/04226.

Appendix A. Spatial filtering
The classical approach to enhancing the SNR of the response is to average many stimulus-aligned epochs of the EEG signal containing the response. This technique is based on the assumption that the EEG signal is composed of two additive components: the invariable response (i.e. structure due to the stimulus that always occurs at the same latency and with the same morphology) and the spontaneous EEG signal, which can be treated as an uncorrelated, stochastic zero-mean process. The first step of our analysis is the synchronized averaging of the signal X c ∈ R k×N×M over the experimental trials.
After performing this step, for each experimental condition we obtain the data matrix Y c ∈ R k×N .
Differences between SSEP and EEG background activity can be made more evident if the multichannel data are substituted for a suitable transformation. The synchronous electrical activity of neurons in a patch of the cerebral cortex can be approximated by a current dipole [39]. Each of such dipoles can be treated as a source, providing a contribution (comp onent) to the measured electrical potential distribution on the surface of the head. Assuming that sources are simultaneously active and statistically independent, the potential measured at the scalp can be modeled as a linear combination of current sources: where rows of matrix Y correspond to the recorded EEG signal, rows s of matrix S represent source activity, matrix W is a mixing matrix and its element w ij describes the contrib ution of the jth source to the signal recorded by the ith electrode. To find independent sources, the following transformation is performed: In practice, we know neither source signal S nor mixing matrix W. Estimating matrix Ŵ −1 and source signal Ŝ from measured signal Y is known as the blind source separation (BSS) problem. Several methods of solving this problem, such as principal component analysis and independent comp onent analysis, are used in electrophysiology [41]. Here we use spatial filtering, known as common spatial patterns [33], which yields a linear combination of original channels so that the variance of the EEG activity is maximized while simultaneously the magnitude of the background activity is minimized.
We exploit the methodology of estimation matrix Ŵ −1 proposed in [33], where an EEG signal recorded during stimulation is modeled as a linear combination of an SSEP signal and noise [15,32,33]. We select from Ŝ the component s(t) with the highest SNR for SSEP response. Each column of matrix Ŵ gives a contribution of the corresponding source to the EEG signal and is called a spatial pattern.

Appendix B. Matching pursuit
The MP algorithm decomposes a signal into a linear combination of functions (commonly called atoms) selected from a predefined set of such functions-the dictionary. The dictionary should be constructed from functions with an appropriate localization in time-frequency space. It should also be redundant, like a base spanning the signal space. The decomposition is performed by means of an iterative procedure and after m iterations, the signal f (t) can be represented as follows [29]: where α i is a weighting factor, g i (t) is the atom selected in the ith iteration and R(t) is the residual signal left after the procedure. MP has the desired property of explaining the signal in terms of time-frequency localized structures. Due to the redundancy of the dictionary, such an explanation is also very flexible. Many different implementations of the MP algorithm exist, but mostly the so-called Gabor functions are used for signal decomposition. Such functions are suitable for describing waveforms present in the signal, and at the same time provide the best possible time-frequency resolution [28]. There are, however, cases when choosing to use Gabor functions is not optimal. Recently it has been shown that signals containing asymmetric components would be insufficiently described by functions from the Gabor family [57]. Taking into account the actual shape of the SSEP we invented a new set of functions described by the following formula: where μ is the position in time, f is the position in frequency, σ is the time scale (time span), and ϕ 0 is the initial phase. The EPs analyzed in our experiment are indeed symmetrical, but they are characterized by steep slopes, so we expected similar difficulties with parametrizing signals with Gabor functions. We therefore decided to use the MP algorithm implementation proposed in [57], as it provides an easy way to change the dictionary. An example of a function described by form ula (B.2) is shown in figure 13. In figure 14, the decomposition of a simulated signal into components is shown. The simulated signal ( figure 14(A)) was not constructed from a function present in the dictionary, but it was created by the multiplication of a 4 s long Tukey window [5] and a sinewave of a frequency of 5 Hz. In the figure 14(B) the time-frequency distribution of signal amplitude is presented. The distribution of amplitude was calculated based on the equation (11) from the article [57].
It is clearly visible that the result of the MP algorithm based on the enlarged dictionary is more adequate. The dictionary composed with Gabor functions is not able to precisely represent the onset of the simulated structure. The fact that only one function was needed to decompose 99% of the energy of the signal is an additional benefit of the enlarged dictionary.