Guest User

Untitled

a guest
Oct 3rd, 2017
110
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2. close all
  3.  
  4. % SETUP PARAMETERS
  5. nSegments = 1000; % Specify number of time iterations
  6. distributionParameter = 0.095; %(for exponential to be similar to simulation)
  7.  
  8. %Bandwidth
  9. BW_PD = 200;
  10.  
  11. %Photodetector Noise amplitude
  12. noiseAmplitude = 0.0001: 0.03: 20;
  13. noiseAmplitude = sqrt(noiseAmplitude.*1e-3);
  14.  
  15. start_windowSize = 150;
  16.  
  17. delay = 60;
  18.  
  19. showcorr = 0;
  20.  
  21. % OPTIONS
  22. % For correlation to use two times the size of the reference (optional)
  23. referenceTwiceSize = 1;
  24. % filter the noise to have same bandwidth as the trace
  25. filterNoise = 1;
  26.  
  27.  
  28. %% Generate trace
  29.  
  30. coarseDelay = fix(delay);
  31. fineDelay = delay*100 - coarseDelay*100;
  32.  
  33. save_var_corr_END = [];
  34. save_var_LS_END = [];
  35.  
  36. % Generate the trace
  37.  
  38. rawData = detrend(normrnd(0,distributionParameter,1000,1),'constant'); % Generate 9000 point trace, force zero mean, and filter it with a gauss shape.
  39. rawData_f = fft(rawData);
  40.  
  41. f = -length(rawData_f)/2:length(rawData_f)/2-1;
  42. gauss_pass = exp(-f.^2/(2*BW_PD^2));
  43. rawData = ifft(ifftshift(fftshift(rawData_f)'.*gauss_pass)).'; % Filtered
  44. rawData=repmat(rawData.',nSegments,1);
  45.  
  46. % Delay the trace
  47.  
  48. fineDelay = fineDelay + 1;
  49. rs_pdTrace = resample(rawData(1:end,:).',100,1).'; %Resample the trace, each point becomes 100.
  50. rawData(1,:) = rs_pdTrace(1,1:100:end); % Save the reference trace
  51. rawData(2:end,:) = rs_pdTrace(2:end,fineDelay:100:end); % Generate the other traces.
  52.  
  53.  
  54. % START SIMULATION
  55. startPosTrace = 0.52 + 0.02*rand; % Pick a window at random (Position is defined so that it finds a perturbation in case of generated trace w/ perturbation).
  56.  
  57. % Allocate storage for the simulated variances.
  58. save_var_corr = zeros(1,length(noiseAmplitude));
  59. save_var_LS = zeros(1,length(noiseAmplitude));
  60. save_snr = zeros(1,length(noiseAmplitude));
  61. idx = 0;
  62.  
  63. % save calculated variances
  64. saveCalcVar_DC = [];
  65. saveCalcVar_LS = [];
  66.  
  67. for nAmp = noiseAmplitude
  68. % Variables to calculate the variances
  69. x_dc = []; y_dc = []; z_dc = [];
  70. x_ls = []; y_ls = []; z_ls = [];
  71.  
  72. % Run it for increasing photodetector noise.
  73. idx = idx+1;
  74. noise = nAmp * randn(size(rawData)); % Generate a matrix of normal distributed values with the amplitude of the noise.
  75.  
  76. if filterNoise == 1
  77. for seg = 1:size(noise,1) % Filter each of the noises using the same bandpass filter as the signal.
  78. noise_f = fft(noise(seg,:));
  79. f = -length(noise(seg,:))/2:length(noise(seg,:))/2-1;
  80. gauss_pass = exp(-f.^2/(2*BW_PD^2));
  81.  
  82. noise_f_filtered = fftshift(noise_f).*gauss_pass;
  83.  
  84. noise(seg,:) = ifft(ifftshift(noise_f_filtered));
  85.  
  86. end
  87. end
  88.  
  89.  
  90. SNR = mean((var(rawData.')) ./ (var(noise.'))); % Calculate SNR as the variance of signal over variance of noise.
  91. save_snr(idx) = SNR;
  92.  
  93. traceMatrix = rawData + noise;
  94.  
  95. if plot_traces == 1
  96. figure(1), subplot(2,2,1),plot(rawData(1,:)) , hold on
  97. figure(1),subplot(2,2,1), plot(noise(1,:)),
  98. hold off
  99. figure(1),subplot(2,2,2), plot(traceMatrix(1,:))
  100. figure(1),subplot(2,2,3), plot(abs(fft(noise(1,:))));
  101. figure(1),subplot(2,2,4),plot(abs(fft(traceMatrix(1,:))));
  102. drawnow;
  103. end
  104.  
  105. windowSize = start_windowSize;
  106. pd_trace = traceMatrix;
  107. windowPos = fix(length(pd_trace(1,:))*startPosTrace);
  108. timeIncrement = 1; % Time increment (only relevant in case there is some changing perturbation to our signal)
  109.  
  110.  
  111.  
  112.  
  113. % Correlation
  114. if referenceTwiceSize
  115. win_t0 = detrend(pd_trace(1,windowPos-windowSize/2:windowPos+windowSize+windowSize/2),'constant');
  116. else
  117. win_t0 = detrend(pd_trace(1,windowPos:windowPos+windowSize),'constant');
  118. end
  119.  
  120. save_lags = [];
  121.  
  122. for i = 2:timeIncrement:nSegments-1 % TDE by direct correlation
  123. win = detrend(pd_trace(i,windowPos+coarseDelay:windowPos+coarseDelay+windowSize),... % Make sure that the window is zero-mean.
  124. 'constant');
  125. cor = xcorr(win_t0,win); % Correlate window_t1 with window_t0
  126.  
  127. [~,max_cor]=max(cor); % Find position of maximum
  128.  
  129. %Save lag values
  130. x_dc = [x_dc cor(max_cor+1)];
  131. y_dc = [y_dc cor(max_cor-1)];
  132. z_dc = [z_dc cor(max_cor)];
  133.  
  134. lag = quadFit_lag(cor,max_cor); % Parabolic fit using nearest neighbours.
  135.  
  136. save_lags = [save_lags lag]; % Save the lag
  137.  
  138.  
  139. if showcorr
  140. figure(2),
  141. if referenceTwiceSize
  142. subplot(2,1,1), plot(fix(-length(cor)/2):1:fix(length(cor)/2),cor,'o-'), hold on;
  143. else
  144. subplot(2,1,1), plot(fix(-windowSize):1:fix(windowSize),cor,'o-'), hold on;
  145. end
  146. end
  147.  
  148. end
  149.  
  150. cov_xy = cov(x_dc,y_dc); cov_zx = cov(z_dc,x_dc); cov_yz = cov(y_dc,z_dc);
  151. X = mean(x_dc); Y = mean(y_dc); Z = mean(z_dc);
  152.  
  153. calc_var_DC = 1.^2 * (2./(X-2*Z+Y)^4 * ((Y-Z).*(Z-X)*cov_xy(3) + ...
  154. (X-Y).*(Y-Z)*cov_zx(3) + (Z-X).*(X-Y)*cov_yz(3) + ...
  155. ((Y - Z)^2/2)*var(x_dc) + ((Z - X)^2/2)*var(y_dc) + ...
  156. ((X - Y)^2/2)*var(z_dc)));
  157.  
  158.  
  159. saveCalcVar_DC = [saveCalcVar_DC calc_var_DC];
  160. save_var_corr(idx) = var(save_lags);
  161.  
  162.  
  163.  
  164.  
  165.  
  166. % LEAST MEAN SQUARES
  167. %
  168. %
  169.  
  170. if referenceTwiceSize % For correlation : Try using a reference window 2x the size.
  171. win_t0 = detrend(pd_trace(1,windowPos-windowSize/2:windowPos+windowSize+windowSize/2),'constant');
  172. else
  173. win_t0 = detrend(pd_trace(1,windowPos:windowPos+windowSize), ... % define a new window t0, force zero-mean (just in case.)
  174. 'constant');
  175. end
  176.  
  177. save_lags_LS = [];
  178.  
  179. for i = 2:timeIncrement:nSegments-1 % define a new window t0, force zero-mean (just in case.)
  180. win = detrend(pd_trace(i,windowPos+coarseDelay:windowPos+coarseDelay+windowSize), ... % define a new window t1, force zero-mean (just in case.)
  181. 'constant');
  182.  
  183. leastSquares = ls_TDE(win_t0,win,windowSize);
  184. if showcorr
  185. figure(2);
  186. subplot(2,1,2), plot(fix(-windowSize/2):1:fix(windowSize/2),leastSquares,'o-'), hold on %title(strcat(num2str(100*i*shift_increment/windowSize),'% of windowSize shift'));
  187. end
  188.  
  189. [~,min_LS]=min(leastSquares); % Find the minimum.
  190.  
  191.  
  192. if min_LS == 1 || min_LS == length(leastSquares) % Just to prevent errors in very low snr
  193. lag = min_LS;
  194. else
  195. lag = quadFit_lag(leastSquares,min_LS); % Parabolic fit using nearest neighbours.
  196. end
  197.  
  198. save_lags_LS = [save_lags_LS lag];
  199.  
  200. if min_LS == 1 || min_LS == length(leastSquares) % This is simply to prevent errors when anomalous detection occurs at very low SNR.
  201. x_ls = [x_ls leastSquares(min_LS)];
  202. y_ls = [y_ls leastSquares(min_LS)];
  203. else
  204. x_ls = [x_ls leastSquares(min_LS+1)];
  205. y_ls = [y_ls leastSquares(min_LS-1)];
  206. end
  207. z_ls = [z_ls leastSquares(min_LS)];
  208.  
  209. end
  210.  
  211. cov_xy = cov(x_ls,y_ls); cov_zx = cov(z_ls,x_ls); cov_yz = cov(y_ls,z_ls);
  212. X = mean(x_ls); Y = mean(y_ls); Z = mean(z_ls);
  213.  
  214. calc_var_LS = 1.^2 * (2./(X-2*Z+Y)^4 * ((Y-Z).*(Z-X)*cov_xy(3) + ...
  215. (X-Y).*(Y-Z)*cov_zx(3) + (Z-X).*(X-Y)*cov_yz(3) + ...
  216. ((Y - Z)^2/2)*var(x_ls) + ((Z - X)^2/2)*var(y_ls) + ...
  217. ((X - Y)^2/2)*var(z_ls)));
  218.  
  219. saveCalcVar_LS = [saveCalcVar_LS calc_var_LS];
  220.  
  221. save_var_LS(idx) = var(save_lags_LS);
  222.  
  223.  
  224. end
  225.  
  226.  
  227. fprintf('\n\n============\n Measured delay : \n \t CORR : %i \n \t LS : %i\n', median(save_lags)-151-225*referenceTwiceSize,median(save_lags_LS)-76)
  228.  
  229.  
  230. % Display Results
  231.  
  232.  
  233. % Theoretical values.
  234. % theocorr = (1/windowSize) * 3/pi.^2 * ( (1 + 2 .* save_snr)./save_snr.^2 + 1);
  235. % theoLS = (1/windowSize) * 3/pi.^2 * ( (2 + 2 .* save_snr)./save_snr.^2);
  236.  
  237. cramer_rao = (1/windowSize) * 3/(pi.^2) * ( (1 + 2 .* save_snr)./save_snr.^2);
  238. save_snr = 10*log10(save_snr);
  239.  
  240. figure(66), semilogy((save_snr),save_var_corr,'.-b',(save_snr),save_var_LS,'.-r',(save_snr),saveCalcVar_DC,'--b',(save_snr),saveCalcVar_LS,'--r',(save_snr),cramer_rao,'-k'),
  241. % hold on
  242. % semilogy((save_snr),theocorr,(save_snr),theoLS);
  243.  
  244.  
  245.  
  246.  
  247.  
  248. % Aux Functions
  249.  
  250. function [ lag ] = quadFit_lag( corr_, coarse )
  251.  
  252. fine = (corr_(coarse+1) - corr_(coarse-1)) / (corr_(coarse+1) + corr_(coarse-1) - 2*corr_(coarse));
  253. lag = -(1/2)*fine + coarse;
  254.  
  255. end
  256.  
  257.  
  258. function [leastSquares] = ls_TDE(win_t0,win, windowSize)
  259. sumResiduals = 0; % Result of the sum of the residuals. Define it as zero for each iteration.
  260. if length(win_t0) == 2*windowSize+1 % IF WE'RE DOING TWICE THE REFERENCE SIZE
  261. for el = 1:windowSize % For position = 0. Just the residual for every pair of points squared.
  262. sumResiduals = sumResiduals + (win(el)-win_t0(el+windowSize/2))^2;
  263. end
  264. else
  265. for el = 1:windowSize % For position = 0. Just the residual for every pair of points squared.
  266. sumResiduals = sumResiduals + (win(el)-win_t0(el))^2;
  267. end
  268. end
  269.  
  270. leastSquares = (1/windowSize * sumResiduals); % Divide it by the number of points.
  271.  
  272. if length(win_t0) == 2*windowSize+1 % IF WE'RE DOING TWICE THE REFERENCE SIZE
  273. for pos = 1:(windowSize/2)
  274. sumResiduals = 0;
  275. for el = 1:windowSize
  276. sumResiduals = sumResiduals + (win(round(el))-win_t0(windowSize/2-pos+el))^2;
  277. end
  278. leastSquares = [(1/(windowSize))*sumResiduals leastSquares];
  279.  
  280. sumResiduals = 0;
  281. for el = 1:windowSize
  282. sumResiduals = sumResiduals + (win(round(el))-win_t0(windowSize/2+pos+el))^2;
  283. end
  284. leastSquares = [leastSquares (1/(windowSize))*sumResiduals];
  285.  
  286. end
  287.  
  288. else
  289. for pos = 1:(windowSize/2) % Shift the position of window t1 up to half the total windowSize.
  290. sumResiduals = 0;
  291. for el = 1:windowSize-pos % Shift to the left, so less one point to calculate.
  292. sumResiduals = sumResiduals + (win(round(pos+el))-win_t0(el))^2; % Add every point from (starting on the amount of points delayed for the wt1, and on zero for wt0).
  293. end
  294.  
  295. leastSquares = [(1/(windowSize-pos))*sumResiduals leastSquares];
  296.  
  297. sumResiduals = 0;
  298. for el = 1:windowSize-pos % Shift to the right.
  299. sumResiduals = sumResiduals + (win(round(el))-win_t0(pos+el))^2; % Add every point from (starting on the amount of points delayed for the wt0, and on zero for wt1).
  300. end
  301. leastSquares = [leastSquares (1/(windowSize-pos))*sumResiduals];
  302. end
  303. end
  304. end
RAW Paste Data