Python中非均匀时间序列的DFT计算与对应频率数组求解

你写的DFT函数是标准离散傅里叶变换实现，仅适用于采样间隔固定的均匀时间序列。你的时间数组t采样间隔不固定，不满足标准DFT的应用前提，因此无法得到正确结果。

1. 改用非均匀离散傅里叶变换（NUDFT）实现 #

标准DFT的e指数项基于均匀采样的时间点t_n = n * dt推导，非均匀采样需要直接代入实际采样时间t计算，正确实现如下：

import numpy as np

def nudft(t, x, f):
    """
    非均匀采样信号的离散傅里叶变换实现
    参数:
        t: 实际采样时间数组，长度N
        x: 对应时间点的信号值（这里就是电场E），长度N
        f: 要计算的目标频率数组，长度M
    返回:
        X: 对应频率点的傅里叶变换结果，长度M
    """
    N = len(t)
    M = len(f)
    # 构造e指数项：维度为[M, N]
    e = np.exp(-2j * np.pi * f.reshape((M, 1)) * t.reshape((1, N)))
    X = np.dot(e, x) / N  # 可根据需求决定是否保留除以N的归一化步骤
    return X

2. 构造对应频率数组 #

你可以按照以下逻辑构造符合需求的频率数组：

• 总采样时长：T = t[-1] - t[0]，默认频率分辨率取df = 1/T
• 最小采样间隔：dt_min = np.min(np.diff(t))，最高分析频率建议不超过对应奈奎斯特频率f_max = 1/(2*dt_min)
• 生成频率数组代码：

t = np.array([1, 5, 6, 8.5, 12, 20, 21.5, 22.3, 27, 30])
T = t[-1] - t[0]
df = 1 / T
dt_min = np.min(np.diff(t))
f_max = 1 / (2 * dt_min)
# 生成0到f_max的频率数组，步长为df
f = np.arange(0, f_max, df)

3. 调用示例 #

# 替换为实际的电场序列E
E = np.array([你的电场测量值])
E_dft = nudft(t, E, f)

内容的提问来源于stack exchange，提问作者Subhadip Saha