testAvg4

 This is a calling routine to test & check out the power spectrum &
 spectrogram routines for unequal segment lengths. In addition, use it 
 to compare with Chronux routines when segments are of equal length.

0001 %
0002 % This is a calling routine to test & check out the power spectrum &
0003 % spectrogram routines for unequal segment lengths. In addition, use it
0004 % to compare with Chronux routines when segments are of equal length.
0005 %
0006 clear all;
0007 
0008 if 0
0009     dir = 'G:\ravi\Chrowser\Pass~ Tioga_0e927741-9673-46e5-9050-ca1d7541bf22\';
0010     xfile = 'Pass~ Tioga_0e927741-9673-46e5-9050-ca1d7541bf22'
0011     %dir = 'G:\ravi\Chrowser\sample~ data_8ef647e3-e5ea-43a6-8c69-fb848b8db7c2\';
0012     %xfile = 'sample~ data_8ef647e3-e5ea-43a6-8c69-fb848b8db7c2'
0013 else
0014     dir = 'Z:\xltekRawData\Wallis~ Terry_c3f44891-afa7-4fa7-a643-55c772a05241\'
0015     xfile = 'Wallis~ Terry_c3f44891-afa7-4fa7-a643-55c772a05241'
0016 end
0017 
0018 % Get header info
0019 % Channels are labelled from C1 through C127 and ''
0020 % total of 128 channels
0021 hdr = eegMex( dir, xfile);
0022 
0023 
0024 gram = 1 ; % 0=spectra, 1=coherence
0025 chronux = 0 ; % 0=no comparison with Chronux; 1=compare with chronux
0026 
0027 %nSamples = 4210;
0028 nChannels = 2; 
0029 nSegments = 1;
0030 movingwin = [25, 25];
0031 
0032 %
0033 % Spectral Parameters
0034 %
0035 params.fpass = [ 0 0.5 ];
0036 params.pad = 2;
0037 params.err = [2 0.05];  % err(1)=0 is no err estimates; err(1)=1 is asymptotic estimates; err(1)=2 is Jacknife
0038 params.trialave = 1;
0039 params.Fs = 1;
0040 
0041 %
0042 % Tapers issues
0043 %
0044 halfBandWidth = 2.5; 
0045 kCap = 2*halfBandWidth - 1;
0046 %params.tapers = [ halfBandWidth, kCap ];
0047 params.tapers = [ halfBandWidth, 2 ];
0048 
0049 %
0050 % Basic checks on inputs
0051 %
0052 if (gram==1) && (nChannels < 2), error( 'Coherence requires at least 2 channels...' ); end 
0053 %if (nSegments==1) && (params.err(1)==2), error( 'Jacknife requires more than 1 segment'); end
0054 
0055 %
0056 % Generate segments endpoints randomly
0057 % myrandint is a 3rd party routine (from matlab site)
0058 %
0059 % Randomly generated segment end points
0060 sMarkers = reshape( sort( myrandint( 2*nSegments, 1, [ ceil(hdr.nSamples/500) : ceil(hdr.nSamples/50) ], 'noreplace' ) )', 2, nSegments )';
0061 %sMarkers = [ ceil(hdr.nSamples/80), ceil(hdr.nSamples/65) ];
0062 
0063 %
0064 % Randomly select a few channels
0065 %
0066 if ~chronux
0067     %chIndices = sort( myrandint( nChannels, 1, [ round(hdr.nChans/4) : round(3*hdr.nChans/4) ], 'noreplace' ) );
0068     chIndices = [ 3 : 3+nChannels-1 ];
0069 else
0070     %chIndices = [ 10 : 10+nChannels-1 ];
0071     chIndices = [ 3, 7 ];
0072 end
0073 
0074 %
0075 % Randomly generate the time series
0076 %
0077 fulldata = eegMex( dir, xfile, chIndices, [ 1 hdr.nSamples/50 1 ] );
0078 mDiscardBits = 0;
0079 conversionFactor = ( 8711 / (2^21 - 0.5) ) * 2^mDiscardBits;
0080 fulldata{:} = fulldata{:} * conversionFactor;
0081 
0082 %
0083 % Create a data matrix with all the segments aligned one after another
0084 %
0085 totalSegmentLength = sum( sMarkers(:,2) - sMarkers(:,1) + 1 );
0086 data = zeros( totalSegmentLength, length(chIndices) ); % preallocate to ensure contiguous memory
0087 newMarkers(1,1) = 1; 
0088 newMarkers(1,2) = sMarkers(1,2) - sMarkers(1,1) + 1;
0089 data( newMarkers(1,1):newMarkers(1,2), : ) = detrend( fulldata{1}( sMarkers(1,1):sMarkers(1,2), :) );
0090 for sg = 2:size( sMarkers, 1 )
0091     newMarkers(sg,1) = newMarkers(sg-1,2) + 1; 
0092     newMarkers(sg,2) = newMarkers(sg,1) + sMarkers(sg,2) - sMarkers(sg,1);
0093     data( newMarkers(sg,1):newMarkers(sg,2), : ) = detrend( fulldata{1}( sMarkers(sg,1):sMarkers(sg,2), :) );
0094 end
0095 
0096 % To ensure that we check results from array indices beyond 1
0097 if nChannels > 1
0098     ix = sort( myrandint( 1, 2, [1:length(chIndices)], 'noreplace' ) ); % Arbitrarily pick two indices from selected channels for testing results
0099     i1=ix(1); i2=ix(2);
0100     % iC = m + (n-1)*(n-2)/2, for elements of the the coherence matrix, Cmn
0101     iC = ix(1) + (ix(2)-1)*(ix(2)-2)/2;
0102 else
0103     ix = sort( myrandint( 1, 1, [1:length(chIndices)], 'noreplace' ) ); % Arbitrarily pick 1 indices from selected channels for testing results
0104     i1=ix(1);
0105 end
0106 
0107 %
0108 % Power spectrum/spectrogram/coherence/coherogram
0109 %
0110 if gram==0
0111     [ S, f, Serr ] = avgSpectrum( data, movingwin, params, newMarkers );
0112     figure; plot( f, 10*log10( S(:,i1) ), 'k', f, 10*log10( Serr(2,:,i1) ), 'g--', f, 10*log10( Serr(1,:,i1)), 'g--' ); title('Avg. Routine:: Spectrum');
0113     %figure; plot( f, 10*log10( S(:,i1) )); title('Avg. Routine:: Spectrum');
0114 elseif gram==1
0115     [Cmn,Phimn,Smn,Smm,f,ConfC,PhiStd,Cerr] = avgCoherence( data, movingwin, params, newMarkers );
0116     %  C(i,j) = Cmn(:,k) where k = j + (1/2)*(i-1)*(i-2)
0117     figure; plot( f, Cmn(:,iC), 'k', f, Cerr(2,:,iC), 'g--', f, Cerr(1,:,iC), 'g--' ); 
0118     title('Avg. Routine:: Coherence'); ylim([0 1])
0119     %figure; plot( f, 10*log10( Cmn(:,iC) ) ); title('Avg. Routine:: Coherence-Magnitude');
0120     %figure; plot( f, phimn(:,iC) ); title('Avg. Routine:: Coherence-Phase');
0121     disp( ['Confidence level for C (confC) at (1-p) level: ' num2str( ConfC(iC)) ] );
0122 end
0123 
0124 
0125 
0126 %
0127 % Use to check against Chronux: only for equal length segments
0128 %
0129 if chronux
0130 
0131     win = floor( newMarkers(1,2) / movingwin(1) );
0132     newMarkers(1,2) = newMarkers(1,2) - mod( newMarkers(1,2), win );
0133     cdata = data( [1:newMarkers(1,2)], i1 );
0134     cdata = detrend( reshape( cdata, [ newMarkers(1,2)/win, win ] ) );
0135     cdata2 = data( [1:newMarkers(1,2)], i2 );
0136     cdata2 = detrend( reshape( cdata2, [ newMarkers(1,2)/win, win ] ) );
0137     params.trialave = 1;
0138     if gram==0
0139         [ cS, cf, cSerr ] = mtspectrumc( cdata, params );
0140         figure; plotvector( cS, cf, 'l', cSerr );
0141         %figure; plot( cf, 10*log10( cS )); title('Chronux:: Spectrum');
0142         figure; plot( cf, 10*log10(cSerr(1,:)), cf, 10*log10(cSerr(2,:)) ); title('Chronux Error-Bar Computations');
0143         figure; plot( cf, 10*log10( cS ) - 10*log10( S(:,i1) )); title('Error in Spectrum = |New Routines - Chronux|');
0144         figure; plot( cf, 10*log10(cSerr(1,:)) - 10*log10(Serr(1,:,i1)), cf, 10*log10(cSerr(2,:)) - 10*log10(Serr(2,:,i1)) );title('Error in Error-Bar Computations = |New Routines - Chronux| ');
0145     elseif gram==1
0146         
0147         [cC,cphi,cS12,cS1,cS2,cf,cconfC,cphistd,cCerr]=coherencyc( cdata, cdata2, params );
0148         
0149         %figure; plotvector( cC(:,1), cf, 'n', cCerr );
0150         figure; plot( cf, cC(:,iC), 'k', cf, cCerr(2,:,iC), 'g--', cf, cCerr(1,:,iC), 'g--' );
0151         title('Chronux:: Coherence'); ylim([0 1])
0152         %figure; plot( cf, 10*log10( cC(:,1) ) ); title('Chronux:: Coherence-Magnitude');
0153         figure; plot( cf, 10*log10( cC(:,1) ) - 10*log10( Cmn(:,iC) ) ); title('Error in Coherence = |New Routines - Chronux|');
0154         % Phase may give a problem of 2pi difference... look into it.
0155         figure; plot( cf, cphi(:,1) -  Phimn(:,iC) ); title('Error in Phase = |New Routines - Chronux|');
0156         %
0157         % Note the remaining quantities do not really need to checked since
0158         % coherence = cross-spectrum/power spectra* power spectra, ie C = S12/(S1*S2)
0159         % so unlikely that S12, S1, S2 are incorrect if C is correct.
0160         if 1
0161             figure; plot( cf, 10*log10( cS1(:,1) ) - 10*log10( Smm(:,ix(1)) ) ); title('Error in Power Spectrogram-1 = |New Routines - Chronux|');
0162             figure; plot( cf, 10*log10( cS2(:,1) ) - 10*log10( Smm(:,ix(2)) ) ); title('Error in Power Spectrogram-2 = |New Routines - Chronux|');
0163         end
0164         %
0165         % Error-Bars & Confidence Levels
0166         disp( ['Confidence levelfor C (confC) at (1-p) level: ' num2str( cconfC) ' (Chronux)' ] );
0167         disp( ['Error in confidence level, confC: ' num2str( ConfC(iC) - cconfC ) ] );     
0168         %figure; plot( cf, cphistd(:,1), f, phistd(:,iC) ); title('Phase-Error-Bar Computations');
0169         figure; plot( cf, cphistd(:,1) -  PhiStd(:,iC) ); title('Error in PhiStd-1');
0170         figure; plot( cf, cphistd(:,1) -  PhiStd(:,iC) ); title('Error in PhiStd-2');
0171         figure; plot( cf, abs(cCerr(1,:,1) -  Cerr(1,:,iC)), cf, abs(cCerr(2,:,1) -  Cerr(2,:,iC))  ); title('Error in Abs(Coherence)-Error-Bar Computations = |New Routines - Chronux|');
0172     end
0173 end
0174 
0175     
0176 
0177 
0178