diff --git a/docs/内积与基函数正交归一化的定义.md b/docs/内积与基函数正交归一化的定义.md new file mode 100644 index 0000000..afd179f --- /dev/null +++ b/docs/内积与基函数正交归一化的定义.md @@ -0,0 +1,88 @@ +### 1. 正交归一化条件 + +首先给出了一组基函数 \(\phi_m(s)\),要求它们满足正交归一化条件: + +\[ +\langle \phi_m(s), \phi_n(s) \rangle = \delta_{mn} \tag{27} +\] + +其中 \(\delta_{mn}\) 是克罗内克δ,表示当 \(m = n\) 时为1,否则为0。 + +--- + +### 2. 内积定义 + +内积被定义为: + +\[ +\langle \phi_m(s), \phi_n(s) \rangle = \frac{1}{2\pi i} \int_{i\mathbb{R}} \phi_m(s) \phi_n^*(s) ds \tag{28} +\] + +这里积分是在虚轴上,\(\phi_n^*(s)\) 表示复共轭。 + +--- + +### 3. 构造第一个基函数 \(\phi_1(s)\) + +首先考虑第一个基函数的归一化: + +\[ +\langle \phi_1(s), \phi_1(s) \rangle += \frac{1}{2\pi i} \int_{i\mathbb{R}} \left| \gamma_1 \right|^2 \frac{1}{(s+a_1)(-s+a_1^*)} ds \tag{29, 30} +\] + +经过计算,归一化条件变为: + +\[ += \frac{|\gamma_1|^2}{a_1 + a_1^*} \tag{31} +\] + +要使其归一化,设 \(Q_1(s) = \gamma_1\),则 \(\gamma_1\) 必须满足: + +\[ +\gamma_1 = \kappa_1 \sqrt{2 \Re(a_1)} +\] + +其中 \(\kappa_1\) 是任意单位模复数(即模为1的复数),最后得到: + +\[ +\phi_1(s) = \kappa_1 \sqrt{2\Re(a_1)} \frac{1}{s + a_1} \tag{32} +\] + +--- + +### 4. 构造第二个基函数 \(\phi_2(s)\) + +第二个基函数需要与第一个正交: + +\[ +\langle \phi_1(s), \phi_2(s) \rangle = \frac{1}{2\pi i} \int_{i\mathbb{R}} \phi_1(s) \phi_2^*(s) ds = 0 \tag{33} +\] + +推导得 \(\phi_2^*(s)\) 在 \(s = -a_1\) 时为零,因此: + +\[ +Q_2(s) = \gamma_2 (s - a_1^*) +\] + +--- + +### 5. 第二个函数的归一化 + +通过归一化条件确定 \(\gamma_2\): + +\[ +\langle \phi_2(s), \phi_2(s) \rangle = \frac{1}{2\pi i} \int_{i\mathbb{R}} \frac{|\gamma_2|^2}{(s + a_2)(-s + a_2^*)} ds = \frac{|\gamma_2|^2}{a_2 + a_2^*} \tag{34, 35} +\] + +--- + +## 总结 + +这张图的内容是关于如何通过积分定义的内积,构造一组在虚轴正交归一化的基函数。具体步骤包括: + +1. 给出正交归一化条件和内积定义。 +2. 构造第一个基函数并归一化。 +3. 构造第二个基函数,要求其与第一个正交,并归一化。 + +这些步骤可用于生成一组满足特定正交归一化条件的函数,常见于信号处理、系统理论或数学物理领域的谱基函数构造。 \ No newline at end of file diff --git a/test/lazy_evaulation.py b/test/lazy_evaulation.py index 446c5cc..c5ab46f 100644 --- a/test/lazy_evaulation.py +++ b/test/lazy_evaulation.py @@ -1,37 +1,46 @@ import numpy as np -class PhiBasisIterator: - """ - 逐步生成有理正交 (类 Muntz-Laguerre) 基函数 φ_p(s) / (或复对产生两列) 的迭代器。 +# ------------------------------------------------------------ +# 按论文公式 (9)(10)(11) 生成 Muntz–Laguerre 正交有理基 (解析形式): +# +# 给定稳定极点集合 {p_k} (Re(p_k)<0)。论文记法中使用 -a_k,其中 Re(a_k)>0。 +# 对应关系:p_k = -a_k ⇒ a_k = -p_k, Re(a_k)= -Re(p_k) >0 +# +# 连续内积意义下(沿 jω 轴积分)这些 φ_k 解析正交。离散频率采样后数值上 +# 可能偏离,可再用加权 QR 做数值再正交(可选)。 +# +# 公式在“稳定极点 p 表达”下的改写: +# (原) 实极点: φ_p(s) = sqrt(2 Re(a_p)) / (s + a_p) * Π (s - a_i^*)/(s + a_i) +# 变换 a_p = -p ⇒ Re(a_p)= -Re(p) = σ >0 且 (s + a_p) = (s - p) +# 且乘积 (s - a_i^*)/(s + a_i) = (s + p_i^*)/(s - p_i) +# ⇒ φ_p(s) = sqrt(-2 Re(p)) / (s - p) * Π_{i0 的 p 作为首: +# (原) φ_p = sqrt(2 Re(a_p)) (s - |a_p|)/[(s + a_p)(s + a_p^*)] * Π(...) +# (原) φ_{p+1}= sqrt(2 Re(a_p)) (s + |a_p|)/[(s + a_p)(s + a_p^*)] * Π(...) +# 代入 a_p=-p: +# Re(a_p)= -Re(p)=σ>0, (s + a_p) = (s - p), (s + a_p^*)=(s - p^*) +# |a_p| = |p| +# 乘积同上 ⇒ Π_{i0) - - 复极点以正虚部在前, 紧跟其共轭: a = -σ + jω, a* = -σ - jω - - 迭代: - 每次 __next__ 返回一个 list: - - 实极点 -> [φ_p(s)] - - 复成对 -> [φ_p^(1)(s), φ_p^(2)(s)] - 使用时可把返回列表中函数 append 到总基函数数组中。 - - 说明: - 这里实现的公式与用户原始代码意图相同, 但修正了: - 1) 原代码 sqrt(2*Re(a)) 在 Re(a)<0 (典型稳定极点) 下产生 NaN, 改为 alpha=-Re(a) >0. - 2) 原来在 __next__ 中使用 yield (非法); 改为标准迭代协议返回 list。 - 3) 复极点的乘积与后续迭代用的累计乘积 product 分开维护。 - 若需严格匹配论文中“多变量正交化”可再加数值正交 (QR) 步骤。 - - 注意: - 本实现假设 poles 已排序; 若不确定, 可预处理把正虚部极点放在其共轭之前。 - """ - def __init__(self, s, poles): - self.s = np.asarray(s, dtype=complex) # (Nf,) 或标量 - self.poles = list(poles) - self.k = 0 # 当前处理到第 k 个极点(或复对的首元素) +class MuntzLaguerreIterator: + def __init__(self, s: np.ndarray, stable_poles: list | np.ndarray): + """ + s: 复频率数组 (Nf,), s = j 2π f + stable_poles: 稳定极点列表 (Re<0). 复共轭对要求正虚部在前 (p, p*). + """ + self.s = np.asarray(s, dtype=complex) + self.poles = list(stable_poles) self.N = len(self.poles) - # 累计乘积 Π_{已完成的块} ( (s - a_i*)/(s + a_i) ) (复对时使用 (s - a)(s - a*) / (s + a)(s + a*)) + self.k = 0 + # 初始化乘积 Π_{i= self.N: raise StopIteration - a = self.poles[self.k] + p = self.poles[self.k] + if np.real(p) >= 0: + raise ValueError(f"极点必须在左半平面: {p}") - # 复对: 需要检查正虚部并且下一个是共轭 - if np.iscomplex(a) and np.imag(a) > 0: - if self.k + 1 >= self.N or not np.isclose(self.poles[self.k + 1], np.conj(a)): - raise ValueError(f"极点序列中复极点 {a} 未紧跟其共轭, 请排序.") - a_conj = self.poles[self.k + 1] - sigma = -np.real(a) # >0 - alpha = sigma - # (s + a)(s + a*) = (s + a)(s + conj(a)) - denom = (self.s + a) * (self.s + a_conj) - # |a| 用于生成两组分子 (参考原代码的构造方式) - r = np.abs(a) - scale = np.sqrt(2 * alpha) + # 复对首 (正虚部) + if np.iscomplex(p) and np.imag(p) > 0: + if self.k + 1 >= self.N: + raise ValueError("复极点缺少共轭") + pc = self.poles[self.k + 1] + if not np.isclose(pc, np.conj(p)): + raise ValueError("复极点未按 (p, p*) 顺序排列 (正虚部在前)") + sigma = -np.real(p) # >0 + scale = np.sqrt(2 * sigma) + r = np.abs(p) + denom = (self.s - p) * (self.s - pc) - numerator1 = scale * (self.s - r) - numerator2 = scale * (self.s + r) + # 两个基函数 + phi_p = scale * (self.s - r) / denom * self.product + phi_pc = scale * (self.s + r) / denom * self.product - phi1 = numerator1 / denom * self.product - phi2 = numerator2 / denom * self.product - - # 更新累计乘积用于后续阶次: - # 对复对的“整体”传统一致写成 ((s - a)(s - a*) / (s + a)(s + a*)) - self.product = self.product * ((self.s - a) * (self.s - a_conj)) / denom + # product 先乘 (s + p^*)/(s - p),再乘 (s + p)/(s - p^*) + self.product = self.product * (self.s + pc) / (self.s - p) + self.product = self.product * (self.s + p) / (self.s - pc) self.k += 2 - return [phi1, phi2] + return [phi_p, phi_pc] - else: - # 实极点 (允许传入 float 或实部为负的 complex(实数)) - a_real = np.real(a) - # 稳定极点: a_real < 0, 令 alpha = -a_real > 0 - if a_real >= 0: - raise ValueError(f"实极点需要在左半平面 (负实部), 得到 {a_real}") - alpha = -a_real - scale = np.sqrt(2 * alpha) - # Laguerre 风格: (s - alpha)/(s + alpha) - numerator = scale * (self.s - alpha) - denom = (self.s + alpha) - phi = numerator / denom * self.product - # 更新累计乘积 - self.product = self.product * (self.s - alpha) / (self.s + alpha) - self.k += 1 - return [phi] + # 复对次 (负虚部) —— 应该被首元素处理,出现表示顺序错误 + if np.iscomplex(p) and np.imag(p) < 0: + raise ValueError("检测到负虚部复极点但其共轭尚未处理,请将正虚部成员放在前面。") -# ------------------ 辅助函数: 批量生成全部基函数 ------------------ -def generate_phi_basis(s, poles): + # 实极点 + sigma = -np.real(p) + if sigma <= 0: + raise ValueError("实极点实部应为负 (稳定)。") + scale = np.sqrt(2 * sigma) + phi = scale / (self.s - p) * self.product + # 更新乘积 + self.product = self.product * (self.s + p) / (self.s - p) + + self.k += 1 + return [phi] + +def generate_muntz_laguerre_basis(s: np.ndarray, stable_poles: list | np.ndarray): """ - 返回一个 list, 其中包含依次生成的所有基函数 (按 real → 1 函数, complex pair → 2 函数). + 生成完整基函数列表: [φ_0=1, φ_1, φ_2, ...] """ - basis = [] - it = PhiBasisIterator(s, poles) - for block in it: + basis = [np.ones_like(s, dtype=complex)] + for block in MuntzLaguerreIterator(s, stable_poles): basis.extend(block) return basis -# ------------------ 示例与自测 ------------------ +# ---------- 可选:离散再正交 (加权 QR) ---------- +# def trapezoid_weights(freqs: np.ndarray): +# if len(freqs) == 1: +# return np.ones(1) +# df = np.diff(freqs) +# w = np.zeros_like(freqs, dtype=float) +# w[0] = 0.5 * df[0] +# w[-1] = 0.5 * df[-1] +# if len(freqs) > 2: +# w[1:-1] = 0.5 * (df[:-1] + df[1:]) +# return w + +def weighted_qr_from_basis(basis_cols: list[np.ndarray], weights: np.ndarray | None = None): + A = np.column_stack(basis_cols) # (Nf, M) + if weights is None: + sw = np.ones(A.shape[0]) + else: + sw = np.sqrt(weights.real) + Aw = sw[:, None] * A + Qw, R = np.linalg.qr(Aw) + Phi = Qw / (sw[:, None] + 1e-30) + return Phi, R # Raw = Phi R + +def check_discrete_orthogonality(Phi: np.ndarray, w: np.ndarray): + G = Phi.conj().T @ (w[:, None] * Phi) + off = np.max(np.abs(G - np.eye(G.shape[0]))) + return G, off + +def verify_orthonormal(Phi: np.ndarray, w: np.ndarray, atol=1e-10, rtol=1e-8): + """ + 返回: + G : Gram 矩阵 (Φ^H W Φ) + max_off : 最大非对角幅值 + diag_err : max |diag(G)-1| + passed : 是否满足阈值 + """ + G = Phi.conj().T @ (w[:, None] * Phi) + I = np.eye(G.shape[0]) + diag_err = np.max(np.abs(np.diag(G) - 1.0)) + max_off = np.max(np.abs(G - I + np.diag(np.diag(G) - 1.0))) + passed = (diag_err <= atol) and (max_off <= rtol) + return G, max_off, diag_err, passed + +def omega_weights(freqs_hz: np.ndarray): + """ + 基于 ω=2πf 的梯形法得到 w_ω = Δω/(2π),使得 + (1/2π) ∫_{-∞}^{∞} → Σ w_ω,k + """ + f = freqs_hz + if len(f) == 1: + return np.ones(1) + df = np.diff(f) + w_f = np.zeros_like(f) + w_f[0] = 0.5 * df[0] + w_f[-1] = 0.5 * df[-1] + if len(f) > 2: + w_f[1:-1] = 0.5 * (df[:-1] + df[1:]) + # dω = 2π df, (1/2π) * dω = df ⇒ 直接 w_f 就是 w_ω + return w_f # 已等价于 Δω/(2π) + +# ------------------ 示例 ------------------ if __name__ == "__main__": - # 示例极点: 1 个实极点 (-1), 一个复对 (-2+0.5j, -2-0.5j), 再一个复对 (-3+1.2j, -3-1.2j) - poles = [-0.005+500j, -0.005-500j, -0.012+1200j, -0.012-1200j] - f = np.linspace(1e9, 1e11, 200) # 频率 (Hz) - s = 1j * 2 * np.pi * f + # 示例稳定极点 (复对正虚部在前) + stable_poles = [ + -0.8e9, + -1.0e9 + 2.5e9j, + -1.0e9 - 2.5e9j, + -2.2e9 + ] + freqs = np.linspace(1e8, 8e9, 400) + s = 1j * 2 * np.pi * freqs - basis_funcs = generate_phi_basis(s, poles) - print(f"生成基函数数量 = {len(basis_funcs)}") # 实极点 1 →1, 两个复对 →2+2 共 5 个 - # 简单数值检查 (避免 NaN/Inf) - for i, bf in enumerate(basis_funcs): - if not np.all(np.isfinite(bf)): - print(f"第 {i} 个基函数存在非有限值") - else: - # print(f"phi[{i}] max|.|={np.max(np.abs(bf)):.3e}, min|.|={np.min(np.abs(bf)):.3e}") - print(f"phi[{i}] = {bf}") + basis = generate_muntz_laguerre_basis(s, stable_poles) + print("解析基函数数量 =", len(basis)) + print("前5个基函数示例 (每个前10个频点):") + for i in range(5): + print(f"φ_{i}:", basis[i][:10]) - # 验证是否为正交基 - for i, bf in enumerate(basis_funcs): - print(f"phi[{i}]*phi[{i}] = {np.vdot(bf, bf)}") # 自身内积 - for j in range(i+1, len(basis_funcs)): - bf2 = basis_funcs[j] - ip = np.vdot(bf, bf2) - print(f"phi[{i}] 与 phi[{j}] 内积 = {ip}") - if np.abs(ip) > 1e-6: - print(f"phi[{i}] 与 phi[{j}] 非正交, 内积={ip}") \ No newline at end of file + w = omega_weights(freqs) + Phi_num, R = weighted_qr_from_basis(basis, w) + Gram, off = check_discrete_orthogonality(Phi_num, w) + print("离散 Gram 最大非对角元素 =", off) + print("R 形状:", R.shape) + # 验证 Raw ≈ Phi R + raw = np.column_stack(basis) + err = np.max(np.abs(raw - Phi_num @ R)) + print("重构误差 ||Raw - Phi R||_∞ =", err) + + # 验证正交性 + print("离散 Gram 矩阵 (前5x5):") + print(Gram[:5, :5]) + + Gcheck, max_off, diag_err, ok = verify_orthonormal(Phi_num, w) + print(f"Diag 误差={diag_err:.3e}, Max off={max_off:.3e}, Orthonormal={ok}") + # 额外: 验证 R + # raw = Φ R => R ≈ Φ^H W raw (因为 Φ^H W Φ = I) + R_alt = Phi_num.conj().T @ (w[:,None] * raw) + print("R 差异 ||R - R_alt||_max =", np.max(np.abs(R - R_alt)))