% TESTACCURACY Test accuracy of Q transform signal recovery % % TESTACCURACY tests the accuracy with which the Q transform recovers % the center time, center frequency, duration, bandwidth, Q, and % signal to noise ratio of sinusoidal Gaussian signal in the presence % of ideal stationary white noise. % % TESTACCURACY repeatedly applies the Q transform to sinusoidal % Gaussian signals injected with random center time, center frequency, % Q, phase, and fixed signal to noise ratio into stationary white % noise. For every injection, TESTACCURACY records the properties of % the most significant template as the measure properties of the % injection. % % As it runs, TESTACCURACY writes the injected center time, center % frequency, Q, and SNR, as well as the measured center time, center % frequency, Q, and SNR, to a log file testaccuracy/testaccuracy.log. % Upon completion, TESTACCURACY also writes results to a MAT file % testaccuracy/testaccuracy.mat and produces the following plots. % % testaccuracy/testaccuracy_1.png % Histrogram of relative time accuracy % % testaccuracy/testaccuracy_2.png % Histogram of relative frequency accuracy % % testaccuracy/testaccuracy_3.png % Histogram of relative duration accuracy % % testaccuracy/testaccuracy_4.png % Histogram of relative bandwidth accuracy % % testaccuracy/testaccuracy_5.png % Histogram of relative Q accuracy % % testaccuracy/testaccuracy_6.png % Histogram of relative SNR accuracy % % See also QTILE and QTRANSFORM. % Shourov K. Chatterji % shourov@ligo.caltech.edu % $Id: testaccuracy.m,v 1.2 2007/04/12 08:49:19 shourov Exp $ % path to qtransform code base path('..', path); % create output directory unix('mkdir -p testaccuracy'); % open log file logFileName = 'testaccuracy/testaccuracy.log'; logFileID = fopen(logFileName, 'w'); % injected snr snr = 10; % number of trials numberOfTrials = 1000; % 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 injectedTimes = []; injectedFrequencies = []; injectedQs = []; injectedPhases = []; injectedSNRs = []; measuredTimes = []; measuredFrequencies = []; measuredQs = []; measuredSNRs = []; % 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 noise noise = normrnd(0, sqrt(tiling.sampleFrequency / 2), 1, length(time)); % construct signal signal = exp(-(time - centerTime).^2 * ... 4 * pi^2 * centerFrequency^2 / q^2) .* ... sin(2 * pi * centerFrequency * (time - centerTime) + phase); % scale signal to desired snr signal = sqrt(snr^2 / 2 - 1) * ... signal / sqrt(sum(signal.^2) / sampleFrequency); % combine signal and noise data = signal + noise; % transform to one-sided frequency domain signal data = fft(data); data = data(1 : length(data) / 2 + 1); % apply q transform calibration = []; outlierFactor = 2.0; transform = qtransform(data, tiling, calibration, outlierFactor); % determine measured signal properties falseEventRate = 1; significants = qthreshold(transform, tiling, 0, falseEventRate, ... [], [], [], [], 1, [], [], [], [], 0); measuredTime = significants{1}.time(1); measuredFrequency = significants{1}.frequency(1); measuredQ = significants{1}.q(1); measuredSNR = sqrt(2 * significants{1}.normalizedEnergy(1)); % report results fprintf(logFileID, ... '%#09.6f %9.3e %9.3e %9.3e %#09.6f %9.3e %9.3e %9.3e\n', ... centerTime, centerFrequency, q, snr, ... measuredTime, measuredFrequency, measuredQ, measuredSNR); % append result to results vectors injectedTimes = [injectedTimes centerTime]; injectedFrequencies = [injectedFrequencies centerFrequency]; injectedQs = [injectedQs q]; injectedPhases = [injectedPhases phase]; injectedSNRs = [injectedSNRs snr]; measuredTimes = [measuredTimes measuredTime]; measuredFrequencies = [measuredFrequencies measuredFrequency]; measuredQs = [measuredQs measuredQ]; measuredSNRs = [measuredSNRs measuredSNR]; % end loop over trials end % save data save testaccuracy/testaccuracy.mat ... tiling numberOfTrials ... injectedTimes injectedFrequencies injectedQs injectedSNRs ... measuredTimes measuredFrequencies measuredQs measuredSNRs % determine durations and bandwidths injectedBandwidths = 2 * sqrt(pi) * injectedFrequencies ./ injectedQs; measuredBandwidths = 2 * sqrt(pi) * measuredFrequencies ./ measuredQs; injectedDurations = 1 ./ injectedBandwidths; measuredDurations = 1 ./ measuredBandwidths; % display relative time accuracy clf; set(gca, 'FontSize', 18); bins = linspace(-0.5, +0.5, 51); relativeError = (measuredTimes - injectedTimes) ./ ... injectedDurations; counts = hist(relativeError, bins); bar(bins, counts); axis([-0.31 +0.31 0 10 * ceil(max(counts) / 10)]); set(gca, 'XTick', -0.3 : 0.1 : +0.3); xlabel('(\tau_{measured} - \tau_{injected}) / \Delta\tau_{injected}'); ylabel('Number of injections'); title(['Relative time error']); aa = axis; handle = text(aa(1) + 0.75 * (aa(2) - aa(1)), ... aa(3) + 0.85 * (aa(4) - aa(3)), ... sprintf(['N = %d\\newline{}' ... '\\mu = %+5.3f\\newline{}' ... '\\sigma = %5.3f'], ... length(relativeError), ... mean(relativeError), ... std(relativeError))); set(handle, 'FontSize', 18); set(handle, 'EdgeColor', 'k'); set(handle, 'LineWidth', 1); set(handle, 'Margin', 5); print -dpng testaccuracy/testaccuracy_1.png unix(['convert -resize 500x500 testaccuracy/testaccuracy_1.png ' ... 'testaccuracy/testaccuracy_1_thumbnail.png']); % display relative frequency accuracy clf; set(gca, 'FontSize', 18); bins = linspace(-0.5, +0.5, 51); relativeError = (measuredFrequencies - injectedFrequencies) ./ ... injectedBandwidths; counts = hist(relativeError, bins); bar(bins, counts); axis([-0.31 +0.31 0 10 * ceil(max(counts) / 10)]); set(gca, 'XTick', -0.3 : 0.1 : +0.3); xlabel('(\phi_{measured} - \phi_{injected}) / \Delta\phi_{injected}'); ylabel('Number of injections'); title(['Relative frequency error']); aa = axis; handle = text(aa(1) + 0.75 * (aa(2) - aa(1)), ... aa(3) + 0.85 * (aa(4) - aa(3)), ... sprintf(['N = %d\\newline{}' ... '\\mu = %+5.3f\\newline{}' ... '\\sigma = %5.3f'], ... length(relativeError), ... mean(relativeError), ... std(relativeError))); set(handle, 'FontSize', 18); set(handle, 'EdgeColor', 'k'); set(handle, 'LineWidth', 1); set(handle, 'Margin', 5); print -dpng testaccuracy/testaccuracy_2.png unix(['convert -resize 500x500 testaccuracy/testaccuracy_2.png ' ... 'testaccuracy/testaccuracy_2_thumbnail.png']); % display relative duration accuracy clf; set(gca, 'FontSize', 18); bins = linspace(-1.5, +1.5, 51); relativeError = log2(measuredDurations ./ injectedDurations); counts = hist(relativeError, bins); bar(bins, counts); axis([-1.05 +1.05 0 10 * ceil(max(counts) / 10)]); set(gca, 'XTick', -1.0 : 0.5 : +1.0); xlabel('log_2(\Delta\tau_{measured} / \Delta\tau_{injected})'); ylabel('Number of injections'); title(['Relative duration error']); aa = axis; handle = text(aa(1) + 0.75 * (aa(2) - aa(1)), ... aa(3) + 0.85 * (aa(4) - aa(3)), ... sprintf(['N = %d\\newline{}' ... '\\mu = %+5.3f\\newline{}' ... '\\sigma = %5.3f'], ... length(relativeError), ... mean(relativeError), ... std(relativeError))); set(handle, 'FontSize', 18); set(handle, 'EdgeColor', 'k'); set(handle, 'LineWidth', 1); set(handle, 'Margin', 5); print -dpng testaccuracy/testaccuracy_3.png unix(['convert -resize 500x500 testaccuracy/testaccuracy_3.png ' ... 'testaccuracy/testaccuracy_3_thumbnail.png']); % display relative bandwidth accuracy clf; set(gca, 'FontSize', 18); bins = linspace(-1.5, +1.5, 51); relativeError = log2(measuredBandwidths ./ injectedBandwidths); counts = hist(relativeError, bins); bar(bins, counts); axis([-1.05 +1.05 0 10 * ceil(max(counts) / 10)]); set(gca, 'XTick', -1.0 : 0.5 : +1.0); xlabel('log_2(\Delta\phi_{measured} / \Delta\phi_{injected})'); ylabel('Number of injections'); title(['Relative bandwidth error']); aa = axis; handle = text(aa(1) + 0.75 * (aa(2) - aa(1)), ... aa(3) + 0.85 * (aa(4) - aa(3)), ... sprintf(['N = %d\\newline{}' ... '\\mu = %+5.3f\\newline{}' ... '\\sigma = %5.3f'], ... length(relativeError), ... mean(relativeError), ... std(relativeError))); set(handle, 'FontSize', 18); set(handle, 'EdgeColor', 'k'); set(handle, 'LineWidth', 1); set(handle, 'Margin', 5); print -dpng testaccuracy/testaccuracy_4.png unix(['convert -resize 500x500 testaccuracy/testaccuracy_4.png ' ... 'testaccuracy/testaccuracy_4_thumbnail.png']); % display relative Q accuracy clf; set(gca, 'FontSize', 18); bins = linspace(-1.5, +1.5, 51); relativeError = log2(measuredQs ./ injectedQs); counts = hist(relativeError, bins); bar(bins, counts); axis([-1.05 +1.05 0 10 * ceil(max(counts) / 10)]); set(gca, 'XTick', -1.0 : 0.5 : +1.0); xlabel('log_2(Q_{measured} / Q_{injected})'); ylabel('Number of injections'); title(['Relative Q error']); aa = axis; handle = text(aa(1) + 0.75 * (aa(2) - aa(1)), ... aa(3) + 0.85 * (aa(4) - aa(3)), ... sprintf(['N = %d\\newline{}' ... '\\mu = %+5.3f\\newline{}' ... '\\sigma = %5.3f'], ... length(relativeError), ... mean(relativeError), ... std(relativeError))); set(handle, 'FontSize', 18); set(handle, 'EdgeColor', 'k'); set(handle, 'LineWidth', 1); set(handle, 'Margin', 5); print -dpng testaccuracy/testaccuracy_5.png unix(['convert -resize 500x500 testaccuracy/testaccuracy_5.png ' ... 'testaccuracy/testaccuracy_5_thumbnail.png']); % display relative SNR accuracy clf; set(gca, 'FontSize', 18); bins = linspace(-1.5, +1.5, 51); relativeError = log2(measuredSNRs ./ injectedSNRs); counts = hist(relativeError, bins); bar(bins, counts); axis([-1.05 +1.05 0 10 * ceil(max(counts) / 10)]); set(gca, 'XTick', -1.0 : 0.5 : +1.0); xlabel('log_2(\rho_{measured} / \rho_{injected})'); ylabel('Number of injections'); title(['Relative SNR error']); aa = axis; handle = text(aa(1) + 0.75 * (aa(2) - aa(1)), ... aa(3) + 0.85 * (aa(4) - aa(3)), ... sprintf(['N = %d\\newline{}' ... '\\mu = %+5.3f\\newline{}' ... '\\sigma = %5.3f'], ... length(relativeError), ... mean(relativeError), ... std(relativeError))); set(handle, 'FontSize', 18); set(handle, 'EdgeColor', 'k'); set(handle, 'LineWidth', 1); set(handle, 'Margin', 5); print -dpng testaccuracy/testaccuracy_6.png unix(['convert -resize 500x500 testaccuracy/testaccuracy_6.png ' ... 'testaccuracy/testaccuracy_6_thumbnail.png']);