-
Notifications
You must be signed in to change notification settings - Fork 0
/
periodicUncertainty.m
80 lines (69 loc) · 11.5 KB
/
periodicUncertainty.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
function sigma = periodicUncertainty(eta,den)
% Estimates the uncertainty in the timing vbetween peaks
% for a discretely sampled sinusoidal signal with gaussian noise
% y(t) = sin(omega*t) + Norm(0,eta)
% where omega = 2*pi / T
% and t=0:dt:Inf (discretely sampled with interval dt)
%
% ====== Input Parameters ==================
%
% eta = noise level = Amplitude / std(noise)
% den = T / dt
%
% ==========================================
%
% Motivation/Implementation:
%
% The uncertainty of was estimated by computation for 1e5 oscillations
% at various combinations of (eta,den) and it is estimated by linear
% interpolation for points off of this precomputed grid
%
% It was found that the distrubution of inter-peak timings is gaussian,
% and centered on the true period used to generate the signals.
% This function provides a good estimate of the standard deviation of
% the gaussian distributions, which depends in a non-linear fashion on
% both of the input parameters.\
eta(eta<=0)=min(eta(eta>0));
global den0 eta0 sigma0
if isempty(den0)
den0=[4 8 16 32 64 128 256 512 724 1024];
eta0=[0 1e-05 0.0204179591836735 0.0408259183673469 0.0612338775510204 0.0816418367346939 0.102049795918367 0.122457755102041 0.142865714285714 0.163273673469388 0.183681632653061 0.204089591836735 0.224497551020408 0.244905510204082 0.265313469387755 0.285721428571429 0.306129387755102 0.326537346938776 0.346945306122449 0.367353265306122 0.387761224489796 0.408169183673469 0.428577142857143 0.448985102040816 0.46939306122449 0.489801020408163 0.510208979591837 0.53061693877551 0.551024897959184 0.571432857142857 0.591840816326531 0.612248775510204 0.632656734693878 0.653064693877551 0.673472653061225 0.693880612244898 0.714288571428571 0.734696530612245 0.755104489795918 0.775512448979592 0.795920408163265 0.816328367346939 0.836736326530612 0.857144285714286 0.877552244897959 0.897960204081633 0.918368163265306 0.93877612244898 0.959184081632653 0.979592040816327 1];
sigma0=[0 0 0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0 0.000167216733187021 0.000538594159026628;0 0 0.00795796493399823 0.0275163460447366 0.0251117357635168 0.0233413930532968 0.0219940642652746 0.020874630208135 0.0203921857785128 0.0199105586222471;0 0 0.0352605279470888 0.0367084699153093 0.0341803120852728 0.0322733622369412 0.0301672168777187 0.0289889636498022 0.0284045150658484 0.02749885708232;0 0.00463301880866539 0.0482529331222217 0.0442160422976132 0.0411726436258915 0.038719007208188 0.0365604817962252 0.0349373878615167 0.0341725460917163 0.0335936833172278;0 0.0194467023041105 0.0552962704172323 0.0504069801098633 0.0469422091600725 0.0442079753700823 0.0420226669460289 0.04013925119215 0.0391062928216625 0.038456862669362;0 0.0378486391806148 0.0607922717735825 0.0557540888382545 0.0520368895368214 0.0490396357567213 0.0464364441606347 0.044540052878321 0.0436909377674061 0.0430017358772549;0 0.0603984539284788 0.0660106427771753 0.0605774128747216 0.0566135836477046 0.0535356694029325 0.0510634091809326 0.0488050282966399 0.0477222399757409 0.0472531625332763;0 0.0815685007760814 0.0721087012801401 0.0652566804911105 0.0609461372212967 0.0574585347246221 0.0546175555379862 0.0529446124902861 0.0520913505779361 0.0511565081249077;0.00449493817804678 0.101218766788684 0.0780298654206908 0.0697795533532848 0.0648829519371232 0.0616626175027078 0.0587470928486827 0.0567760954510013 0.0556330479665272 0.0551782414103948;0.0114604675473153 0.122382709339387 0.0853035829686103 0.0747058751566651 0.0691351336413619 0.065145558603258 0.0626179531578431 0.0602906756993893 0.0592325393791448 0.0586851868029179;0.0233453840996579 0.13504212465591 0.093883849957335 0.0794988788037953 0.0724802115020946 0.0683922683269488 0.0655979126491223 0.0634949596746868 0.061944385765617 0.0616619856884449;0.039924439220234 0.151245753159067 0.101251860283555 0.0834358526552035 0.0755684141666431 0.0718257476120522 0.0687359022863134 0.0667560907072045 0.0653769641388523 0.0644008829907055;0.0597740736798335 0.1630969768885 0.112195739978361 0.0878590491256766 0.0790696259695254 0.0750446039451841 0.0715346815049466 0.0692796248948373 0.0682743285555554 0.06722418622887;0.0845389912246113 0.174583338415192 0.115979068078214 0.0921058853995532 0.0822345381518018 0.0777296615537455 0.0746652444077066 0.0719774275925962 0.0709716348301568 0.0693085287255602;0.107033520837682 0.185006336149924 0.122694459769382 0.0942762619849246 0.0849632545959282 0.0806956684253898 0.0769664523900384 0.0745621157064068 0.0735233180735542 0.0718684767045576;0.135591599506387 0.195932585521976 0.128066661333644 0.0964328628104407 0.0874775667779204 0.0830961758850663 0.0799010036988856 0.0769584006683248 0.075418218717068 0.0742191189256464;0.160821890410338 0.205064374761007 0.136102391933557 0.0995730556396398 0.0897907637688793 0.0860242246692556 0.0828760698754304 0.079346848227837 0.078412897776597 0.0772158974771005;0.187232348792869 0.2146358824571 0.138293656332542 0.100805940016515 0.0923155219472297 0.0883092132526536 0.0847952566876096 0.0815824723521121 0.0802791671863134 0.0790946388176888;0.209187052860506 0.218840746991764 0.142725699000972 0.10303163262593 0.0947981389743969 0.0910344336992703 0.0872890554514067 0.0838569179943918 0.0828112147007577 0.0812106776681071;0.234581938379756 0.230704069277405 0.148154721264338 0.10510792599248 0.0963957075991864 0.0923705329108983 0.0891302274191235 0.0860354289346616 0.0847562514247498 0.0832822231996734;0.258082306221332 0.238023032266927 0.149822464708314 0.107934254884741 0.0988181227641825 0.0946270323906431 0.091396911312797 0.0884342317437982 0.0866633689411956 0.0850908660145102;0.275814071499036 0.245482388585262 0.15400533111076 0.108733230436346 0.100878180072452 0.0966675988766992 0.0928798948222617 0.0902781999107741 0.0890091283370911 0.0871326824426818;0.294175827170229 0.252461145747259 0.152123158091731 0.111074979482578 0.102788344199196 0.0987041028837534 0.0951470674789337 0.0916647075237698 0.0901985720559396 0.0891831126282349;0.313483404198309 0.254153794190094 0.155863292562463 0.112637801382163 0.104831576677142 0.100618581299727 0.0969426496592525 0.0942796223875339 0.0921839696657216 0.0910207727991452;0.325540681695771 0.259662339174044 0.158545903342636 0.114136712667267 0.106352704247069 0.102374497933186 0.0992241896933977 0.0957183961408276 0.0941228327091231 0.0929540838763703;0.34070563892979 0.262049637638878 0.159069667269538 0.115142560813691 0.107645363665396 0.104363553068262 0.101022537500487 0.0978949050011616 0.0957674466526638 0.0942505844910534;0.353821569743488 0.268097724500007 0.161304577669176 0.116799416980614 0.110220068603446 0.106040123218354 0.102643504480658 0.0991074951871595 0.097762506021243 0.0965788044419196;0.36734548900081 0.273750213611813 0.163960431302483 0.11855801195844 0.111502885681747 0.10741872381817 0.104103344761781 0.100673257057166 0.0994421124443855 0.0978608899635706;0.374810559099222 0.274501970842274 0.164270186158587 0.120015882175679 0.112926902366568 0.10902993387388 0.105673079490261 0.102404157905211 0.100686637528415 0.099559484018666;0.384280458168218 0.27743282935198 0.166614482577709 0.121209270927315 0.114314275031918 0.110587525141772 0.10698840093021 0.103983189728397 0.102627249422125 0.100840169909179;0.390977946558293 0.281493057343607 0.164814169623821 0.122596535579722 0.116036611622507 0.112254405054321 0.108835189744902 0.105494650754082 0.10379748559834 0.102256735458359;0.399224423274957 0.28119201537678 0.167097266861813 0.124200709884631 0.117194439123809 0.113664245828873 0.110006269339114 0.107289281287306 0.105244265294734 0.10399142548791;0.406723558599659 0.284744547576882 0.167343667061751 0.125807742209499 0.119245473486518 0.115318409565918 0.111277079207954 0.108296655242168 0.10663842321568 0.10490293274276;0.411622164994793 0.287459251491073 0.169263185330808 0.126304313347819 0.120084179668847 0.116413863468462 0.113033757380587 0.109525957470476 0.107878168602047 0.106234326354155;0.417682127280332 0.28744643522828 0.172468806367494 0.127973509922196 0.121463118753201 0.117950781261437 0.114592253767719 0.111174626505401 0.109202892669052 0.107877946586128;0.422101184327038 0.291479686250745 0.172369812526624 0.12955932083003 0.123285327169418 0.118663118354147 0.11583661667287 0.112491129658531 0.111121903402865 0.108947200660061;0.424430005429608 0.289163862537616 0.174251720678334 0.131288640340588 0.124204736631231 0.120567000001841 0.116395608079371 0.113499241134999 0.112234527364015 0.110399534818735;0.428955318732811 0.290966086342794 0.174174121205226 0.132553430271991 0.125528813770043 0.121750185085015 0.118166437427975 0.115134058010545 0.113468260923711 0.1115565929181;0.432483691306632 0.29310603946419 0.175444774768931 0.134314116078118 0.127234958502551 0.123156785621216 0.119185635281957 0.115966701059858 0.11444922406307 0.112783835074126;0.433814744647033 0.29453775049305 0.17800921665173 0.136002380315529 0.128543415141938 0.124581283055288 0.12048011489362 0.117327246102255 0.115586810437847 0.114263010186321;0.434063812646694 0.296603835611339 0.178180971139982 0.136899366949859 0.129516231046485 0.125872497351834 0.122010845361782 0.118718467577744 0.11701724216938 0.115224106707609;0.44064103252091 0.296405315159288 0.178089977393662 0.138485986493343 0.130955244415031 0.126581294139553 0.123436300443306 0.119232351491648 0.118337303063296 0.116268350800995;0.439897412647765 0.297320880648902 0.180331814604975 0.139465921331888 0.132385034970545 0.127871381953999 0.124146589986805 0.120373635936551 0.118686397741644 0.117454643480668;0.441820074401881 0.298149125287685 0.181305228297702 0.141321425843668 0.133014135810169 0.128798510731085 0.125594116137686 0.12177904458047 0.119653902128978 0.118747247841502;0.440359440173477 0.298036512572929 0.183747728385398 0.14197623721344 0.135200660728465 0.130454172455509 0.1265050669098 0.122864394867843 0.121446974917723 0.119909387353336;0.445288649608908 0.297294984326996 0.184003351100748 0.144382145554597 0.136497110134418 0.131562698414185 0.127366906793172 0.123776510953841 0.122265006515625 0.12053974064904;0.442676363392054 0.297672685328428 0.184342738954514 0.145248457465936 0.137064589664427 0.132248188632748 0.128464612310818 0.125151510162406 0.123565708877633 0.121859554852159;0.443553574187659 0.298524432976607 0.18601027726785 0.146078423871254 0.138076019857517 0.133848606376707 0.129877494149706 0.126204830895979 0.124642454989757 0.123548443872008;0.444164043774572 0.298687857267062 0.187455425113084 0.148561597214282 0.139409036942638 0.134800362615567 0.130425437270955 0.127526298919826 0.125884071690797 0.124549947910363;0.447086739534848 0.298161115199698 0.187477704350284 0.15024904717843 0.14057051201911 0.135378402059176 0.131384755126234 0.128341175765186 0.126843836677155 0.125049959106927];
end
function sig = calcUncert(n,d)
d=min(den0(end),d);
n=min(eta0(end),n);
denA=find(den0<=d,1,'last');
denB=find(den0>=d,1);
if denA==denB
alpha=0;
else
alpha=(d-den0(denA))./abs(den0(denB)-den0(denA));
end
etaA=find(eta0<=n,1,'last');
etaB=find(eta0>=n,1);
if etaA==etaB
beta=0;
else
beta=(n-eta0(etaA))./abs(eta0(etaB)-eta0(etaA));
end
try
sig=(((1-alpha)*sigma0(etaA,denA)+alpha*sigma0(etaA,denB))*(1-beta))+(beta*(((1-alpha)*sigma0(etaB,denA)+alpha*sigma0(etaB,denB))));
catch e
disp(e)
end
end
if length(den)~=1
N=length(den(:));
sigma=zeros(size(den));
for i=1:N
try
sigma(i)=calcUncert(eta(i),den(i));
catch e
getReport(e)
end
end
else
sigma=calcUncert(eta,den);
end
end