% TESTMISMATCH Test template placeent algorithm and Q transform % % TESTMISMATCH tests that the template bank generated by QTILE % sufficiently covers the target signal space. % % TESTMISMATCH first tiles the signal space using QTILE and then % repeatedly applies the Q transform on sinusoidal Gaussian signals % with random center time, center frequency, Q, and phase within the % the targeted signal space. All of the signals are injected with % unity energy and without any detector noise. For every injection, % TESTMISMATCH records the fractional error in recovering the injected % energy. % % As it runs, TESTMISMATCH writes the signal center time, center % frequency, Q, and measured fractional energy loss to a log file % testmismatch/testmismatch.log. TESTMISMATCH can also produce % intermediate plots of the cumulate fractional energy loss % distribution after each injection if the hardcoded watchFlag is set % to a non-zero value. Upon completion, TESTMISMATCH also writes % results to a MAT file testmismatch/testmismatch.mat. % % Upon completion, TESTMISMATCH also produces the following four % plots. % % testmismatch/testmismatch_1.png % Cumulative distribution of fractional energy loss % % testmismatch/testmismatch_2.png % Q-frequency scatter plot of fractional energy loss % % testmismatch/testmismatch_3.png % Time-frequency scatter plot of fractional energy loss % % testmismatch/testmismatch_4.png % Time-Q scatter plot of fractional energy loss % % See also QTILE and QTRANSFORM. % Shourov K. Chatterji % shourov@ligo.caltech.edu % $Id: testmismatch.m,v 1.2 2007/04/11 11:30:35 shourov Exp $ % path to qtransform code base path('..', path); % create output directory unix('mkdir -p testmismatch'); % open log file logFileName = 'testmismatch/testmismatch.log'; logFileID = fopen(logFileName, 'w'); % number of trials numberOfTrials = 1000; % boolean flag to produce intermediate plots watchFlag = 0; % signal space parameters timeRange = 64; qRange = [sqrt(11) 100]; frequencyRange = [0 Inf]; sampleFrequency = 4096; % desired tiling accuracy maximumMismatch = 0.2; % tile the signal space tiling = qtile(timeRange, qRange, frequencyRange, sampleFrequency, ... maximumMismatch); % time vector for signal construction time = 0 : 1 / sampleFrequency : timeRange - 1 / sampleFrequency; % initialize random number generator rand('state', sum(100 * clock)); randn('state', sum(100 * clock)); % initialize results vectors centerTimes = []; centerFrequencies = []; qs = []; phases = []; energyLosses = []; % begin loop over trials for trialNumber = 1 : numberOfTrials, % select random center time centerTime = unifrnd(0 + tiling.transientDuration, ... tiling.duration - tiling.transientDuration, ... 1, 1);; % select random quality factor q = 2.^unifrnd(log2(tiling.minimumQ), ... log2(tiling.maximumQ)); % select random center frequency minimumFrequency = 50 * q / (2 * pi * tiling.duration); maximumFrequency = tiling.sampleFrequency / (2 * (1 + sqrt(11) / q)); centerFrequency = 2.^unifrnd(log2(max(minimumFrequency, ... tiling.highPassCutoff)), ... log2(min(maximumFrequency, ... tiling.lowPassCutoff))); % select random phase phase = unifrnd(0, 2 * pi, 1, 1); % construct signal signal = exp(-(time - centerTime).^2 * ... 4 * pi^2 * centerFrequency^2 / q^2) .* ... sin(2 * pi * centerFrequency * (time - centerTime) + phase); % normalize signal to unity amplitude and energy signal = signal / sqrt(sum(signal.^2) / sampleFrequency); % transform to one-sided frequency domain signal signal = fft(signal); signal = signal(1 : length(signal) / 2 + 1); % apply q transform calibration = []; outlierFactor = 2.0; transform = qtransform(signal, tiling, calibration, outlierFactor); % determine measured signal energy measuredEnergy = 0; for plane = 1 : tiling.numberOfPlanes, for row = 1 : tiling.planes{plane}.numberOfRows, measuredEnergy = max([measuredEnergy ... transform{1}.planes{plane}.rows{row}.energies]); end end % determine energy loss energyLoss = 1 - measuredEnergy; % report results fprintf(logFileID, '%#09.6f %9.3e %9.3e %9.3e\n', ... centerTime, centerFrequency, q, energyLoss); % append result to results vectors centerTimes = [centerTimes centerTime]; centerFrequencies = [centerFrequencies centerFrequency]; qs = [qs q]; energyLosses = [energyLosses energyLoss]; % display cumulative distribution of energy loss if watchFlag ~= 0, if trialNumber > 1, threshold = sort(energyLosses); fraction = (length(threshold) : -1 : 1) / length(threshold); clf; set(gca, 'FontSize', 18); handle = semilogy(threshold, fraction, 'b.-'); axis([0 tiling.maximumMismatch ... 10^floor(log10(1/length(threshold))) 1]); grid on; set(handle, 'MarkerSize', 20); set(handle, 'LineWidth', 2); xlabel('Fractional energy loss threshold'); ylabel('Fraction of trials exceeding threshold'); title('Cumulative distribution of energy loss'); drawnow; end end % end loop over trials end % save data save testmismatch/testmismatch.mat ... tiling numberOfTrials centerTimes centerFrequencies qs energyLosses % order by increasing energy loss [ignore, indices] = sort(energyLosses); centerTimes = centerTimes(indices); centerFrequencies = centerFrequencies(indices); qs = qs(indices); energyLosses = energyLosses(indices); % display cumulative distribution of energy loss threshold = sort(energyLosses); fraction = (length(threshold) : -1 : 1) / length(threshold); clf; set(gca, 'FontSize', 18); handle = semilogy(threshold, fraction, 'b.-'); axis([0 tiling.maximumMismatch ... 10^floor(log10(1/length(threshold))) 1]); grid on; set(handle, 'MarkerSize', 20); set(handle, 'LineWidth', 2); xlabel('Fractional energy loss threshold'); ylabel('Fraction of trials exceeding threshold'); title('Cumulative distribution of energy loss'); drawnow; print -dpng testmismatch/testmismatch_1.png unix(['convert -resize 500x500 testmismatch/testmismatch_1.png ' ... 'testmismatch/testmismatch_1_thumbnail.png']); % saturate colormap at maximumMismatch energyLosses(find(energyLosses > tiling.maximumMismatch)) = ... tiling.maximumMismatch; centerTimes = [0 centerTimes 0]; centerFrequencies = [0 centerFrequencies 0]; qs = [0 qs 0]; energyLosses = [0 energyLosses tiling.maximumMismatch]; % display q-frequency distirbution of energy loss clf; set(gca, 'FontSize', 18); handle = scatter(qs, centerFrequencies, 16, energyLosses, 'filled'); set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); axis([tiling.minimumQ tiling.maximumQ ... tiling.highPassCutoff tiling.lowPassCutoff]); set(gca, 'XTick', ... 2.^(floor(log2(tiling.minimumQ)) : ... ceil(log2(tiling.maximumQ)))); set(gca, 'XTickLabel', ... 2.^(floor(log2(tiling.minimumQ)) : ... ceil(log2(tiling.maximumQ)))); set(gca, 'YTick', ... 2.^(floor(log2(tiling.highPassCutoff)) : ... ceil(log2(tiling.lowPassCutoff)))); set(gca, 'YTickLabel', ... 2.^(floor(log2(tiling.highPassCutoff)) : ... ceil(log2(tiling.lowPassCutoff)))); set(gca, 'XGrid', 'on'); set(gca, 'YGrid', 'on'); set(gca, 'XMinorTick', 'off'); set(gca, 'YMinorTick', 'off'); set(gca, 'XMinorGrid', 'off'); set(gca, 'YMinorGrid', 'off'); xlabel('Q []'); ylabel('Center frequency [Hz]'); title('Q-frequency distribution of energy loss'); handle = colorbar; set(handle, 'FontSize', 18); set(get(handle, 'YLabel'), 'String', 'Fractional energy loss'); set(get(handle, 'YLabel'), 'FontSize', 18); drawnow; print -dpng testmismatch/testmismatch_2.png unix(['convert -resize 500x500 testmismatch/testmismatch_2.png ' ... 'testmismatch/testmismatch_2_thumbnail.png']); % display time-frequency distirbution of energy loss clf; set(gca, 'FontSize', 18); handle = scatter(centerTimes, centerFrequencies, 16, energyLosses, ... 'filled'); set(gca, 'XScale', 'linear'); set(gca, 'YScale', 'log'); axis([tiling.transientDuration tiling.duration - tiling.transientDuration ... ... tiling.highPassCutoff tiling.lowPassCutoff]); set(gca, 'YTick', ... 2.^(floor(log2(tiling.highPassCutoff)) : ... ceil(log2(tiling.lowPassCutoff)))); set(gca, 'YTickLabel', ... 2.^(floor(log2(tiling.highPassCutoff)) : ... ceil(log2(tiling.lowPassCutoff)))); set(gca, 'XGrid', 'on'); set(gca, 'YGrid', 'on'); set(gca, 'XMinorTick', 'off'); set(gca, 'YMinorTick', 'off'); set(gca, 'XMinorGrid', 'off'); set(gca, 'YMinorGrid', 'off'); xlabel('Center time [seconds]'); ylabel('Center frequency [Hz]'); title('Time-frequency distribution of energy loss'); handle = colorbar; set(handle, 'FontSize', 18); set(get(handle, 'YLabel'), 'String', 'Fractional energy loss'); set(get(handle, 'YLabel'), 'FontSize', 18); drawnow; print -dpng testmismatch/testmismatch_3.png unix(['convert -resize 500x500 testmismatch/testmismatch_3.png ' ... 'testmismatch/testmismatch_3_thumbnail.png']); % display time-q distirbution of energy loss clf; set(gca, 'FontSize', 18); handle = scatter(centerTimes, qs, 16, energyLosses, 'filled'); set(gca, 'XScale', 'linear'); set(gca, 'YScale', 'log'); axis([tiling.transientDuration tiling.duration - tiling.transientDuration ... ... tiling.minimumQ tiling.maximumQ]); set(gca, 'YTick', ... 2.^(floor(log2(tiling.minimumQ)) : ... ceil(log2(tiling.maximumQ)))); set(gca, 'YTickLabel', ... 2.^(floor(log2(tiling.minimumQ)) : ... ceil(log2(tiling.maximumQ)))); set(gca, 'XGrid', 'on'); set(gca, 'YGrid', 'on'); set(gca, 'XMinorTick', 'off'); set(gca, 'YMinorTick', 'off'); set(gca, 'XMinorGrid', 'off'); set(gca, 'YMinorGrid', 'off'); xlabel('Center time [seconds]'); ylabel('Q []'); title('Time-Q distribution of energy loss'); handle = colorbar; set(handle, 'FontSize', 18); set(get(handle, 'YLabel'), 'String', 'Fractional energy loss'); set(get(handle, 'YLabel'), 'FontSize', 18); drawnow; print -dpng testmismatch/testmismatch_4.png unix(['convert -resize 500x500 testmismatch/testmismatch_4.png ' ... 'testmismatch/testmismatch_4_thumbnail.png']);