MATLAB FFT 频率的顺序

计算科学 matlab 傅立叶分析
2021-12-10 08:13:59

该 wikibook指出 MATLAB 的输出FFT与按以下顺序排列的波数相对应:

k={0,1,...,n2,n2+1,n2+2,...,1}

但是,在同一页面上的示例代码中,波数编码为

k = [0:n/2-1 0 -n/2+1:-1];

这与第一个相同,但与n/2-wavenumber(“最大wavemnumber”)替换为0. 似乎很奇怪,他们会包括0两次。

如wikibook中所述,通过傅里叶变换获取导数似乎需要正确的顺序。哪一个是正确的,MATLAB 是否在任何地方记录了这一点?

4个回答

我想扩展我的评论并以比原始示例更易于理解的方式重新编写您引用的示例,并解释为什么fft以这种方式返回系数。

作为参考,示例的 fft 部分是:

Nx = size(x,2);
k = 2*pi/(b-a)*[0:Nx/2-1 0 -Nx/2+1:-1];
dFdx = ifft(1i*k.*fft(f));
d2Fdx2 = ifft(-k.^2.*fft(f));

我在它下面直接添加了另一段代码:

Nx = size(x,2);
k = 2*pi/(b-a)*(-Nx/2:Nx/2-1);
dFdxp = ifft(ifftshift(1i*k.*fftshift(fft(f))));
d2Fdx2p = ifft(ifftshift(-k.^2.*fftshift(fft(f))));

并将两段代码都包装在 a 中以tic; toc进行粗略的计时。在更易读的格式中,第二种方法使用:

ckf = fftshift(fft(f));
ckdf = 1i*k.*ckf;
df = ifft(ifftshift(ckdf));

第一个区别是第二个例子更直观k这是第二个例子的主要优点,因为 k 现在是我们考虑它们的形式。在第二行和第三行我不得不fftshift在 call to 周围添加一个 call to fft,然后ifftshift直接在 call to 里面添加一个 call to ifft这些额外的函数调用将系数从计算机使用它们所需的系数重新排序为人类通常对它们的看法。

第二个例子的问题是,虽然k对我们来说更直观,但这留下了用于求解和反转fft形式的内部矩阵并不那么有利。因此,要么我们必须通过调用fftswitchand来切换顺序,ifftswitch要么必须将其硬编码到fft函数中。这不太容易从用户那里出错(假设他们不熟悉 fft 的工作原理,就像许多人一样),但是您在运行时付出了代价。

正如我之前所说,我在两个块周围添加了计时调用以进行比较,并运行了多个 N。计时结果是:

N =     1000,  Ex1 = 0.000222 s,   Ex2 = 0.007072 s
N =    10000,  Ex1 = 0.001576 s,   Ex2 = 0.003506 s
N =   100000,  Ex1 = 0.023857 s,   Ex2 = 0.034051 s
N =  1000000,  Ex1 = 0.213816 s,   Ex2 = 0.406250 s
N = 10000000,  Ex1 = 4.555143 s,   Ex2 = 7.102348 s

如您所见,来回切换值的行为会大大减慢该过程,尤其是在低 N 时(慢 30 倍)。这只是一个示例,您的计算机可能会根据内存速度、处理器内核/速度等因素显示出略有不同的趋势,但这只是说明了这一点。fft输出令人困惑的原因是因为它为您节省了相当多的计算时间。

您关于波数“替换”的问题相当棘手。一般来说,这种波数修改并不是为了挽救失败,正如一些人在这里所建议的那样,而是为了尊重例如某些微分算子的分析特性而设计的。我很惊讶在 Trefethen 的 Spectral Methods 中找不到相关的讨论。为了继续,我假设您关心基于 FFT 的谱微分,并且您正在对偶基数域执行变换。

对于奇数导数,经验法则是设置

k = [0:n/2-1 0 -n/2+1:-1];

甚至对于导数,设置

k = [0:n/2 -n/2+1:-1];

如果您正在处理任何填充,则奈奎斯特频率的处理同样繁琐。这里有一篇优秀的文章证明和评论这些主题

似乎两者都是正确的;这就是为什么。n/2模式别名为等间距网格上的零模式n点(因为该模式在所有网格点处精确地通过零)。因此两者是无法区分的。我不确定为什么在代码中(以及在 Trefethen 的书中)他们费心替换n/20; 也许是因为这样可以挽救一些失败。

有关别名的简短说明和演示,请参阅我的 IPython 笔记本

如果N是偶数(我们通常使用 fft 来快速),为什么要设置N/2模式归零?好吧,如果信号是真实的 xR,取傅里叶逆变换,模 N/2可能产生一个虚部,因为没有共轭 N/2组件取消 N/2模式贡献,因此我们将其设置为零。