The MEG recordings were performed in a magnetically shielded room (MSR) (Vacuum-Schmelze, Hanau, Germany) [19, 20]. A portable MEG system was developed using 10 OPM sensors, a NI 9205 Data Acquisition Unit, and a notebook computer. The OPM sensors (Gen-2) were provided by QuSpin (QuSpin. Inc, Louisville, CO, USA) [4, 24]. The technical details have been described in previous reports [4, 24]. The field sensitivity of the sensor was 15 fT/$\surd$Hz (Z-direction) in the 3–100 Hz band (typically 7–10 fT/$\surd$Hz). The dynamic range was $\pm$5 nT. A helmet was developed for holding the OPMs. The main frame of the helmet was 3D printed (Fig. 1).

Fig. 1.

A photo showing the helmet, coils (fiducial points), source (dipole), and optically pumped magnetometer (OPM). The coil is used to localize the fiducial points relative to the position and orientation of sensors. The source generates predefined dipolar signals for developing and verifying methods (e.g., the location of the source is known). The OPMs are used to detect the signal from the dipole to assess the method for artifact reduction.

The position of each OPM was adjustable for flexible detection of neuromagnetic signals from a desired brain region (e.g., the auditory cortex). Three coils were attached to the nasion and to the right and left pre-auricular points of each subject to locate the subject’s position head relative to the sensor position and orientation. Synchronized data acquisition from the OPMs was accomplished with a National Instruments (NI) card cDAQ-9171 (Austin, TX, USA). The OPMs (magnetometers) were linked to the NI card, which was then connected to a notebook computer (Surface 2.0, Microsoft Corporation, Redmond, WA, USA). A software package was developed to acquire data, perform dynamic artifact reduction and real-time signal processing. To identify system and environmental artifacts, two MEG datasets without phantom or subject were recorded just before MEG tests (empty room recordings).

We used a current dipole phantom, which came with our CTF MEG system (CTF, Coquitlam, BC), to assess the performance of artifact reduction for the portable OPM MEG prototype [8]. A helmet was placed on the outside surface of the current dipole phantom for testing. The current dipole phantom was a spherical container filled with a conducting saline solution which contained a current source and sink. The globe had a 130 mm inner diameter. The current dipole itself was constructed of two gold spheres about 2.0 mm in diameter, separated by 9.0 mm centre to centre. The position of the current dipole could be adjusted within the globe. The location of the dipole was recorded relative to the center of the sphere (0, 0, 0 mm). The signal from the phantom was a sine waveform at 23 Hz. Two datasets were recorded without turning on the source (current dipole) before the phantom experiments to assess magnetic artifacts.

Six healthy subjects (age: 12–55 years, mean age: 30 $\pm$ 18 years; 3 males and 3 females) were recruited for this study. A written informed consent, approved by the Institutional Review Board at Cincinnati Children’s Hospital Medical Center (CCHMC), was obtained from each participant prior to testing. Similar to previous reports [26, 27], inclusion criteria for participation were: (1) healthy, without history of neurological disorder or brain injury; (2) normal hearing, vision, and hand movement. Exclusion criteria for participation included: (1) inability to sit comfortably for about 30 minutes; and (2) unidentifiable magnetic noise. Participants were instructed to listen to a 500 Hz square wave tone cue. The stimuli consisted of 200 trials of tones; 100 tones per ear, which were presented randomly through plastic tubes and earphones. Stimulus presentation and response recording were performed using BrainX software (Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA), which was based on DirectX (Microsoft Corporation, Redmond, WA, USA) [26, 27]. Participants were instructed to sit comfortably and natural small movement was allowed. Two datasets without hearing tones were recorded before the experiments for analyzing magnetic artifacts.

According to our pilot data from empty room recordings, each OPM had intrinsic artifacts itself (e.g., from the laser, vapor, and heater). Since the frequency ranges of the artifacts varied among OPMs and dynamically changed over recordings, it was difficult to remove them in the time domain (e.g., a band-stop filter requires the high and low frequency edges). However, in Fourier space, we could easily remove the artifact by processing spectral data in each frequency bin, the smallest frequency unit in spectrogram. To precisely remove artifacts in Fourier space, we transformed MEG data from the time domain to the frequency domain using the Fourier transform. Artifacts could be analyzed and identified by comparing experimental spectrograms (e.g., from phantom recordings) with empty room spectrograms. Once artifacts were identified and removed, we then transformed the frequency domain data back into the time domain for other data analyses [15, 19]. Before Fourier transformation, MEG data were high pass filtered with a corner frequency (CF) to remove slow variations at frequencies CF for time-frequency analysis. The numbers of points (L) of the data sample were matched to the chosen CF, otherwise the MEG signal would be distorted and its amplitude diminished. A suitable matching condition was L $\ge$Sampling-Rate/CF or frequency resolution $\le$CF [15, 19]. Instead of using spectral power (magnitude) for artifact reduction, the present study used the real and imaginary portion of data for artifact reduction so as to maintain the polarity information of the MEG data [19, 20].

Common mode artifacts (e.g., powerline noise, head movement) typically affect MEG data from all OPMs. Assuming MEG data include both brain signals and artifacts, we can partition the data into separate subspaces that originate from inside the head (brain signals) and outside the head (artifacts) [8, 11]. Fig. 2 shows the spatial relationship between brain signals and artifacts. Mathematically, we can determine basis sets used to distinguish between external signals and internal signals. In this study, MEG data from the OPM were therefore divided into outside MEG data (whose sources are out of the head) and inside MEG data (whose sources are within the head). The two types of data are supposed to belong to separate subspaces (Fig. 2).

Fig. 2.

A diagram illustrating the spatial relationship between sensors, sources, and artifacts. All sensors are optically pumped magnetometers (OPMs) which can be placed and adjusted to cover a specific brain region (e.g., “1”, “2” and “3”). Sources are the signals of interest, which are from the brain in the head subspace (e.g., “S”). Artifacts are the noise and interferences from the environment and/or from the sensors (e.g., “A”), which are from the subspace out of the head.

Building on previous reports on SSS and DSSP [11, 18], the relationship between the measured signals B and artifacts can be described with the following equation.

(1)$$B(t)=Bi(t)+Bo(t)$$

where t represents time slices. Bi(t) represents inside signals from the brain. Bo(t) represents outside signals from artifacts. Brain signals can be estimated with lead field within the brain, which is the subspace of the head.

(2)$$Bi(t)=LQ(t)$$

(3)$$Q(t)=Bi(t)L^{+}$$

Where L is the lead field of sources (Q). Since the orthonormal signals of the source matrix Q(t) are assumed to be not related to source signals, these orthonormal signals are considered to be artifacts (A) [9]. The orthogonality of the brain signal (or the lead field) and artefact are mathematically calculated with cross-products (e.g., a lead field vector or a matrix of multi-lead field vectors). Consequently, MEG signals can be classified as two subspaces: brain signals and environmental artifacts. The brain signals can be obtained with the following equation.

(4)$$\widehat{Bl(t)}=B(t)\left(I-A{A}^T\right)$$

Where $\widehat{Bl(t)}$ represents estimated brain signals and B(t) represents measured signals. A${}^T$ represents the translated matrix of artifacts. $I$ represents identity matrix. From beamforming point of view, $\widehat{Bl(t)}$ represents the subspace of the sources (signal of interest); AA${}^T$ represents the subspace of the artifacts (interference). In the present study, the separation of sources and artifacts is called signal space classification (SSC).

Building on previous reports on artifact reduction in the time domain, we have derived our artifact reduction methods (Eqns. 3,4) [9, 11, 21]. In the present study, we moved one step further by removing artifacts in the frequency domain. Specifically, before artifact removal, we transferred MEG data from the time domain to the frequency domain. The removal procedure is based on a general input-output model with spectral data. Assuming a signal from a primary sensor contains both artifacts and brain signals; the relationship between the signal, brain activity and artifact can be described as the following equation.

(5)$$\widehat{Bl(f)}=B(f)\left(I-A(f)A{(f)}^T\right)$$

where B(f) represents the raw signal in frequency domain, $\widehat{Bl(f)}$ represents the brain signal in the frequency domain, and A(f) represents the matrix of artifacts in the frequency domain. $I$ represents identity matrix. Our method includes a few novel options for artifact removal. Specifically, signals in certain frequency bands could be predefined for removal by analyzing spectrograms (e.g., artifacts identified in empty room recordings). For signals with variable frequency bands, the frequency bands were dynamically identified with SSC [15, 16]. A(f)A(f)${}^T$ represents the artifacts in the subspace of external signals, which are orthogonal to the subspace of interest signals (e.g., the brain) defined by lead field. Since the method enabled the removal of artifacts in a specific frequency band and subspace, this method was named frequency specific signal space classification (FSSSC). The method has been implemented in MEG Processor, which is publically available (https://sites.google.com/site/braincloudx/).

Shielding factor was used to quantify the level of suppression of artifacts. The algorithm for computing shielding factor for artifact reduction has been described in a previous report [28]. For waveform data, the shielding factor was estimated by calculating the ratio of the Frobenius norms of the measured and processed magnetic data. For frequency data, the shielding factor was the ratio of the magnitude (power) at a specific frequency band. Mathematically, shielding factor was computed with raw data and processed data. The raw data were OPM data, whose artifacts had not been removed with FSSSC. The processed data were OPM data, whose artifacts had been removed with FSSSC. The shielding factor can be obtained with the following equation.

(6)$$SF=\frac{\parallel Bm\parallel }{\parallel Bp\parallel }$$

where SF represents shielding factor; $\parallel Bm\parallel$ represents the norm of the recorded raw data; $\parallel Bp\parallel$ represents the norm of the processed data. The norm of data is the square root of the sum of each signal squared. From a MEG data analysis point of view, the shielding factor is equivalent to the signal-to-noise ratio.

The waveforms of MEG data at sensor levels were visually inspected for distortion [19, 20]. MEG data were filtered with a high-pass filter of 3 Hz and low pass filter of 100 Hz before artifact identification. Frequency data were analysed in terms of power density. The covariance matrix was used to estimate the noise or signal correlation among sensors [19, 20]. The covariance matrix is an N$\times$N matrix, where N is the number of sensors (channels) of the MEG recordings (7$\times$7). The covariance matrix summarizes the spatial distribution in the sensor space of the noise or signal power (diagonal entries) and the spatial correlations between all MEG sensor recordings (off-diagonal entries).

Source localization was accomplished with beamforming [19, 20]. The grid spacing for the beamforming was 3 millimetres. We measured the location of the maximum activity and the distance from the maximum activity to the true sources (e.g., the current dipole in the phantom). The resulting artifacts and signals were used to estimate the source accuracy, which was quantified by directly calculating the distance between the estimated sources and the ground truth sources in a current dipole phantom.

The normality of all values was checked using the Kolmogorov-Smirnov test. The comparisons of measurements (e.g., waveform amplitude) between raw data and processed data were performed using the Student t-test. The Bonferroni correction was applied for multiple comparisons. The significance level was set at p 0.05. All statistical analyses were performed in R 3.6.1 and SPSS software package version 22.0.0.0 (IBM Corp, Armonk, NY, USA).

MEG datasets from empty room recordings without a subject showed environmental artifacts (background noise) and intrinsic artifacts associated with OPMs. Fig. 3 shows an example of waveforms from 7 OPM sensors. The amplitude of signals varied slightly among sensors and was dynamically changing with time. Since neuromagnetic signals from the brain are typically very weak (in a few hundreds of femtotesla), MEG signals with large amplitude (e.g., a few picotesla) are considered to be potential artifacts. The amplitude of one sensor (OPM2) was significantly higher than that of another sensor (OPM 3) (982 $\pm$ 171 fT vs. 704 $\pm$ 103 fT; p 0.02) (see Fig. 3). After applying FSSSC, the processed data showed waveforms with low amplitude, which seemed to be stable and clean. Statistically, the difference between OPM2 and OPM3 was not significant after applying FSSSC (615 $\pm$ 42 fT vs. 590 $\pm$ 46 fT; p $>$ 0.05). The shielding factor analyses revealed that frequency specific artifact reduction could significantly decrease the artifact in OPM MEG data (837 $\pm$ 158 fT vs. 638 $\pm$ 69 fT; p 0.01).

Fig. 3.

A block diagram of frequency specific signal space classification (FSSSC) for artifact reduction. The experimental time-series data include both brain signal and artifacts while the empty room time-series data mainly include artifacts. After performing time-frequency analysis and signal space classification, artifacts are removed and time-series data are computed for waveform analysis and source localization.

Time frequency analyses revealed artifacts in multi-frequency bands. The magnitude (spectral power) and frequency bands of some artifacts varied significantly among OPMs. The artifacts specifically associated with individual OPMs were considered to be intrinsic artifacts from the electronic constructing of the OPM (e.g., vapor, laser, heater, and wire connections). Fig. 4 shows one example of an intrinsic artifact around 8 Hz. After applying FSSSC, the processed data showed that the intrinsic artifact could be completely removed (Fig. 4). In comparison to spectrograms of raw MEG data, the spectrograms of processed MEG data showed less variation among sensors, and lower spectral power. The shielding factor analyses revealed that FSSSC could significantly decrease frequency-specific artifacts in OPMs. In the present study, the artifact reduction for signals around 7–9 Hz could attain a shielding factor of 160 (Fig. 3).

Fig. 4.

Segments of MEG waveforms from an empty room recording. The “Raw Data” show the waveforms of signals recorded from optically pumped magnetometers (OPMs) before artifact reduction. The “Processed Data” show the waveforms of signals after artifact reduction. In comparison to raw waveforms, processed waveforms have low amplitude and less variation among sensors (e.g., “OPM2” vs. “OPM3”). All waveforms are filtered with a bandpass filter of 3–100 Hz.

Fig. 5 shows an example of the artifacts in the frequency domain. We noted at least two types of artifacts. One type of artifacts appeared in a consistent frequency band in all OPMs (e.g., power line noise at 60 Hz and its harmonics). Another type of artifacts appeared in a few sensors (e.g., intrinsic artifacts). Artifacts identified on all OPMs appeared to be environmental (external) artifacts; the artifacts identified on a specific OPM seemed to be intrinsic (internal) artifacts. FSSSC developed in this study could remove both types of artifacts. Fig. 5 shows the removal of frequency specific artifacts in OPM data. In comparison to raw data, data processed with FSSSC showed a significant reduction in magnetic artifacts for all sensors. The shielding factor for reducing the artifacts reached 140 for artifact around 60 Hz.

Fig. 5.

Spectrograms showing the removal of frequency specific artifacts. The “Raw Data” shows the time-frequency representation of MEG signals before artifact reduction. The “Processed Data” shows the time-frequency representation of MEG signals after artifact reduction. A frequency specific artifact (“A”) identified in raw data is completely removed in the processed data using FSSSC. Of note, the signals in other frequency bands are almost intact. The color bar indicates the color coding of spectral power, which is identical for the two spectrograms.

Covariance analyses revealed that some artifacts were correlated among OPMs. The correlated artifacts in OPM pairs were interpreted to be environmental artifacts. We noted that FSSSC removed sensor correlated artifacts and changed the covariance patterns. Fig. 6 shows an example of covariance matrices before and after the removal of correlated artifacts.

Fig. 6.

Power density revealing external and internal artifacts in OPM MEG data. The power density computed with raw data (“Raw Data”) shows two artifacts. The internal artifact (“Internal Artifact”) is predominantly identified in OPM2 but not OPM3. This observation implies the artifact is associated with the particular sensor. The external artifact (“External artifact”) is identified in all sensors (e.g., both “OPM2” and “OPM3”), which implies the artifact is from the environment and affects all sensors. The power density computed with processed data (“Processed Data”) show that frequency specific artifacts can be completely removed by FSSSC.

The results of waveforms, spectrograms, and covariance of empty room recordings revealed that artifacts in OPM data included external (environmental) and intrinsic artifacts. FSSSC removed both external (e.g., the artifact around 60 Hz) and intrinsic (e.g., 7–9 Hz) artifacts in certain frequency bands and subspaces (Figs. 2,3,4).

Source waveforms from the current dipole in the phantom were predefined as sine waves around 23 Hz. Fig. 7 shows the source waveforms and the recorded magnetic signals during experiments. The recorded magnetic waveforms were the source signals mixed with artifacts. There were significant differences between the source signals and the recorded signals in terms of the shape and phase of the waveforms. This was anticipated because artifacts from medical devices within a hospital setting could contaminate the experimental data (Fig. 7). The correlation between the source signals and the recorded signals in raw data was 0.51 $\pm$ 0.14. After processing with FSSSC, the signals of processed data were significantly improved and were similar to the source waveform (Fig. 7). The correlation between the source signals and processed signals was 0.76 $\pm$ 0.08. In comparison to raw data without artifact removal, processed data with artifact removal showed significantly stronger correlation with the source signals (0.76 vs. 0.61, p 0.001). Since the source signals are the ground truth signals, the significant increase of correlation indicated that the FSSSC normalized and recovered the distorted signals in OPM MEG data.

Fig. 7.

Covariance data showing the changes of the spatial distribution of signals for all sensor pairs. For data from the empty room recordings, the artifact reduction method significantly decreases the amplitude. For data from the phantom and a human subject, FSSSC significantly changes the patterns of covariance matrixes. The change of covariance pattern contributes to the improvement of source localization.

Fig. 7 shows an example of covariance matrices of phantom data. Covariance matrix analyses revealed that FSSSC changed the spatial distribution. Fig. 8 shows the waveforms from a phantom source. Time frequency analyses revealed that the raw data had distorted frequency components, which were not from a single source signal. Fig. 9 shows an example of spectrograms computed for raw data and processed data. We noted that raw spectrograms included source signals and artifacts. After artifact removal (or reduction) processing, the processed spectrograms revealed a clear single frequency component which was at around 23 Hz.

Fig. 8.

MEG Waveforms from a phantom source. The raw signal (“Raw Signal”) shows the source waveform distorted by artifacts. The processed signal (“Processed signal”) shows the source waveform with reduced artifacts. In comparison to the raw signal, the processed signal is much cleaner and is similar to the source signal in terms of wave shape and phase.

Fig. 9.

Spectrograms showing the time-frequency representation of phantom source and artifact. The spectrogram of raw data (upper panel, “Raw Data”) shows the time-frequency representation of the source signal distorted with artifacts. The spectrogram of processed data (bottom panel, “Processed Data”) shows the time-frequency representation of source signals after removing artifacts. A comparison of the raw and processed data indicates that removal of artifacts can recover and normalize source signals in the time-frequency domain. The color bar indicates the color coding of power density, which is identical for the two spectrograms.

The accuracy of source localization with raw OPM MEG data from the current dipole phantom was 4.97 $\pm$ 0.34 mm. The accuracy of source localization with processed OPM MEG data was 3.89 $\pm$ 0.072 mm. The source localized with processed OPM MEG data was significantly closer to the ground truth location, as compared to the source localized with raw OPM MEG data (p 0.02). The results indicated that FSSSC significantly improved the accuracy of source localization.

MEG data recorded from human participants revealed that OPMs recorded auditory evoked magnetic fields (AEFs). Fig. 10 shows the AEFs in raw MEG data and in processed MEG data. The strongest response in AEFs in SQUID MEG data was around 100 ms [27]. We also identified the strongest response in AEFs in OPM MEG data around 100 ms. This strongest response in AEFs was named M2 ($>$150 fT). In comparison to the raw MEG data, the processed MEG data also revealed a clear deflection around 70 ms (M1) and another deflection around 150 ms (M3). The M1 and M3 were well-studied in a previous report with a commercial SQUID MEG system [27]. We noted that the strongest AEF deflection was around 100 ms, which was clearly identifiable in raw MEG data. However, the weak M1 and M3 in AEFs (150 fT) were only identified in processed MEG data. This observation indicated that the relatively weak M1 and M3 in AEFs in raw MEG data were overshadowed and distorted by artifacts. Time frequency analyses revealed one component in the raw data. However, time-frequency analyses revealed three components in the FSSC processed data (see Fig. 11). Figs. 10,11 demonstrate that FSSSC reveals these weak AEF components by removing or reducing artifacts.

Fig. 10.

MEG Waveforms from a subject hearing a tone. The waveform of raw data (“Raw Data”) shows one deflection (M2) while M1 is completely obscured by artifacts. M3 is also lightly distorted by artifacts. The waveform of processed data shows three deflections (M1, M2 and M3) clearly.

Fig. 11.

Spectrograms showing the time-frequency representation of auditory evoked magnetic fields (AEFs) and artifact. The spectrogram of raw data (upper panel, “Raw Data”) shows the time-frequency representation of the AEFs distorted with artifacts. The spectrogram of processed data (bottom panel, “Processed Data”) shows the time-frequency representation of AEFs after removing artifacts. A comparison of the raw and processed data indicates that removal of artifacts can reveal three components (M1, M2 and M3). The color bar indicates the color coding of power density, which is identical for the two spectrograms.

Shielding factor analyses revealed that FSSSC removed artifacts by a factor of 8 for M2. The shielding factors for M1 and M2 were undetermined because the waveforms of M1 and M2 in the raw data were not clearly identifiable.