使用FFT处理ode23离散ODE解数据：PSD计算与绘图疑问

我正在求解离散化ODE系统，采用N=256个网格点，域长L=256，空间网格为等间距的-L/2:1:L/2-1。使用Matlab的ode23求解器后，得到的解矩阵u尺寸为2177 x 256，根据文档说明，u的每一行对应列向量t中的一个时间点，即每一行是一个时间切片。

我希望计算功率谱密度$|\hat{u}(t)|^2$，并将其绘制在范围为$-\pi$到$\pi$的离散波数$\kappa$上。

根据fft(X,n)的文档：

若X为矩阵，则每列按向量情况处理；若X为向量，则fft(X)返回该向量的傅里叶变换。

我有以下疑问：

1. 是否应设置n = length(t) = 2177，并使用uhat=fft(u,n)？此时功率谱密度为PSD = uhat.*conj(uhat)/n（除以n进行归一化），还是应该使用fft(u',n)？
2. 如何将功率谱密度与波数进行绘图？

问题1：FFT的正确调用方式 #

你需要对每个时间切片（即矩阵u的每一行）做空间傅里叶变换——因为目标是获取每个时刻下空间分布的功率谱，而非时间维度的频谱。

• 不要用n=length(t)：length(t)=2177是时间点数量，你的空间网格点为256个，FFT应针对空间维度执行。
• 正确调用方式：Matlab的fft默认对列操作，要对每一行做FFT需指定维度参数：uhat = fft(u, [], 2);。其中第二个参数留空表示使用当前维度长度（256），第三个参数2表示对矩阵的第二个维度（行）执行FFT。
• 归一化规则：功率谱需除以空间网格点数量N=256，而非时间点数量。计算式为：PSD = abs(uhat).^2 / N;（abs(uhat).^2与uhat.*conj(uhat)效果完全一致）。
• 禁止使用fft(u',n)：转置矩阵后对列操作，得到的是每个空间点的时间频谱，与需求完全不符。

问题2：功率谱与波数的绘图 #

步骤1：生成对应-π到π的离散波数向量 #

基于你的域长和网格点参数，波数间隔为$\Delta k = 2\pi/L = \pi/128$，可直接生成目标范围的波数向量：

N = 256;
L = 256;
dk = 2*pi/L;
k = dk*(0:N-1) - pi; % 生成范围为-π到π的波数

或通过fftshift调整默认FFT波数范围：

k_default = dk*(0:N-1); % 默认0到2π的波数
k = fftshift(k_default) - pi;

步骤2：单时刻功率谱绘图 #

选择某个时间点（比如第100个时刻）绘制功率谱：

idx = 100; % 指定时间点索引
plot(k, PSD(idx,:));
xlabel('波数 \kappa');
ylabel('功率谱密度 |\hat{u}|^2');
xlim([-pi, pi]);

步骤3：全时间维度功率谱热力图 #

若要展示所有时刻的功率谱分布，可使用pcolor或imagesc：

pcolor(t, k, PSD'); % 转置PSD使波数对应y轴、时间对应x轴
shading flat;
xlabel('时间 t');
ylabel('波数 \kappa');
colorbar;

内容的提问来源于stack exchange，提问作者KZ-Spectra