-
Notifications
You must be signed in to change notification settings - Fork 4
/
sdmse.m
150 lines (130 loc) · 4.31 KB
/
sdmse.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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
function [eta,mseO,mseR,mseC,fO,fR,fC]=sdmse(SN,eta,TH,L,phi)
% [eta,mseO,mseR,mseC,fO,fR,fC]=sdmse(SN,eta,TH,L,phi)
%
% Calculates the colatitudinal mean square error AVERAGED over the
% (complementary) region of observation or over the entire sphere
% for spherical harmonic expansions, by Gauss-Legendre averaging
%
% INPUT:
%
% SN Signal-to-noise ratio, scalar
% eta Damping parameter (0 <= eta <= 1) or
% A single number identifiying a linearly spaced set
% TH Angular extent of the spherical caps, in degrees
% L Bandwidth (maximum angular degree)
% phi Perhaps you want to put in a particular longitude (default: not)
%
% OUTPUT:
%
% eta Damping parameter
% mseO Mean square error averaged over the sphere
% mseR Mean square error averaged over the region
% mseC Mean square error averaged over the complementary region
% fO Eta of minimum mean square error over the sphere
% fR Eta of minimum mean square error over the region
% fC Eta of minimum mean square error over the complementary region
%
% Last modified by fjsimons-at-alum.mit.edu, 09/23/2023
%
% See also SDETA, SDMSK, SDMSE2
defval('SN',10)
defval('eta',100)
defval('TH',10)
defval('L',45)
defval('phi',NaN);
% Signal and noise variances
S=1;
N=S/SN;
% (Number of) Damping or truncation levels tried
if eta>1 & length(eta)==1
eta=linspace(0,1,eta);
end
% For the double polar cap
sord=2;
% Start doing colatitudes only
% Gauss-Legendre points ! It's a sum of squared polynomials...
intvO=[-1 1];
[wO,xO]=gausslegendrecof(2*L,[],[-1 1]);
intvR=cos([180-TH TH]*pi/180);
[wR,xR]=gausslegendrecof(2*L,[],intvR);
intvC=[cos(TH*pi/180) 1];
[wC,xC]=gausslegendrecof(2*L,[],intvC);
% Specify integration points
thetaO=acos(xO);
thetaR=acos(xR);
thetaC=acos(xC);
% The area of measurement is the belt
srt='belt';
% Load precomputed Slepian functions?
fnpl=sprintf('%s/SDERRGL-%i-%i-%i.mat',...
fullfile(getenv('IFILES'),'SDERR'),TH,L,phi);
if exist(fnpl,'file')==2
eval(sprintf('load %s',fnpl))
disp(sprintf('load %s',fnpl))
else
% Calculate eigenfunctions and eigenvalues
[GO,Vo,jk1,jk2,jk3,KA,K]=galpha(TH,L,sord,thetaO,phi,srt);
% Calculate eigenfunctions and eigenvalues
GR=galpha(TH,L,sord,thetaR,phi,srt);
% Calculate eigenfunctions and eigenvalues
GC=galpha(TH,L,sord,thetaC,phi,srt);
eval(sprintf(['save %s '...
'GO GR GC Vo thetaO thetaR thetaC KA K'],...
fnpl))
end
% Load precomputed variance results?
fnpl=sprintf('%s/SDERRGL-%i-%i-%i-%i-%i.mat',...
fullfile(getenv('IFILES'),'SDERR'),SN,length(eta),TH,L,phi);
if exist(fnpl,'file')==2 & eta(1)==0 & eta(end)==1
eval(sprintf('load %s',fnpl))
disp(sprintf('load %s',fnpl))
else
% DAMPED SPHERICAL HARMONIC INVERSIONS
mseO=repmat(NaN,1,length(eta));
mseR=repmat(NaN,1,length(eta));
mseC=repmat(NaN,1,length(eta));
for index=1:length(eta)
% Variance matrix
V=sdvar(N,eta(index),Vo);
% Bias matrix
B=sdbias(S,eta(index),Vo);
% Error matrix for the damped spherical harmonic case
% Diagonal terms only, at identical locations
errSHO=diag(GO'*(V+eta(index)^2*B)*GO);
errSHR=diag(GR'*(V+eta(index)^2*B)*GR);
errSHC=diag(GC'*(V+eta(index)^2*B)*GC);
% Mean square error over the entire sphere
% Integral over the entire sphere, divided by area
mseO(index)=wO(:)'*errSHO/range(intvO);
% Integral over the belt
mseR(index)=wR(:)'*errSHR/range(intvR);
% Integral over both caps (factor two cancels)
mseC(index)=wC(:)'*errSHC/range(intvC);
end
% Check that these things add up correctly
difer(mseO-mseR*range(intvR)/2-2*mseC*range(intvC)/2,[],0);
% So much for the calculation, now for the prediction
for und=1:length(eta)
da=Vo+eta(und)*(1-Vo);
fO(und)=eta(und)-sum(Vo.*(1-Vo)./da.^3)...
./sum(Vo.*(1-Vo).^2./da.^3)/SN;
fR(und)=eta(und)-sum(Vo.^2.*(1-Vo)./da.^3)...
./sum(Vo.^2.*(1-Vo).^2./da.^3)/SN;
fC(und)=eta(und)-sum(Vo.*(1-Vo).^2./da.^3)...
./sum(Vo.*(1-Vo).^3./da.^3)/SN;
end
% Find minimum variance and compare to theory
[a,fO]=min(abs(fO));
difer(mseO(fO)-min(mseO))
fO=eta(fO);
[c,fR]=min(abs(fR));
difer(mseR(fR)-min(mseR))
fR=eta(fR);
[e,fC]=min(abs(fC));
difer(mseC(fC)-min(mseC))
fC=eta(fC);
if eta(1)==0 & eta(end)==1
eval(sprintf([...
'save %s eta mseO mseR mseC fO fR fC'],fnpl))
end
end