-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbandpass_err_fcn.m
214 lines (180 loc) · 16.8 KB
/
bandpass_err_fcn.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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
function [bandpassed_array,trend_array,nan_ind] = bandpass_err_fcn(input_array,n_dim,delta_dim,low_bound,high_bound,steepness_factor,trend_handling_opt,nan_handling_opt,edge_handling_opt,uneven_edge_handling_opt)
% delta_dim: spacing of grid in dimension of bandpassing
% low_bound: low frequency-wavenumber bound
% high_bound: high frequency-wavenumber bound
% n_dim: number of dimension to bandpass along
% steepness_factor: how rapidly the bandpass filter cuts off at the
% boundaries (must be positive; higher values have a steeper cutoff, and
% possibly more spectral ringing)
% trend_handling_opt: 0 = do not add trend back to bandpassed array, 1 = add trend
% back to bandpassed array
% nan_handling_opt: 0 = do not put NaNs in array, 1 = include NaNs in array
% (trend array and NaN indices are included as outputs regardless)
% edge_handling_opt: 0 = do not replace edge values with NaNs, 1 = replace
% edge values with NaNs
% uneven_edge_handling_opt: 0 = treat edges as even, with zeros for NaNs;
% 1 = adjust padding positions for uneven edges
dim_length = size(input_array,n_dim);
f_vec = (1/(delta_dim*(2*dim_length)))*(mod((0:1:((2*dim_length) - 1))' + dim_length,(2*dim_length)) - dim_length);
% n_avg_ind = round(((1/8)*exp(-mean(log([low_bound high_bound]))))/delta_dim);
n_avg_ind = round(((1/8)*exp(-mean(log([max([low_bound (1/(delta_dim*dim_length))]) min([high_bound (1/(2*delta_dim))])]))))/delta_dim);
% bandpass_filter = 2*0.5*(erf(steepness_factor*(log(abs(f_vec)) - log(low_bound))) - erf(steepness_factor*(log(abs(f_vec)) - log(high_bound))));
bandpass_filter = 0.5*(erf(steepness_factor*(log(abs(f_vec)) - log(low_bound))) - erf(steepness_factor*(log(abs(f_vec)) - log(high_bound))));
% G = [ones(dim_length,1) (1:1:dim_length)'];
% reg_operator_trend = G*(((G')*G)^(-1))*(G');
padding_ind_inside = floor(dim_length/2) + (1:1:dim_length)';
if abs(uneven_edge_handling_opt) < 1e-5
padding_ind_begin = (1:1:floor(dim_length/2))';
padding_ind_end = ((floor(dim_length/2) + dim_length + 1):1:(2*dim_length))';
end
size_array = size(input_array);
prod_size_not_inc_bp_dim = prod(size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]));
input_array_permuted = reshape(permute(input_array,[n_dim (1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]),[dim_length prod_size_not_inc_bp_dim]);
nan_ind_permuted = union(union(find(isnan(input_array_permuted) == 1),find(isinf(input_array_permuted) == 1)),find(abs(input_array_permuted) < (1e-15)*max(abs(input_array_permuted(isinf(input_array_permuted) == 0)))));
clear input_array
input_array_permuted(nan_ind_permuted) = 0;
nan_mask_input_permuted = ones(size(input_array_permuted));
nan_mask_input_permuted(nan_ind_permuted) = 0;
% trend_array_permuted = reshape(reg_operator_trend*(reshape(input_array_permuted,[dim_length prod_size_not_inc_bp_dim])),size_array([n_dim (1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]));
% remove means and trends
index_array_nomean = repmat((1:1:size(input_array_permuted,1))',[1 prod_size_not_inc_bp_dim]) - repmat((sum(nan_mask_input_permuted.*repmat((1:1:size(input_array_permuted,1))',[1 prod_size_not_inc_bp_dim]),1)./sum(nan_mask_input_permuted,1)),[size(input_array_permuted,1) 1]);
input_array_permuted_nomean = input_array_permuted - repmat((sum(nan_mask_input_permuted.*input_array_permuted,1)./sum(nan_mask_input_permuted,1)),[size(input_array_permuted,1) 1]);
trend_array_permuted = (index_array_nomean.*repmat((1./sum((nan_mask_input_permuted.*index_array_nomean).^2,1)).*sum(nan_mask_input_permuted.*index_array_nomean.*input_array_permuted_nomean,1),[size(input_array_permuted,1) 1]));
trend_array_permuted = trend_array_permuted + (input_array_permuted - input_array_permuted_nomean);
input_padded = zeros([(2*dim_length) prod_size_not_inc_bp_dim]);
input_array_permuted_minus_trend = input_array_permuted - trend_array_permuted;
clear input_array_permuted
input_array_permuted_minus_trend(nan_ind_permuted) = 0;
% input_array_variance = sum(nan_mask_input_permuted.*(input_array_permuted_minus_trend.^2),1)./sum(nan_mask_input_permuted,1);
if abs(edge_handling_opt - 1) < 1e-5
test_lags = [1; 2; 3; round(exp((1.2:0.3:log(dim_length))'))];
N = length(test_lags);
autocorr_at_lags = NaN([N prod_size_not_inc_bp_dim]);
integral_timescale = NaN([1 prod_size_not_inc_bp_dim]);
n_lag_step = 1;
complete_flag = 0;
while ((complete_flag == 0) && (n_lag_step <= N))
lag_ind = test_lags(n_lag_step);
autocorr_at_lags(n_lag_step,:) = sum((nan_mask_input_permuted(1:(dim_length - lag_ind),:) & nan_mask_input_permuted((1 + lag_ind):dim_length,:)).*input_array_permuted_minus_trend(1:(dim_length - lag_ind),:).*input_array_permuted_minus_trend((1 + lag_ind):dim_length,:),1)./((sum((nan_mask_input_permuted(1:(dim_length - lag_ind),:) & nan_mask_input_permuted((1 + lag_ind):dim_length,:)).*(input_array_permuted_minus_trend(1:(dim_length - lag_ind),:).^2),1).^(1/2)).*(sum((nan_mask_input_permuted(1:(dim_length - lag_ind),:) & nan_mask_input_permuted((1 + lag_ind):dim_length,:)).*(input_array_permuted_minus_trend((1 + lag_ind):dim_length,:).^2),1).^(1/2)));
if n_lag_step == N
first_neg_ind = find(isnan(integral_timescale) == 1);
else
first_neg_ind = find((autocorr_at_lags(n_lag_step,:) < 0) & (isnan(integral_timescale) == 1));
end
if isempty(first_neg_ind) == 0
first_neg_ind_mask_2D = zeros([1 prod_size_not_inc_bp_dim]);
first_neg_ind_mask_2D(first_neg_ind) = 1;
first_neg_ind_mask = repmat(first_neg_ind_mask_2D,[n_lag_step 1]);
% first_neg_ind_mask(n_lag_step,:) = 0.5*first_neg_ind_mask(n_lag_step,:);
curr_integral_timescale = delta_dim*(first_neg_ind_mask(1,:) + (2*sum(repmat(0.5*(test_lags([(2:1:n_lag_step)'; n_lag_step]) - [0; test_lags(1:(n_lag_step - 1))]),[1 prod_size_not_inc_bp_dim]).*first_neg_ind_mask.*autocorr_at_lags(1:n_lag_step,:),1)));
integral_timescale(first_neg_ind) = curr_integral_timescale(first_neg_ind);
end
if sum(sum(isnan(integral_timescale))) == 0
complete_flag = 1;
end
n_lag_step = n_lag_step + 1;
end
% input_array_permuted_lag1_autocorr = sum((nan_mask_input_permuted(1:(dim_length - 1),:).*input_array_permuted_minus_trend(1:(dim_length - 1),:)).*(nan_mask_input_permuted(2:dim_length,:).*input_array_permuted_minus_trend(2:dim_length,:)),1)./(((sum(nan_mask_input_permuted(1:(dim_length - 1),:).*(input_array_permuted_minus_trend(1:(dim_length - 1),:).^2),1)).^(1/2)).*((sum(nan_mask_input_permuted(2:dim_length,:).*(input_array_permuted_minus_trend(2:dim_length,:).^2),1)).^(1/2)));
% freq_dominant_approx = 1./(4*delta_dim*((1/(1 - input_array_permuted_lag1_autocorr)) - 1));
freq_dominant_approx = 1./(pi*integral_timescale);
end
input_padded(padding_ind_inside,:) = input_array_permuted_minus_trend;
clear autocorr_at_lags input_array_permuted_minus_trend
% pad with erf slopes on either end of the range
if abs(uneven_edge_handling_opt - 1) < 1e-5
first_good_ind = ones([1 prod_size_not_inc_bp_dim]);
last_good_ind = dim_length*ones([1 prod_size_not_inc_bp_dim]);
cum_forward_nan_mask_input_permuted = cumsum(nan_mask_input_permuted,1);
cum_reverse_nan_mask_input_permuted = flipdim(cumsum(flipdim(nan_mask_input_permuted,1)),1);
first_good_ind(abs(cum_forward_nan_mask_input_permuted(dim_length,:)) > 1e-5) = (mod(union(find((abs(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - 1) < 1e-5) & (abs(cum_forward_nan_mask_input_permuted - 1) < 1e-5)),find((abs([((-1)*ones([1 prod_size_not_inc_bp_dim])); cum_forward_nan_mask_input_permuted(1:(dim_length - 1),:)]) < 1e-5) & (abs([((-1)*ones([1 prod_size_not_inc_bp_dim])); cum_forward_nan_mask_input_permuted(2:dim_length,:)] - 1) < 1e-5))) - 1,dim_length) + 1)';
last_good_ind(abs(cum_reverse_nan_mask_input_permuted(1,:)) > 1e-5) = (mod(union(find((abs(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - dim_length) < 1e-5) & (abs(cum_reverse_nan_mask_input_permuted - 1) < 1e-5)),find((abs([cum_reverse_nan_mask_input_permuted(2:dim_length,:); ((-1)*ones([1 prod_size_not_inc_bp_dim]))]) < 1e-5) & (abs([cum_reverse_nan_mask_input_permuted(1:(dim_length - 1),:); ((-1)*ones([1 prod_size_not_inc_bp_dim]))] - 1) < 1e-5))) - 1,dim_length) + 1)';
end
if n_avg_ind > 0
if abs(uneven_edge_handling_opt) < 1e-5
pad_slope_begin = (mean(input_padded(floor(dim_length/2) + ((n_avg_ind + 1):(2*n_avg_ind)),:),1) - mean(input_padded(floor(dim_length/2) + (1:n_avg_ind),:),1))/n_avg_ind;
pad_val_begin = mean(input_padded(floor(dim_length/2) + (1:n_avg_ind),:),1) - (((n_avg_ind - 1)/2)*pad_slope_begin);
pad_slope_end = (mean(input_padded(floor(dim_length/2) + ((dim_length - n_avg_ind + 1):dim_length),:),1) - mean(input_padded(floor(dim_length/2) + ((dim_length - (2*n_avg_ind) + 1):(dim_length - n_avg_ind)),:),1))/n_avg_ind;
pad_val_end = mean(input_padded(floor(dim_length/2) + ((dim_length - n_avg_ind + 1):dim_length),:),1) + (((n_avg_ind - 1)/2)*pad_slope_end);
else
along_dim_ind_array = repmat((1:1:(2*dim_length))',[1 prod_size_not_inc_bp_dim]);
avg_range_mask_1 = zeros(size(input_padded));
avg_range_mask_1((along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 >= 1) & (along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 <= n_avg_ind)) = 1;
avg_range_mask_2 = zeros(size(input_padded));
avg_range_mask_2((along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 >= n_avg_ind + 1) & (along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 <= (2*n_avg_ind))) = 1;
pad_slope_begin = ((sum(avg_range_mask_2.*input_padded,1)./sum(avg_range_mask_2,1)) - (sum(avg_range_mask_1.*input_padded,1)./sum(avg_range_mask_1,1)))/n_avg_ind;
pad_val_begin = (sum(avg_range_mask_1.*input_padded,1)./sum(avg_range_mask_1,1)) - (((n_avg_ind - 1)/2)*pad_slope_begin);
avg_range_mask_1 = zeros(size(input_padded));
avg_range_mask_1((along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 <= -1) & (along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 >= -n_avg_ind)) = 1;
avg_range_mask_2 = zeros(size(input_padded));
avg_range_mask_2((along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 <= -(n_avg_ind + 1)) & (along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 >= -(2*n_avg_ind))) = 1;
pad_slope_end = (-(sum(avg_range_mask_2.*input_padded,1)./sum(avg_range_mask_2,1)) + (sum(avg_range_mask_1.*input_padded,1)./sum(avg_range_mask_1,1)))/n_avg_ind;
pad_val_end = (sum(avg_range_mask_1.*input_padded,1)./sum(avg_range_mask_1,1)) + (((n_avg_ind - 1)/2)*pad_slope_end);
clear avg_range_mask_*
end
else
if abs(uneven_edge_handling_opt) < 1e-5
pad_slope_begin = input_padded(floor(dim_length/2) + 1,:)/0.5;
pad_val_begin = input_padded(floor(dim_length/2) + 1,:);
pad_slope_end = -input_padded(floor(dim_length/2) + dim_length,:)/0.5;
pad_val_end = input_padded(floor(dim_length/2) + dim_length,:);
else
along_dim_ind_array = repmat((1:1:(2*dim_length))',[1 prod_size_not_inc_bp_dim]);
avg_range_mask = zeros(size(input_padded));
avg_range_mask(abs(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1])) < 1e-5) = 1;
pad_slope_begin = (sum(avg_range_mask.*input_padded,1)./sum(avg_range_mask,1))/0.5;
pad_val_begin = sum(avg_range_mask.*input_padded,1)./sum(avg_range_mask,1);
avg_range_mask = zeros(size(input_padded));
avg_range_mask(abs(along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1])) < 1e-5) = 1;
pad_slope_end = (sum(avg_range_mask.*input_padded,1)./sum(avg_range_mask,1))/0.5;
pad_val_end = sum(avg_range_mask.*input_padded,1)./sum(avg_range_mask,1);
clear avg_range_mask
end
end
if abs(uneven_edge_handling_opt) < 1e-5
input_padded(padding_ind_begin,:) = repmat(pad_val_begin,[length(padding_ind_begin) 1]).*(1 + erf(((pi^(1/2))/2)*repmat(abs(pad_slope_begin./pad_val_begin),[length(padding_ind_begin) 1]).*(repmat(reshape(padding_ind_begin,[length(padding_ind_begin) 1]),[1 prod_size_not_inc_bp_dim]) - min(padding_ind_inside))));
input_padded(padding_ind_end,:) = repmat(pad_val_end,[length(padding_ind_end) 1]).*(1 - erf(((pi^(1/2))/2)*repmat(abs(pad_slope_end./pad_val_end),[length(padding_ind_end) 1]).*(repmat(reshape(padding_ind_end,[length(padding_ind_end) 1]),[1 prod_size_not_inc_bp_dim]) - max(padding_ind_inside))));
else
begin_padding_array = repmat(pad_val_begin,[(2*dim_length) 1]).*(1 + erf(((pi^(1/2))/2)*repmat(abs(pad_slope_begin./pad_val_begin),[(2*dim_length) 1]).*(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]))));
input_padded(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) < 0) = begin_padding_array(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) < 0);
end_padding_array = repmat(pad_val_end,[(2*dim_length) 1]).*(1 - erf(((pi^(1/2))/2)*repmat(abs(pad_slope_end./pad_val_end),[(2*dim_length) 1]).*(repmat((1:1:(2*dim_length))',[1 prod_size_not_inc_bp_dim]) - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]))));
input_padded(along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) > 0) = end_padding_array(along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) > 0);
clear along_dim_ind_array
end
input_padded((isnan(input_padded) == 1) | (isinf(input_padded) == 1)) = 0;
% padded_array_variance = sum(input_padded.^2,1)/size(input_padded,1);
bandpassed_array_padded = ifft(repmat(bandpass_filter,[1 prod_size_not_inc_bp_dim]).*fft(input_padded,[],1),[],1);
mean_input_padded = mean(input_padded,1);
clear input_padded
% % make variance adjustment
% bandpassed_array_padded = repmat((input_array_variance./padded_array_variance).^(1/2),[(2*dim_length) 1]).*bandpassed_array_padded;
% add trend back to array, depending on option given
if abs(trend_handling_opt - 0) < 1e-5
bandpassed_array_not_padded = bandpassed_array_padded(padding_ind_inside,:);
elseif abs(trend_handling_opt - 1) < 1e-5
bandpassed_array_not_padded = bandpassed_array_padded(padding_ind_inside,:) + reshape(trend_array_permuted,[dim_length prod_size_not_inc_bp_dim]) + repmat(mean_input_padded,[dim_length 1]);
end
clear bandpassed_array_padded
bandpassed_array_not_padded(nan_ind_permuted) = NaN;
% % pad edges with NaNs where there is low precision about bandpassed values
% n_edge_err = floor(((1/high_bound)/3)/delta_dim);
if abs(edge_handling_opt - 1) < 1e-5
denom_edge_scale = (delta_dim.*high_bound).*((1 - ((low_bound./high_bound).^0.5)).^1.25);
n_edge_err = round(max([min([((((0.9*(low_bound./((pi/2)*freq_dominant_approx))) - (3/40))./(delta_dim.*((((pi/2)*freq_dominant_approx).*low_bound).^(1/2)))) + ((0.45 - (1.8*(low_bound./((pi/2)*freq_dominant_approx))))./denom_edge_scale)); (0.15./(delta_dim.*((((pi/2)*freq_dominant_approx).*low_bound).^(1/2))))],[],1); (0.3./repmat(denom_edge_scale,[1 size(freq_dominant_approx,2)]))],[],1));
% bandpassed_array_not_padded(1:n_edge_err,:) = NaN;
% bandpassed_array_not_padded((dim_length - n_edge_err + 1):dim_length,:) = NaN;
if abs(uneven_edge_handling_opt) < 1e-5
bandpassed_array_not_padded(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - repmat(1 + n_edge_err,[dim_length 1]) < 0) = NaN;
bandpassed_array_not_padded(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - repmat(dim_length - n_edge_err,[dim_length 1]) > 0) = NaN;
else
bandpassed_array_not_padded(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - repmat(first_good_ind + n_edge_err,[dim_length 1]) < 0) = NaN;
bandpassed_array_not_padded(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - repmat(last_good_ind - n_edge_err,[dim_length 1]) > 0) = NaN;
end
end
bandpassed_array = real(permute(reshape(bandpassed_array_not_padded,[dim_length size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]));
trend_array = permute(reshape(trend_array_permuted,[dim_length size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
clear bandpassed_array_not_padded trend_array_permuted
nan_ind = find(isnan(bandpassed_array) == 1);
% remove NaNs from array, depending on option given
if abs(nan_handling_opt - 0) < 1e-5
bandpassed_array(nan_ind) = 0;
end