fft 相位谱_Matlab中的FFT应用：导数逼近的谱方法

在上一篇文章中，我们介绍了 Matlab 中 FFT 计算的表达式和基本操作。

其中最值得注意的就是「指标的习惯」

我们也介绍了 FFT 可以用来帮助我们快速计算「傅里叶系数的数值近似」

不过这部分内容讲的不够详细，因此本文先对此进行一部分的补充。

回顾

周期函数傅里叶系数的计算

将积分区间

进行等距剖分

使用矩形公式（梯形公式同理）进行数值离散

重点

周期函数

「采样」时 「取左端点不取右端点」

该部分的「主题」是详细测试一下：

Matlab 中 FFT 在使用的过程中的结果是否和我预想相同。

如果不同，我需要调整一下自己对 DFT 的理解。

混叠效应

如果我们的「采样频率低于两倍的目标频率」的话，将会发生混叠效应。

我们从离散傅里叶变换的公式上来看一下这个性质

其中

所以对于下面两个周期函数

和

它们的采样

是相同的。

「具体的例子」如下，取

灰色的线条为

，蓝色的线条为

从「等距采样节点」（红色点）来看，我们无法区分两者。

所以，从这个角度，我们采样频率应该至少等于目标频率。那为什么要两倍呢？

比如，对于周期函数

我们在

上采样 40 个点，即取

，通过

以及

来生成采样向量 X，然后用

Y=fft(X) 来计算对应的离散傅里叶变换。

按照之前的分析，向量

的第

个分量应该是

的一个「数值离散近似」。因为

有两个模态，

• 一个频率为
• 一个频率为

我们写个程序测试一下，运行结果如何呢？

N = 40;
x = linspace(0,2*pi,N+1);
x = x(1:end-1);
y = exp(1i*30*(x))+exp(-1i*10*(x));
frequency = fft(y);
figure;
set(gcf,'unit','centimeters','position',[10,10,24,15])
plot(1/N*real((frequency)),'-o','LineWidth');

结果「居然只有一个波峰」，而且强度为

解释

之前我们讲过离散傅里叶变换得到的

的周期性

因此，在

的时候，

的模态和

的混合了。

我们可以测试一下

的情形

可以看到

• 的模态给出了在
  的一个波峰
• 的模态，等价于
  ，给出了
  处的波峰

因此，为了能够区分出

和

这两个模态（以及和所有

），我们需要选择

至少为两倍的目标频率。

也因此，我们需要使用 fftshift 函数将得到的频率超过

的部分等价转移到

部分。

平移与其他周期

思考这个问题的「起源」

有些书本中的采样舍去了左端点（而不是右端点）

因此，利用

区间的等距剖分（不包含左端点），我们有

所以对这个向量执行离散傅里叶变换的话 Y1 = fft(X1)

给出的其实是

「向左平移」

之后得到函数的傅里叶变换系数

如果我们记向左平移得到的函数为

，那么

所以 Y1 这个向量逼近的是积分

特别地，我们取

则

特殊的，我们取

，这样，多出来的因子为

这个时候我们再绘制 Y1 的图像，代码如下

N = 60;
x = linspace(0,2*pi,N+1);
x = x(2:end);
y = exp(1i*30*x);
frequency = fft(y);
figure;
set(gcf,'unit','centimeters','position',[10,10,24,15])
plot(1/N*real((frequency)),'-o','LineWidth',2,'MarkerIndices',31);

结果如图

可以发现

• 结果符合我们的分析，「对应模态的系数」变成了
• 如果我们使用 fft 来进行积分的数值计算，那么采样过程需要「舍去右端点」，而不是左端点。否则计算的是平移过后函数的傅里叶系数，会额外出现一个相位因子

其他周期

假如我们的函数

是

周期的，我们通过仿射变换将

映射到单个周期

则

是

周期的。

因此我们对它采样，得到

则

得到的其实是

上的等距剖分（保留左端点）

下面我们将讲到本文的重点。如何利用 fft 来对导数进行「具有谱精度的近似」。

导数近似

对函数

在

处的导数，我们可以通过差分来近似，例如下面的两种

这些都叫做「局部的近似」，它们都只用到了

附近的点。

下面我们介绍一种对周期函数导数的全局近似，这是一种谱方法（收敛很快）

我们仍然假定函数

是

周期的函数，在

上有等距剖分

我们记这些点上的函数的采样值为

「导数逼近的思路」如下：

• 先利用
  来构造函数
  的一个逼近函数
• 使用
  来作为
  的近似

逼近函数

的构造

首先，因为

是

周期函数，假定有足够的光滑性，那么它有傅里叶级数展开

回忆离散傅里叶变换（也就是 FFT）做的事情是通过

来计算了

的近似值。

因此，我们可以「将傅里叶展开截断」。令

作为

的一个近似。

此时，我们可以计算

的导数来获得

导数的近似

特别地，我们关心

在采样点

上的导数值，近似表示为

这是什么呢？对应上了 ifft （逆快速傅里叶变换）

怎么个对应关系呢？利用

得到

拆分原因

上面的式子之所以这么拆分是因为这样，

的

的取值就和

ifft 有了对应，即

因此，结合上 fft 的运算，我们总结一下「导数离散的计算思路」：

• 首先周期内采样，计算采样点的函数值
• 通过 fft 计算各个模态的傅里叶系数的数值近似
• 乘上对应的模态
• 逆变换 ifft

具体计算过程中会涉及到

• 常数因子
  的乘除（最后都抵消了）
• 模态
  的排序

我们使用下面的代码片段来测试一下

N = 100;
y = linspace(0,2*pi,N+1);
y = y(1:end-1);
X = exp(sin(y));
Y = fft(X);
% 计算数值导数
deri_Y = ifft(Y.*[1i*(0:N/2), 1i*(-N/2+1:-1)]);
% 计算精确导数
deri_exact = cos(y).*X;
% 绘图
figure;
plot(y,deri_Y,'r-o'); hold on;
plot(y,deri_exact,'b-.');
legend('数值导数','精确导数');

这显示了谱方法逼近导数具有非常好的效果。

「两个注意点」

• 「第一个注意点」：程序的第三行，我们的采样点是从左端点开始。这里修改成

y = y(2:end);

同样没有问题。

我们之前说，在使用 fft 计算频率的模值的时候，为了对应上具体的公式，我们需要从左端点开始计数，否则会出现一个相位因子。

但是此处由于是计算导数的近似，我们在 fft 之后还要进行一个 ifft，所以这个因子会消去。

同样道理，我们在上述代码中，在 fft 之后也没有像之前那样除掉

，它和本来在

ifft 中要处理的

一样抵消了。

我们展示一下最终得到的结果，最左侧的点有一个微小的平移。

• 「第二个注意点」

ifft(Y.*[1i*(0:N/2), 1i*(-N/2+1:-1)]);

由于离散傅里叶变换公式的特殊性

在

是实数的情况下，如果

• ，则计算出来
  是实数
• 此外，

所以

。 举个例子就是

所以在 ifft 的时候

前两项加起来正好是实数，只剩下最后一项。

我们希望它不起贡献。

实际上

的实数部分在采样节点上确实导数为

，因此我们可以修改我们的近似。

在程序中表现为

deri_Y = ifft(Y.*[1i*(0:N/2-1),0, 1i*(-N/2+1:-1)]);

总结

本文补充说明了使用 fft 来作为傅里叶级数数值离散计算工具的一些细节和例子。

然后讨论了使用 fft 和 ifft 来进行导数谱精度离散的具体操作与注意事项。

下期，我们将使用这一套离散来求解微分方程——谱方法。