diff --git a/contrastBeforeAfterAmplitude.m b/contrastBeforeAfterAmplitude.m index 81b1b0d1680dbcebf3c70eda38914b51670ba698..14aeb1b63621ac9bde8c79fced51390317814443 100644 --- a/contrastBeforeAfterAmplitude.m +++ b/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) diff --git a/fullspectBoth.m b/fullspectBoth.m index 23b099f7c26d93548c687926d224735152d6e9c2..f7b5057b8a696b039b94316bb741ac9a2d6bbc14 100644 --- a/fullspectBoth.m +++ b/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) diff --git a/fullspectDiff.m b/fullspectDiff.m index edf5d507c201b0b678469f2102202d6c52dc9ec0..aa51f916e3a24f1aa4eee4ef9271d0541e4e0304 100644 --- a/fullspectDiff.m +++ b/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)