This is a static copy of a profile reportHome
calc_loadcorrect_final (1 call, 2.219 sec)
Generated 05-Aug-2011 13:00:59 using cpu time.
function in file /home/LeechJ/cbass_analysis/reduc/calc_loadcorrect_final.m
Copy to new window for comparing multiple runs
Parents (calling functions)
Lines where the most time was spent
| Line Number | Code | Calls | Total Time | % Time | Time Plot |
| 144 | load_shift = circshift(load_re... | 1000 | 1.563 s | 70.4% |  |
| 158 | [x,f_val] = fminsearch(test_fu... | 10 | 0.328 s | 14.8% | |
| 26 | dcut = cutObs(d, 'blank', 'onl... | 1 | 0.197 s | 8.9% | |
| 145 | value(shiftIndex+1) = (dot((sk... | 1000 | 0.044 s | 2.0% | |
| 165 | load_res = circshift(load_res,... | 10 | 0.033 s | 1.5% | |
| All other lines | | | 0.055 s | 2.5% | |
| Totals | | | 2.219 s | 100% | |
Children (called functions)
Code Analyzer results
| Line number | Message |
| 35 | FOR might not be aligned with its matching END (line 176). |
| 110 | Using ISEMPTY is usually faster than comparing LENGTH to 0. |
| 114 | Using ISEMPTY is usually faster than comparing LENGTH to 0. |
| 145 | The variable 'value' appears to change size on every loop iteration. Consider preallocating for speed. |
| 149 | The value assigned here to 'H' appears to be unused. Consider replacing it by ~. |
| 152 | The variable 'shift' appears to change size on every loop iteration. Consider preallocating for speed. |
| 158 | The value assigned to variable 'f_val' might be unused. |
| 160 | The variable 'K_best' appears to change size on every loop iteration. Consider preallocating for speed. |
| 173 | The variable 'K_best' appears to change size on every loop iteration. Consider preallocating for speed. |
| 174 | The variable 'shift' appears to change size on every loop iteration. Consider preallocating for speed. |
Coverage results
[ Show coverage for parent directory ]
| Total lines in function | 201 |
| Non-code lines (comments, blank lines) | 109 |
| Code lines (lines that can run) | 92 |
| Code lines that did run | 56 |
| Code lines that did not run | 36 |
| Coverage (did run/can run) | 60.87 % |
Function listing
time calls line
1 function [K_best,shift,dorig] = calc_loadcorrect_final(d)
2
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 %
5 % function [K_best,shift,dorig] = calc_loadcorrect_final(d)
6 %
7 % To calculate the correction factor to remove the effect
8 % of fluctuations in the cold load
9 %
10 % Do this on each of the 6 channels after alpha corrs and before r-factor corrections
11 % Returns arrays of best fit amplitude and shift (offset of thermal data) for each channel
12 %
13 % Use this function on a nice bit of static,clean data --> the factors should then be
14 % useable on any other data.
15 %
16 % Should check this works over various timescales/ conditions before
17 % releasing as final
18 %
19 % act
20 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21
22 % This will work on ony the selected 'static' dataT (r-factored and
23 % corrected)
24
25 % Do this on blank, static data
0.20 1 26 dcut = cutObs(d, 'blank', 'only');
27
28 % change the flags so that they match the data
0.01 1 29 flags = convertFlagToDataSize(dcut.flags, ...
30 size(dcut.antenna0.receiver.data,2));
31
1 32 numDims = size(dcut.antenna0.receiver.data,2);
33
34
1 35 for chan=1:numDims
36
10 37 if(isfield(dcut.antenna0.receiver, 'dataT'))
38 isT = 1;
39 temp = [dcut.antenna0.receiver.dataT(:,chan), dcut.antenna0.thermal.coldLoad(:)];
10 40 else
10 41 isT = 0;
10 42 temp = [dcut.antenna0.receiver.data(:,chan), dcut.antenna0.thermal.coldLoad(:)];
10 43 end
44
45 % apply flags
10 46 ind = flags.fast(:,chan);
10 47 temp(ind,:) = nan;
48
49 % Now we should remove the mean of the load temperature
0.01 10 50 load_res = temp(:,2)-nanmean(temp(:,2));
51
52 % Already r-factored so don't need to remove mean from the main data set
10 53 sky_res = temp(:,1) - nanmean(temp(:,1));
54
55 % we don't want the nan values to contribute, so set them to zero.
10 56 load_res(isnan(load_res)) = 0;
10 57 sky_res(isnan(sky_res)) = 0;
58
59 % At this point have de-trended , mean removed data in sky_res and load_res
60 % all flagged data masked (not in these arrays)
61
62 % Now let's try just finding the factor by which we need to multiply the
63 % load to subtract it from the sky signal
64 % C = sum(sky.load)
65 % Minimise via factor K such that C_min = sum((sky-Kload).load)
66
67 % First we need to shift load_res around in time to match the sky_res phase
68 % - thermal lag??
69
70 % want over a full period that is not flagged
71
72 % 17-Feb-2011 (MAS): This has the potential to be really, really slow.
73 % for mm=1:length(ind)-100
74 % if(ind(mm)==1)
75 % consecutiveLength(mm) = 0;
76 % else
77 % ff = find(ind(mm:mm+100)==1,1);
78 % if(isempty(ff))
79 % consecutiveLength(mm) = 100;
80 % else
81 % consecutiveLength(mm) = ff;
82 % end
83 % end
84 % end
85 % fstart = find(consecutiveLength == max(consecutiveLength),1);
86 % fend = fstart + max(consecutiveLength)-1;
87
88
89 % % This is a lot faster.
10 90 diffInd = cat(1,0,diff(ind));
10 91 startInd = find(diffInd < 0);
10 92 stopInd = find(diffInd > 0);
93
10 94 length(startInd);
10 95 length(stopInd);
96
97 % check they're the same length.
10 98 if(length(startInd)~=length(stopInd))
99 if(length(startInd)<length(stopInd))
100 if(startInd(1)>stopInd(1))
101 stopInd(1) = [];
102 end
103 else
104 if(startInd(2)<stopInd(1))
105 startInd(1) = [];
106 end
107 end
108 end
109
10 110 if (length(startInd)>0)
111
112 % First - if all the data is ok i.e no unflagged data
113
10 114 if (length(stopInd)==0)
115 fstart = startInd(1);
116 fend = fstart + 99;
10 117 else
118 % This bit assumes that there is both flagged and unflagged data
10 119 if startInd(1) > stopInd(1)
120 % keyboard
10 121 startInd = cat(1,1,startInd);
10 122 stopInd = cat(1,stopInd,length(ind));
10 123 firstLong = find((stopInd - startInd) >= 100,1);
10 124 fstart = startInd(firstLong);
10 125 fend = fstart + 99;
126 else
127 if startInd(end) > stopInd(end)
128 stopInd = cat(1,startInd(1)-1,stopInd,length(ind));
129 startInd = cat(1,1,startInd);
130 firstLong = find((stopInd - startInd) >= 100,1);
131 fstart = startInd(firstLong);
132 fend = fstart + 99;
133 else
134 firstLong = find((stopInd - startInd) >= 100,1);
135 fstart = startInd(firstLong);
136 fend = fstart+99;
137 end
138 end
10 139 end
140
141
10 142 shiftIndex = 0;
10 143 for mm=fstart:fend
1.56 1000 144 load_shift = circshift(load_res, shiftIndex);
0.04 1000 145 value(shiftIndex+1) = (dot((sky_res - load_shift), load_shift));
1000 146 shiftIndex = shiftIndex+1;
1000 147 end
148
10 149 [H,I]=max(value);
150
10 151 required_shift = I-1;
10 152 shift(chan) = (required_shift);
153
0.01 10 154 load_res = circshift(load_res,required_shift);
155
10 156 test_function = @(x) (std((sky_res - x*load_res)));
157
0.33 10 158 [x,f_val] = fminsearch(test_function, 1);
159
10 160 K_best(chan) = x;
161
162
163 %Put thermal load data back to how it was originally
164
0.03 10 165 load_res = circshift(load_res,-required_shift);
166
167 else
168 % And what if all the data is flagged
169 if (chan ==1)
170 disp('All your blank data is flagged') % put warning in once on first channel
171 disp('Will not apply a correction')
172 end
173 K_best(chan) = 0; %% Once we know what the stability of this
174 shift(chan) = 0; %% correction is we can use logged values
175 end
10 176 end
177
178
179
180 % Apply correction to dcut data and look at power spectrum to check it
181 % worked
1 182 for chan = 1:numDims
0.01 10 183 shifted_load = circshift(load_res,shift(chan));
0.01 10 184 correction = K_best(chan)*(shifted_load - mean(shifted_load));
185
10 186 if(isT)
187 dcut.antenna0.receiver.dataT(:,chan) = ...
188 (dcut.antenna0.receiver.dataT(:,chan))- correction;
10 189 else
10 190 dcut.antenna0.receiver.data(:,chan) = ...
191 (dcut.antenna0.receiver.data(:,chan))- correction;
10 192 end
10 193 end
194
195
196 % For checking - keep an array which is the original data for comparison
197
1 198 dorig = d;
199
200
1 201 return;
Other subfunctions in this file are not included in this listing.