-
Notifications
You must be signed in to change notification settings - Fork 0
/
lib_filters_obp.py
256 lines (221 loc) · 10.3 KB
/
lib_filters_obp.py
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
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
import numpy as np
import scipy.signal as sg
#import scipy.signal as sgw # this is for older versions of scipy
import scipy.signal.windows as sgw
import matplotlib.pyplot as plt
import scipy.constants as consts
def G_field_oneway(theta):
# This antenna pattern function was fitted on measured antenna patterns by Alejandro Bohe (hence the "strange" 2.25 power)
# this is for the field (hence the sqrt, because the fit was on the "power") . This was corrected on 2024/05/26
return np.sqrt(np.exp(-np.abs(theta/(np.radians(.107)/2))**2.25*np.log(2)))
def get_obp_filter(L_filt = 1, sampling_in = 0.025, f_axis = None, plot_flag = True, kernel="parzen"):
"""
Get the kernel shape and the spectral response of the filter type specified in the arguments
Arguments:
L_filt, Kernel length in km
sampling_in, Spatial spacing (in km) of the filter kernel
f_axis, Frequency axis where to compute the spectral response of the filter
plot_flag, Plot filter kernel and spectral response
kernel, Kernel to compute ('parzen' or 'bharris' -for Blackman-Harris)
"""
Nparzen = int(np.round(L_filt/sampling_in))
if Nparzen%2 == 0:
Nparzen += 1
#print("Nb of points OBP kernel: %d" % Nparzen)
if kernel == "bharris":
w_obp = sgw.blackmanharris(Nparzen) # parzen of 41 points if input sampling is 25 m (kernel length is ~1 km)
if kernel == "parzen":
w_obp = sgw.parzen(Nparzen) # parzen of 41 points if input sampling is 25 m (kernel length is ~1 km)
if kernel == "sinc2":
xn=(np.arange(0,Nparzen-1,1)-Nparzen//2)*2*np.pi
w_obp = (np.sinc(xn*L_filt*0.00222))**2
if kernel == "alejandro_azptr_beam0":
PRF = 4332
N_burst = 9
v_sat = 7310.
H = 875.0e3
f_s = 300.0e6
v_ground = v_sat/(1+H/6377.0e3)
lambda_c = (consts.c/35.75e9)
incidence_deg = 3.5 # just used for the ground/range conversion for xtrack
d = N_burst*v_ground/PRF # along track posting before OBP averaging
w_ptr = np.zeros(512)
w_ptr[:401] = np.sinc(2*v_sat*N_burst/(lambda_c*H*PRF)*np.arange(-200,201)*d)**2
w_ptr /= np.sum(w_ptr)*d
w_obp = np.zeros(512)
y_tab = np.arange(-200,201)*d+1.0e-3
ant = G_field_oneway(y_tab/H)**4
w_obp[:len(y_tab)] = np.sin(2*np.pi*v_sat*N_burst/(lambda_c*H*PRF)*y_tab)**2/np.sin(2*np.pi*v_sat/(lambda_c*H*PRF)*y_tab)**2*ant
w_obp /= np.sum(w_obp)*d
sampling_in=d/1000.
if kernel == "alejandro_azptr":
PRF = 4332
N_burst = 9
v_sat = 7310.
H = 875.0e3
f_s = 300.0e6
v_ground = v_sat/(1+H/6377.0e3)
lambda_c = (consts.c/35.75e9)
incidence_deg = 3.5 # just used for the ground/range conversion for xtrack
d = N_burst*v_ground/PRF # along track posting before OBP averaging
w_ptr = np.zeros(512)
w_ptr[:401] = np.sinc(2*v_sat*N_burst/(lambda_c*H*PRF)*np.arange(-200,201)*d)**2
w_ptr /= np.sum(w_ptr)*d
w_obp = np.zeros(512)
y_tab = np.arange(-200,201)*d+1.0e-3
ant = G_field_oneway(y_tab/H)**4
# w_obp[:len(y_tab)] = np.sin(2*np.pi*v_sat*N_burst/(lambda_c*H*PRF)*y_tab)**2/np.sin(2*np.pi*v_sat/(lambda_c*H*PRF)*y_tab)**2*ant
# Defines the PTR for each of the 9 beams
w_ale9 = np.zeros(512)
w_ptrantana_b = np.zeros((9,512))
w_ptrantana_c = np.zeros((9,512))
b_vals = np.arange(-4,5)
arg_b = 0.8/N_burst*b_vals[:,None]+2*v_sat/(lambda_c*H*PRF)*y_tab[None,:]
w_ptrantana_b[:,:len(y_tab)] = np.sin(np.pi*N_burst*arg_b)**2/np.sin(np.pi*arg_b)**2*(G_field_oneway(y_tab/H)**4)[None,:]
w_ptrantana_b = w_ptrantana_b/(np.sum(w_ptrantana_b, axis=1)[:,None]*d)
weights = np.array([0.07019711, 0.09767803, 0.12220272, 0.13814296, 0.14355835,0.13814296, 0.12220272, 0.09767803, 0.07019711])
for ib in range(9) :
w_ale9=w_ale9+weights[ib]*np.roll(w_ptrantana_b[ib,:],14*(ib-4))
w_ptrantana_c[ib,:] = np.roll(w_ptrantana_b[ib,:],14*(ib-4))
# these weights actually vary along the track. This is one set of weights obtained by averaging along one (?) cycle
w_obp[:len(y_tab)]=w_ale9[:len(y_tab)]
w_obp /= np.sum(w_obp)*d
sampling_in=d/1000.
w_obp /= np.sum(w_obp)
x_axis = np.arange(-(len(w_obp)-1)/2, (len(w_obp)-1)/2 + 1)*sampling_in
# OBP spectrum
if type(f_axis) != np.ndarray:
if f_axis == None:
f_obp, S_obp = sg.freqz(w_obp, fs=1/sampling_in)
else:
f_obp, S_obp = sg.freqz(w_obp, fs=1/sampling_in, worN=f_axis)
else:
f_obp, S_obp = sg.freqz(w_obp, fs=1/sampling_in, worN=f_axis)
S_obp = np.abs(S_obp)**2
if plot_flag:
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
ax[0].plot(x_axis, w_obp/np.max(w_obp), ".-")
ax[1].plot(f_obp, 10*np.log10(S_obp))
ax[0].grid()
ax[0].set_xlabel("[km]")
ax[0].set_title("%s window (not normalized)" % kernel)
ax[1].grid()
ax[1].set_xlabel("freq [1/km]")
ax[1].set_ylabel("dB [m**2.km]")
plt.show()
return x_axis, w_obp, f_obp, S_obp
def get_obp_filter9(L_filt = 1, sampling_in = 0.025, f_axis = None, plot_flag = True, kernel="parzen"):
"""
Get the kernel shape and the spectral response of the filter type specified in the arguments
Arguments:
L_filt, Kernel length in km
sampling_in, Spatial spacing (in km) of the filter kernel
f_axis, Frequency axis where to compute the spectral response of the filter
plot_flag, Plot filter kernel and spectral response
kernel, Kernel to compute ('parzen' or 'bharris' -for Blackman-Harris)
"""
Nparzen = int(np.round(L_filt/sampling_in))
if Nparzen%2 == 0:
Nparzen += 1
print("Nb of points OBP kernel: %d" % Nparzen)
weights = np.array([0.07019711, 0.09767803, 0.12220272, 0.13814296, 0.14355835,0.13814296, 0.12220272, 0.09767803, 0.07019711])
PRF = 4332
N_burst = 9
v_sat = 7310.
H = 875.0e3
f_s = 300.0e6
v_ground = v_sat/(1+H/6377.0e3)
lambda_c = (consts.c/35.75e9)
incidence_deg = 3.5 # just used for the ground/range conversion for xtrack
d = N_burst*v_ground/PRF # along track posting before OBP averaging
w_ptr = np.zeros(512)
y_tab = np.arange(-200,201)*d+1.0e-3
ant = G_field_oneway(y_tab/H)**4
b_vals = np.arange(-4,5)
arg_b = 0.8/N_burst*b_vals[:,None]+2*v_sat/(lambda_c*H*PRF)*y_tab[None,:]
w_ptrantana_b = np.zeros((9,512))
w_ptrantana_c = np.zeros((9,512))
w_obp = np.zeros(512)
# Defines the PTR for each of the 9 beams
w_ptrantana_b[:,:len(y_tab)] = np.sin(np.pi*N_burst*arg_b)**2/np.sin(np.pi*arg_b)**2*(G_field_oneway(y_tab/H)**4)[None,:]
w_ptrantana_b = w_ptrantana_b/(np.sum(w_ptrantana_b, axis=1)[:,None]*d)
for ib in range(9) :
w_obp[:len(y_tab)]=w_ptrantana_b[ib,:len(y_tab)]
sampling_in=d/1000.
w_obp /= np.sum(w_obp)
x_axis = np.arange(-(len(w_obp)-1)/2, (len(w_obp)-1)/2 + 1)*sampling_in
# OBP spectrum
if type(f_axis) != np.ndarray:
if f_axis == None:
f_obp, S_obp = sg.freqz(w_obp, fs=1/sampling_in)
else:
f_obp, S_obp = sg.freqz(w_obp, fs=1/sampling_in, worN=f_axis)
else:
f_obp, S_obp = sg.freqz(w_obp, fs=1/sampling_in, worN=f_axis)
if ib== 0:
S_obpsum=weights[ib]*S_obp
else :
S_obpsum=S_obpsum+weights[ib]*S_obp
S_obp = np.abs(S_obpsum)**2
if plot_flag:
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
ax[0].plot(x_axis, w_obp/np.max(w_obp), ".-")
ax[1].plot(f_obp, 10*np.log10(S_obp))
ax[0].grid()
ax[0].set_xlabel("[km]")
ax[0].set_title("%s window (not normalized)" % kernel)
ax[1].grid()
ax[1].set_xlabel("freq [1/km]")
ax[1].set_ylabel("dB [m**2.km]")
plt.show()
return x_axis, w_obp, f_obp, S_obp
def compute_aliased_spectrum_2D(fx_in, fy_in, S_in, fsx, fsy, nrep=2):
"""
Given the main replica of a spectrum (S_in), compute its aliased version
for spatial sampling frequencies (fsx, fsy)
Arguments:
fx_in, horizontal frequency axis of the main spectrum replica
fy_in, vertical frequency axis of the main specumtr replica
S_in, main spectrum replica
fsx, target sampling frequency of the horizontal axis
fsy, target sampling frequency of the vertical axis
nrep, number of alias to compute (alias replicas at +-fs, +-2*fs... +-nrep*fs)
in1side, the spectrum provided
"""
fx_2side = fx_in
fy_2side = fy_in
S_2side = S_in
S = S_2side.copy()
Ly, Lx = np.shape(S_2side)
fs_idx = np.array([ np.argmin(np.abs(fy_2side - fsy)),\
np.argmin(np.abs(fx_2side - fsx)) ])
central_idx = np.array([ np.argmin(np.abs(fy_2side - 0)), \
np.argmin(np.abs(fx_2side - 0)) ])
shift_idx = fs_idx - central_idx
# fig, ax = plt.subplots()
# ax.plot(f_2side, 10*np.log10(S_2side), label='fundamental')
NX, NY = np.meshgrid( np.arange(-nrep, nrep+1), np.arange(-nrep, nrep+1))
for nx, ny in zip(NX.flatten(), NY.flatten()):
if (nx == ny == 0):
continue
shift_idx_nx = np.abs(nx)*shift_idx[1]
shift_idx_ny = np.abs(ny)*shift_idx[0]
if (shift_idx_nx >= Lx) or (shift_idx_ny >= Ly):
continue
#print("Adding replica %d, %d ..." % (nx, ny))
if nx >= 0:
col_indices_replica = slice(shift_idx_nx, Lx)
col_indices_in = slice(0, Lx-shift_idx_nx)
else: # inverse indices
col_indices_replica = slice(0, Lx-shift_idx_nx)
col_indices_in = slice(shift_idx_nx, Lx)
if ny >= 0:
row_indices_replica = slice(shift_idx_ny, Ly)
row_indices_in = slice(0, Ly-shift_idx_ny)
else: # inverse indices
row_indices_replica = slice(0, Ly-shift_idx_ny)
row_indices_in = slice(shift_idx_ny, Ly)
# print( row_indices_replica[0], row_indices_replica[-1], col_indices_replica[0], col_indices_replica[-1])
# print( row_indices_in[0], row_indices_in[-1], col_indices_in[0], col_indices_in[-1])
S[row_indices_replica][:, col_indices_replica] += S_2side[row_indices_in][:, col_indices_in]
return fx_2side, fy_2side, S