Home > NoiseTools > EXAMPLE > example5.m

# example5

## PURPOSE

This example uses DSS to isolate narrowband power near 16 Hz in MEG data.

## SYNOPSIS

This is a script file.

## DESCRIPTION

     This example uses DSS to isolate narrowband power near 16 Hz in MEG data.
     The presence of strong narrowband power near 10Hz would mask this
     component without DSS.

     The best component is dominated by 16 Hz, but the best sensor (ie the sensor
     most strongly containing this component) is dominated mainly by 10 Hz.

     Uses nt_bias_filter(), nt_dss0().

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 % This example uses DSS to isolate narrowband power near 16 Hz in MEG data.
    0002 % The presence of strong narrowband power near 10Hz would mask this
    0003 % component without DSS.
    0004 %
    0005 % The best component is dominated by 16 Hz, but the best sensor (ie the sensor
    0006 % most strongly containing this component) is dominated mainly by 10 Hz.
    0007 %
    0008 % Uses nt_bias_filter(), nt_dss0().
    0009
    0010 clear;
    0011 disp(mfilename);
    0012 help(mfilename)
    0013
    0014 % load data
    0015 FNAME=[fileparts(which('nt_version')), '/DATA/example_data.mat'];
    0016 if ~exist(FNAME);
    0017     disp('file ''../DATA/example_data.mat'' not found, get it at http://cognition.ens.fr/Audition/adc/NoiseTools/DATA/');
    0018     return
    0019 end
    0020 load(FNAME)  % loads 'meg', 'sr'
    0021 % excerpt of data from
    0022 % Duncan, K.K., Hadjipapas, A., Li, S., Kourtzi, Z., Bagshaw, A., Barnes, G., 2009. Identifying
    0023 % spatially overlapping local cortical networks with MEG. Hum. Brain Mapp. 7,
    0024 % 1003?1016.
    0025 % With thanks to authors.
    0026
    0027 % first remove 50 Hz & harmonics (see example4)
    0028 disp('remove 50 Hz & harmonics...');
    0029 [c0,c1]=nt_bias_fft(meg,[50, 100, 150]/sr, 512);
    0030 [todss,pwr0,pwr1]=nt_dss0(c0,c1);
    0031 p1=pwr1./pwr0; % score, proportional to power ratio of 50Hz & harmonics to full band
    0032 z=nt_mmat(meg,todss);
    0033 NREMOVE=20;
    0034 meg=nt_tsr(meg,z(:,1:NREMOVE,:)); % regress out to get clean data
    0035
    0036 % downsample
    0037 DSRATIO=3;
    0038 meg=nt_dsample(meg,DSRATIO);
    0039 sr=sr/DSRATIO;
    0040
    0041
    0042
    0043 % bias filter is second-order resonator
    0044 disp('DSS to isolate 16 Hz components...');
    0045 FPEAK=16; % Hz
    0046 Q=8; % determines width
    0047 [b,a]=nt_filter_peak(FPEAK/(sr/2),Q);
    0048
    0049 % covariance matrices of full band (c0) and filtered (c1)
    0050 [c0,c1]=nt_bias_filter(meg,b,a);
    0051
    0052 % DSS matrix
    0053 [todss,pwr0,pwr1]=nt_dss0(c0,c1);
    0054 p1=pwr1./pwr0; % score, proportional to power ratio of 50Hz & harmonics to full band
    0055
    0056 % DSS components
    0057 z=nt_mmat(meg,todss);
    0058
    0059
    0060 % plot bias score
    0061 figure(1); clf; set(gcf,'color', [1 1 1]);
    0062 plot(p1, '.-'); xlabel('component'); ylabel('score'); title('DSS score');
    0063
    0064
    0065 % plot spectra of DSS components
    0066 figure(2); clf; set(gcf,'color', [1 1 1]);
    0067 nt_spect_plot2(nt_normcol(z(:,1:30,:)),512,[],[],sr);
    0068 title('spectra of first 30 DSS components'); ylabel('component')
    0069
    0070 % plot spectra of best DSS component and best sensor
    0071 figure(3); clf; set(gcf,'color', [1 1 1]);
    0072 nt_spect_plot(nt_normcol(z(:,1,:)),512,[],[],sr);
    0073 hold on;
    0074 [~,idx]=max(abs(nt_xcov(z(:,1,:),nt_normcol(meg))));
    0075 nt_spect_plot(nt_normcol(meg(:,idx,:)),512,[],[],sr);
    0076 nt_linecolors([], [3 1]);
    0077 legend('best component', 'best sensor'); legend boxoff; axis tight
    0078 title('bias to 16 Hz')
    0079
    0080