Home > NoiseTools > EXAMPLE > example4.m

# example4

## PURPOSE

This example uses DSS to remove power-line noise from MEG data.

## SYNOPSIS

This is a script file.

## DESCRIPTION

     This example uses DSS to remove power-line noise from MEG data.
     Finds components that are most dominated by 50 Hz and
     harmonics, regresses them out to obtain clean data.

     Uses nt_bias_fft(), nt_dss0(), nt_tsr().

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 % This example uses DSS to remove power-line noise from MEG data.
    0002 % Finds components that are most dominated by 50 Hz and
    0003 % harmonics, regresses them out to obtain clean data.
    0004 %
    0005 % Uses nt_bias_fft(), nt_dss0(), nt_tsr().
    0006
    0007 clear;
    0008 disp(mfilename);
    0009 help(mfilename)
    0010
    0011 % load data
    0012 FNAME=[fileparts(which('nt_version')), '/DATA/example_data.mat'];
    0013 if ~exist(FNAME);
    0014     disp(['file ', FNAME, ' not found, get it at http://cognition.ens.fr/Audition/adc/NoiseTools/DATA/']);
    0015     return
    0016 end
    0017 load(FNAME)  % loads 'meg', 'sr'
    0018 % excerpt of data from
    0019 % Duncan, K.K., Hadjipapas, A., Li, S., Kourtzi, Z., Bagshaw, A., Barnes, G., 2009. Identifying
    0020 % spatially overlapping local cortical networks with MEG. Hum. Brain Mapp. 7,
    0021 % 1003?1016.
    0022 % With thanks to authors.
    0023
    0024 % use DSS to isolate 50 Hz & harmonics
    0025 disp('50 Hz & harmonics DSS...');
    0026
    0027 % covariance matrices of full band (c0) and filtered to 50 Hz & harmonics (c1)
    0028 [c0,c1]=nt_bias_fft(meg,[50, 100, 150]/sr, 512);
    0029
    0030 % DSS matrix
    0031 [todss,pwr0,pwr1]=nt_dss0(c0,c1);
    0032 p1=pwr1./pwr0; % score, proportional to power ratio of 50Hz & harmonics to full band
    0033
    0034 % DSS components
    0035 z=nt_mmat(meg,todss);
    0036
    0037 % first components are most dominated by 50Hz & harmonics
    0038 NREMOVE=20;
    0039 clean=nt_tsr(meg,z(:,1:NREMOVE,:)); % regress them out
    0040
    0041 % plot bias score
    0042 figure(1); clf; set(gcf,'color', [1 1 1]);
    0043 plot(p1, '.-'); xlabel('component'); ylabel('score'); title('DSS score');
    0044
    0045 % plot spectra of DSS components
    0046 figure(2); clf; set(gcf,'color', [1 1 1]);
    0047 nt_spect_plot2(nt_normcol(z(:,1:30,:)),512,[],[],sr);
    0048 title('spectra of first 30 DSS components'); ylabel('component')
    0049
    0050 % plot spectra of data before and after removal of power line components
    0051 figure(3); clf; set(gcf,'color', [1 1 1]);
    0052 subplot 121;
    0053 nt_spect_plot(meg,1024,[],[],sr);
    0054 set(gca,'ygrid','on');
    0055 hold on
    0056 nt_spect_plot(clean,1024,[],[],sr);
    0057 nt_linecolors([],[3 1]);
    0058 title('power spectra, average over channels');
    0059 legend('before','after'); legend boxoff
    0060 set(gca,'ygrid','on');
    0061 subplot 122
    0062 nt_spect_plot(meg-clean,1024,[],[],sr);
    0063 title('noise power (removed)');
    0064 set(gca,'ygrid','on');
    0065
    0066
    0067