Skip to content
This repository has been archived by the owner on May 31, 2024. It is now read-only.

Commit

Permalink
wip(propagate_DFT): adjust GUM_iDFT to match np.irfft(..., n=...)
Browse files Browse the repository at this point in the history
related to #310
  • Loading branch information
mgrub committed Mar 7, 2023
1 parent 318b56a commit 9060b40
Showing 1 changed file with 18 additions and 15 deletions.
33 changes: 18 additions & 15 deletions src/PyDynamic/uncertainty/propagate_DFT.py
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,7 @@
is_2d_matrix,
is_vector,
number_of_rows_equals_vector_dim,
trimOrPad_ND,
)


Expand Down Expand Up @@ -388,47 +389,49 @@ def GUM_iDFT(
ValueError
If Nx is not smaller than dimension of UF - 2
"""
# (complex) input length
N_in = F.size // 2

# default output length, assumes even output length
N = UF.shape[0] - 2

if Nx is None:
Nx = N
elif Nx > UF.shape[0] - 2:
raise ValueError(
f"GUM_iDFT: Nx if provided is expected to be smaller or equal to number "
f"of rows and columns of UF - 2, but UF is of shape {UF.shape} and Nx ="
f" {Nx}."
)


# TODO: continue here, length=Nx probably wrong
F = trimOrPad_ND(F, length=Nx, real_imag_type=True)
UF = trimOrPad_ND(UF, length=Nx, real_imag_type=True)

beta = 2 * np.pi * np.arange(Nx) / N

# calculate inverse DFT; Note: scaling factor 1/N is accounted for at the end
x = np.fft.irfft(F[: N // 2 + 1] + 1j * F[N // 2 + 1 :])[:Nx]
x = np.fft.irfft(F[: N_in] + 1j * F[N_in :], n=Nx)
if not isinstance(Cc, np.ndarray): # calculate sensitivities
Cc = np.zeros((Nx, N // 2 + 1))
Cc = np.zeros((Nx, N_in))
Cc[:, 0] = 1.0
Cc[:, -1] = np.cos(np.pi * np.arange(Nx))
for k in range(1, N // 2):
Cc[:, k] = 2 * np.cos(k * beta)

if not isinstance(Cs, np.ndarray):
Cs = np.zeros((Nx, N // 2 + 1))
Cs = np.zeros((Nx, N_in))
Cs[:, 0] = 0.0
Cs[:, -1] = -np.sin(np.pi * np.arange(Nx))
for k in range(1, N // 2):
Cs[:, k] = -2 * np.sin(k * beta)

# calculate blocks of uncertainty matrix
if len(UF.shape) == 2:
RR = UF[: N // 2 + 1, : N // 2 + 1]
RI = UF[: N // 2 + 1, N // 2 + 1 :]
II = UF[N // 2 + 1 :, N // 2 + 1 :]
RR = UF[: N_in, : N_in]
RI = UF[: N_in, N_in :]
II = UF[N_in :, N_in :]
# propagate uncertainties
Ux = np.dot(Cc, np.dot(RR, Cc.T))
Ux = Ux + 2 * np.dot(Cc, np.dot(RI, Cs.T))
Ux += np.dot(Cs, np.dot(II, Cs.T))
else:
RR = UF[: N // 2 + 1]
II = UF[N // 2 + 1 :]
RR = UF[: N_in]
II = UF[N_in :]
Ux = np.dot(Cc, _prod(RR, Cc.T)) + np.dot(Cs, _prod(II, Cs.T))

if returnC:
Expand Down

0 comments on commit 9060b40

Please sign in to comment.