在MATLAB中对离散数值解y(t)做离散时间傅里叶变换（DTFT）

Question

运用四阶龙格库塔求解得到微分方程组的离散数值解y（t），现在需要在matlab中实现对y（t）做离散时间傅里叶变换求出y*（Ω）的代码，其计算定义如图所示。

Answer 1

基于Monster 组和GPT的调写：

下面是matlab代码示例，假设y(n)保存在向量y中，N为向量长度，Δt为时间步长：

% 计算DFT
Y = fft(y);
% 频移操作
Yshift = fftshift(Y);
% 计算频率向量
N = length(y);
fs = 1 / (N * Δt); % 采样频率
fnyquist = fs / 2; % 奈奎斯特频率
freq = linspace(-fnyquist, fnyquist, N);
% 计算离散时间傅里叶变换
y_star = Yshift;

其中，linspace函数用于生成频率向量，fft函数用于计算DFT。注意，如果要对y*(Ω)进行可视化，一般需要对结果进行幅值归一化，即y*(Ω)除以N，以便于比较不同信号的能量。

• 根据定义，离散时间傅里叶变换可以用离散傅里叶变换（DFT）实现，具体步骤如下：

• 对时间序列y(t)进行离散化，得到长度为N的向量y(n)，其中n=0,1,...,N-1，即y(n)=y(tn)，tn=nΔt，Δt为时间步长。

• 对y(n)进行N点DFT，得到频域向量Y(k)，其中k=0,1,...,N-1，即Y(k)=F[y(n)]。

• Y(k)的计算公式为： Y(k) = ∑ y(n) * e^(-j2πnk/N) ，其中n=0,1,...,N-1。

• 对Y(k)进行频移操作，将零频分量移到中心，以便于可视化和分析，得到频域向量Yshift(k)，其中

• Yshift(k) = Y(k-N/2)，当k < N/2时，Yshift(k) = Y(k+N/2)。

• 这里我们使用matlab的fftshift函数来实现频移操作。

• 计算频率向量ω(k)，其中ω(k) = 2πk/N，即ω(k)为频率k对应的角频率。

• 计算y*(Ω)，其中Ω为角频率，Ω=ω(k)。即：

• y*(Ω) = F[y(t)] = Yshift(k)。