testAvg3

 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 gram = 0 ; % 0=spectra, 1=coherence
0009 chronux = 0; % 1=compare, 0=don't
0010 line = 0;
0011 
0012 nSamples = 4210; 
0013 nChannels = 1; 
0014 nSegments = 1;
0015 movingwin = [50, 50];
0016 
0017 %
0018 % Spectral Parameters
0019 %
0020 params.fpass = [ 0 0.5 ];
0021 params.pad = 2;
0022 params.err = [1 0.85];  % err(1)=0 is no err estimates; err(1)=1 is asymptotic estimates; err(1)=2 is Jacknife
0023 params.trialave = 1;
0024 params.Fs = 1;
0025 
0026 %
0027 % Tapers issues
0028 %
0029 halfBandWidth = 2.5; 
0030 kCap = 2*halfBandWidth - 1;
0031 %params.tapers = [ halfBandWidth, kCap ];
0032 params.tapers = [ halfBandWidth, 1 ];
0033 
0034 %
0035 % Basic checks on inputs
0036 %
0037 if (gram==1) && (nChannels < 2), error( 'Coherence requires at least 2 channels...' ); end 
0038 if (nSegments==1) && (params.err(1)==2), error( 'Jacknife requires more than 1 segment'); end
0039 
0040 %
0041 % Generate segments endpoints randomly
0042 % myrandint is a 3rd party routine (from matlab site)
0043 %
0044 if ~chronux
0045     % Randomly generated segment end points
0046     sMarkers = reshape( sort( myrandint( 2*nSegments, 1, [ 1 : nSamples ], 'noreplace' ) )', 2, nSegments )';
0047 else
0048     % Equal length segments (to compare with Chronux)
0049     sMarkers = ones( nSegments, 2 ); 
0050     sMarkers( :, 2 ) = round( nSamples/2 ); 
0051 end
0052 
0053 %
0054 % Randomly generate the time series
0055 %
0056 fulldata = randn( nSamples, nChannels );  
0057 if line % add line harmonic for testing purposes
0058     f1 = 0.45; a1 = 0.20;
0059     f2 = 0.25; a2 = 0.15;
0060     for c=1:size(fulldata,2)
0061         mx = max(fulldata(:,c));
0062         fulldata(:,c) = fulldata(:,c) + a1*mx*sin( f1*2*pi*[1:size(fulldata,1)]') + a2*mx*sin( f2*2*pi*[1:size(fulldata,1)]') ;
0063     end
0064 end
0065 
0066 %
0067 % Randomly select a few channels
0068 %
0069 chIndices = sort( myrandint( ceil(nChannels/1.5), 1, [ 1 : nChannels ], 'noreplace' ) );
0070 %chIndices = [ 1 : nChannels ];
0071 
0072 %
0073 % Create a data matrix with all the segments aligned one after another
0074 %
0075 totalSegmentLength = sum( sMarkers(:,2) - sMarkers(:,1) + 1 );
0076 data = zeros( totalSegmentLength, length(chIndices) ); % preallocate to ensure contiguous memory
0077 newMarkers(1,1) = 1; 
0078 newMarkers(1,2) = sMarkers(1,2) - sMarkers(1,1) + 1;
0079 data( newMarkers(1,1):newMarkers(1,2), : ) = fulldata( sMarkers(1,1):sMarkers(1,2), chIndices(:));
0080 for sg = 2:size( sMarkers, 1 )
0081     newMarkers(sg,1) = newMarkers(sg-1,2) + 1; 
0082     newMarkers(sg,2) = newMarkers(sg,1) + sMarkers(sg,2) - sMarkers(sg,1);
0083     data( newMarkers(sg,1):newMarkers(sg,2), : ) = fulldata( sMarkers(sg,1):sMarkers(sg,2), chIndices(:));
0084 end
0085 
0086 % To ensure that we check results from array indices beyond 1
0087 if nChannels > 1
0088     ix = sort( myrandint( 1, 2, [1:length(chIndices)], 'noreplace' ) ); % Arbitrarily pick two indices from selected channels for testing results
0089     i1=ix(1); i2=ix(2);
0090     % iC = m + (n-1)*(n-2)/2, for elements of the the coherence matrix, Cmn
0091     iC = ix(1) + (ix(2)-1)*(ix(2)-2)/2;
0092 else
0093     ix = sort( myrandint( 1, 1, [1:length(chIndices)], 'noreplace' ) ); % Arbitrarily pick 1 indices from selected channels for testing results
0094     i1=ix(1);
0095 end
0096 
0097 %
0098 % Power spectrum/spectrogram/coherence/coherogram
0099 %
0100 if gram==0
0101     [ S, f, Serr ] = ueSpectrogram( data, movingwin, params, newMarkers );
0102     figure; plotvector( S(:,i1), f, 'l', squeeze(Serr(:,:,i1)) );
0103     %figure; plot( f, 10*log10( S(:,i1) )); title('Avg. Routine:: Spectrum');
0104 elseif gram==1
0105     [Cmn,Phimn,Smn,Smm,f,I,ConfC,PhiStd,Cerr] = ueCoherence( data, params, newMarkers );
0106     %  C(i,j) = Cmn(:,k) where k = j + (1/2)*(i-1)*(i-2)
0107     figure; plot( f, 10*log10( Cmn(:,iC) ) ); title('Avg. Routine:: Coherence-Magnitude');
0108     %figure; plot( f, phimn(:,iC) ); title('Avg. Routine:: Coherence-Phase');
0109     disp( ['Confidence levelfor C (confC) at (1-p) level: ' num2str( ConfC(iC)) ] );
0110 end
0111 
0112 %
0113 % Use to check against Chronux: only for equal length segments
0114 %
0115 if chronux
0116     cdata = repmat( fulldata( sMarkers(1,1):sMarkers(1,2), chIndices(i1) ), 1, nSegments ); % round(nSamples/2) x nSegments
0117     params.trialave = 1;
0118     if gram==0
0119         [ cS, cf, cSerr ] = mtspectrumc( cdata, params );
0120         figure; plot( cf, 10*log10( cS )); title('Chronux:: Spectrum');
0121         figure; plot( cf, 10*log10(cSerr(1,:)), cf, 10*log10(cSerr(2,:)) ); title('Chronux Error-Bar Computations');
0122         figure; plot( cf, 10*log10( cS ) - 10*log10( S(:,i1) )); title('Error in Spectrum = |New Routines - Chronux|');
0123         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| ');
0124     elseif gram==1
0125         cdata2 = repmat( fulldata( sMarkers(1,1):sMarkers(1,2), chIndices(i2) ), 1, nSegments ); % round(nSamples/2) x nSegments
0126         params.trialave = 1;
0127         
0128         [cC,cphi,cS12,cS1,cS2,cf,cconfC,cphistd,cCerr]=coherencyc( cdata, cdata2, params );
0129         
0130         figure; plot( cf, 10*log10( cC(:,1) ) ); title('Chronux:: Coherence-Magnitude');
0131         figure; plot( cf, 10*log10( cC(:,1) ) - 10*log10( Cmn(:,iC) ) ); title('Error in Coherence = |New Routines - Chronux|');
0132         % Phase may give a problem of 2pi difference... look into it.
0133         figure; plot( cf, cphi(:,1) -  Phimn(:,iC) ); title('Error in Phase = |New Routines - Chronux|');
0134         %
0135         % Note the remaining quantities do not really need to checked since
0136         % coherence = cross-spectrum/power spectra* power spectra, ie C = S12/(S1*S2)
0137         % so unlikely that S12, S1, S2 are incorrect if C is correct.
0138         if 1
0139             figure; plot( cf, 10*log10( cS1(:,1) ) - 10*log10( Smm(:,ix(1)) ) ); title('Error in Power Spectrogram-1 = |New Routines - Chronux|');
0140             figure; plot( cf, 10*log10( cS2(:,1) ) - 10*log10( Smm(:,ix(2)) ) ); title('Error in Power Spectrogram-2 = |New Routines - Chronux|');
0141         end
0142         %
0143         % Error-Bars & Confidence Levels
0144         disp( ['Confidence levelfor C (confC) at (1-p) level: ' num2str( cconfC) ' (Chronux)' ] );
0145         disp( ['Error in confidence level, confC: ' num2str( ConfC(iC) - cconfC ) ] );     
0146         %figure; plot( cf, cphistd(:,1), f, phistd(:,iC) ); title('Phase-Error-Bar Computations');
0147         figure; plot( cf, cphistd(1,:,1) -  PhiStd(1,:,iC) ); title('Error in PhiStd-1');
0148         figure; plot( cf, cphistd(2,:,1) -  PhiStd(2,:,iC) ); title('Error in PhiStd-2');
0149         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|');
0150     end
0151 end
0152 
0153     
0154 
0155 
0156