Updated and made statistics

contrastBeforeAfterAmplitude.m

@@ -13,19 +13,29 @@ afterRaw = load('S6H1_Hann_40_2000feedforward_bins.txt');
 % Our object is to average all the calibrated amplitude rows corresponding
 % to a given SFT number.
 
-beforeL = unique(beforeRaw(:,1));
-afterL = unique(afterRaw(:,1));
+% The text file is zero-indexed
+beforeL = sort(unique(beforeRaw(:,1)));
+afterL = sort(unique(afterRaw(:,1)));
 
 before = ones(size(beforeL));
 after = ones(size(afterL));
+
+
 % Average over all the bins corresponding to a given SFT number,
 % then multiply by sqrt(8/3) because of the Hann window normalization factor.
+% Have to shift index on the right-hand side, to compensate for
+% Matlab's one-indexing vs the C-generated text file's zero-indexing
 for ii = 1:length(before)
-    before(ii) = sqrt(8/3) * mean(beforeRaw(beforeRaw(:,1) == ii, 4));
+    before(ii) = sqrt(8/3) * mean(beforeRaw(beforeRaw(:,1) == ii - 1, 4));
 end
 for ii = 1:length(after)
-    after(ii) = sqrt(8/3) * mean(afterRaw(afterRaw(:,1) == ii, 4));
+    after(ii) = sqrt(8/3) * mean(afterRaw(afterRaw(:,1) == ii - 1, 4));
 end
+% Strip any NaN -- but there should not be any
+before(isnan(before)) = [];
+after(isnan(after)) = [];
+clear beforeRaw
+clear afterRaw
 
 figure(1)
 semilogy(beforeL, before,afterL, after)
@@ -43,6 +53,7 @@ differenceClean = difference;
 differenceClean(isnan(difference)) = [];
 arithmean = mean(differenceClean)*ones(size(difference));
 arithmeanString = horzcat('Arithmetic mean of difference: ', num2str(arithmean(1)));
+length(difference)
 
 figure(2)
 plot(beforeL, difference, beforeL, arithmean)
@@ -57,3 +68,30 @@ ylim([-1e-24 1e-24])
 print('-dpdf', 'HannAmpDiff.pdf')
 print('-dpng', 'HannAmpDiff.png')
 close(2)
+
+% Try some statistics
+HoftHistVector = 1e-24*(50:100);
+beforeHist = hist(before, HoftHistVector);
+afterHist = hist(after, HoftHistVector);
+arithmeanBefore = mean(before);
+arithmeanAfter = mean(after);
+arithmeanDifference = arithmeanBefore - arithmeanAfter;
+harmmeanBefore = harmmean(before);
+harmmeanAfter = harmmean(after);
+harmmeanDifference = harmmeanBefore - harmmeanAfter;
+harmmeanString = 'Harmonic mean: before, after, difference';
+
+figure(3)
+plot(HoftHistVector, beforeHist, HoftHistVector, afterHist)
+xlabel('Hoft')
+ylabel('Histogram count')
+legend('Before feedforward', 'After feedforward')
+title({'Hann windowed calibrated amplitude histogram';...
+    harmmeanString;...
+    num2str(harmmeanBefore);
+    num2str(harmmeanAfter);
+    num2str(harmmeanDifference)})
+grid on
+print('-dpdf', 'HannHist.pdf')
+print('-dpng', 'HannHist.png')
+close(3)

fullspectBoth.m

@@ -5,7 +5,7 @@ hold off
 runlist = ['S6'];
 run = ['S6'];
 %windowlist = ['Hann ';'Tukey']; 
-windowlist = 'Tukey';
+windowlist = 'Hann';
 
 %ifolist = ['H1'; 'L1'];
 ifolist = ['H1'];
@@ -54,10 +54,10 @@ irun = 1;
      color = strtrim(plotcolor(iifo,:));
      fnameroot = sprintf('%s%s',run,ifo);
      fname1 = sprintf('%s_%s_40_2000test.txt',fnameroot,strtrim(windowlist(irun,:)));
-     fnameFull1 = strcat('~gmeadors/2012/10/24/AMPS/', fname1);
+     fnameFull1 = strcat('~gmeadors/2012/11/21/AMPS/', fname1);
      data1 = load(fnameFull1);
      fname2 = sprintf('%s_%s_40_2000feedforward.txt',fnameroot,strtrim(windowlist(irun,:)));
-     fnameFull2 = strcat('~gmeadors/2012/10/24/AMPS/', fname2);
+     fnameFull2 = strcat('~gmeadors/2012/11/21/AMPS/', fname2);
      data2 = load(fnameFull2);
      % Try resampling to smooth random variation.
      % This combines neighboring bins.
@@ -69,11 +69,12 @@ irun = 1;
      resampleFilterSize = 225;
      % Theoretically, the frequency array is linear so it could be resampled
      % using a less sophisticated algorithm. Yet why not try out resample here?
+     % Multiply by sqrt(8/3) for the Hann window normalization
      freq{irun,iifo} = resample(data1(:,1), 1,...
          resampleFactor, resampleFilterSize);
-     amppsd{irun,iifo} = resample(data1(:,5), 1,...
+     amppsd{irun,iifo} = sqrt(8/3)*resample(data1(:,5), 1,...
          resampleFactor, resampleFilterSize);
-     amppsdwt{irun,iifo} = resample(data2(:,5), 1,...
+     amppsdwt{irun,iifo} = sqrt(8/3)*resample(data2(:,5), 1,...
          resampleFactor, resampleFilterSize);
      sprintf('Looping over single-IFO bands...')
      for iband = 1:length(bandlolist)

fullspectDiff.m

@@ -5,7 +5,7 @@ hold off
 runlist = ['S6'];
 run = ['S6'];
 %windowlist = ['Hann ';'Tukey']; 
-windowlist = 'Tukey';
+windowlist = 'Hann';
 
 %ifolist = ['H1'; 'L1'];
 ifolist = ['H1'];
@@ -54,10 +54,10 @@ irun = 1;
      color = strtrim(plotcolor(iifo,:));
      fnameroot = sprintf('%s%s',run,ifo);
      fname1 = sprintf('%s_%s_40_2000test.txt',fnameroot,strtrim(windowlist(irun,:)));
-     fnameFull1 = strcat('~gmeadors/2012/06/29/AMPS/', fname1);
+     fnameFull1 = strcat('~gmeadors/2012/11/21/AMPS/', fname1);
      data1 = load(fnameFull1);
      fname2 = sprintf('%s_%s_40_2000feedforward.txt',fnameroot,strtrim(windowlist(irun,:)));
-     fnameFull2 = strcat('~gmeadors/2012/06/29/AMPS/', fname2);
+     fnameFull2 = strcat('~gmeadors/2012/11/21/AMPS/', fname2);
      data2 = load(fnameFull2);
      % Try resampling to smooth random variation.
      % This combines neighboring bins.
@@ -69,11 +69,12 @@ irun = 1;
      resampleFilterSize = 900;
      % Theoretically, the frequency array is linear so it could be resampled
      % using a less sophisticated algorithm. Yet why not try out resample here?
+     % Because we are using Hann windows, multiply by sqrt(8/3)
      freq{irun,iifo} = resample(data1(:,1), 1,...
          resampleFactor, resampleFilterSize);
-     amppsd{irun,iifo} = resample(data1(:,3) - data2(:,3), 1,...
+     amppsd{irun,iifo} = sqrt(8/3)*resample(data1(:,3) - data2(:,3), 1,...
          resampleFactor, resampleFilterSize);
-     amppsdwt{irun,iifo} = resample(data1(:,5) - data2(:,5), 1,...
+     amppsdwt{irun,iifo} = sqrt(8/3)*resample(data1(:,5) - data2(:,5), 1,...
          resampleFactor, resampleFilterSize);
      sprintf('Looping over single-IFO bands...')
      for iband = 1:length(bandlolist)
@@ -99,7 +100,7 @@ irun = 1;
 	%plot(freqdesign,sensdesign,'color','blue');
         grid
 	legend('Average amplitude PSD','Weighted average amplitude PSD','Location','North')
-	titlestr = sprintf('%s %s Average Spectra (%d-%d Hz)',run,ifo,round(bandlo),round(bandhi));
+	titlestr = sprintf('%s %s Average Spectra, Before - After (%d-%d Hz)',run,ifo,round(bandlo),round(bandhi));
         title(titlestr)
 	fnamepdf = sprintf('%s_%s.pdf',fnameroot,bandname)
         print('-dpdf',fnamepdf)