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

Generated on Wed 29-Nov-2017 23:17:18 by m2html © 2005