与傅里叶空间微分的差异答案

【问题标题】：Discrepancy with differentiation in Fourier space与傅里叶空间微分的差异
【发布时间】：2016-11-21 19:30:46
【问题描述】：

我正在使用 Matlab 求解傅里叶空间中的微分方程。但是，我遇到了一个问题：在区分我的真实信号后，我得到了复杂的答案（这是不正确的）。

考虑一个对x 进行微分的示例（在傅立叶空间中乘以ik）：

 a=rand(6,1).';
 fr=fftshift(-3:1:2);
 ifft(1i*fr.*fft(a))

输出很复杂。我知道为什么会这样：我们的频谱是-3,-2,-1,0,1,2。因此，最高频率没有对（我们有 -3，但没有 3）。我想知道如何解决它。

如果我们考虑一下，从技术上讲，最高频率的贡献是非零的。如果频率 -3 上的傅立叶幅度为 c0，这意味着实际上我们在频率 -3 和 3 上具有幅度c0/2，因此经过微分我们得到：

(c0/2)*i*(-k)*exp(-ikx)+(c0/2)*i*(k)*exp(ikx)=kc0*sin(kx)

我很好奇，如何实现正确的微分。我的问题是 2D，所以我使用 fft2 和 ifft2。但问题是同源的。

谢谢

【问题讨论】：

标签： matlab signal-processing fft dft

【解决方案1】：

你需要考虑三件事：

时间微分对应于乘以 DFT 乘以1-exp(-1j*2*pi/N*fr)，其中N 是信号周期，fr = 0:N-1 是频率样本。这源于 DFT 的时移特性（例如参见 here）。
这种区分应该是循环的（参见上面的链接），因为 DFT 固有地假定时间信号是周期性的。所以时域中的第一个微分样本是a(1)-a(N)，第二个是a(2)-a(1)，以此类推。
由于浮点精度，您可能会得到一个非常小的虚部。

所以，代码应该是：

a = rand(6,1).';
N = numel(a);
fr = 0:N-1;
a_diff_fr = ifft((1-exp(-1j*2*pi/N*fr)).*fft(a));

检查：

>> a_diff_fr % imag part should be small
a_diff_fr =
  Columns 1 through 5
  -0.5490 - 0.0000i   0.3169 - 0.0000i  -0.5662 + 0.0000i   0.6851 + 0.0000i  -0.5155 - 0.0000i
  Column 6
   0.6287 + 0.0000i

>> real(a_diff_fr) % real part only
ans =
   -0.5490    0.3169   -0.5662    0.6851   -0.5155    0.6287

>> a([1 2:N])-a([N 1:N-1]) % circular differentiation
ans =
   -0.5490    0.3169   -0.5662    0.6851   -0.5155    0.6287

【讨论】：

谢谢！我总是处理域 [-N,N] 中的频率。这种方法可能有效。我只是不明白，为什么你必须乘以1-exp(i*2pi*n/N)。你有解释这个的链接吗？
@Mikhail 查看答案中的链接。它告诉您信号的（循环）移位版本的 DFT。为了区分你需要减去 1

【解决方案2】：

我认为您只是在其中缺少 2 pi。这里是一个使用x = 2 cos(2*pi*t)的例子

Tp = 10; % sample length
deltaTime = Tp / 200; % time step
time = 0:deltaTime:Tp; % time
x = 2*cos(2*pi*time); % function
plot(time, x)
fMax = 1/deltaTime/2; % maximum frequency
fMin = 1/Tp; % lowest observable frequency
xfft = fft(x) ./ (length(time)/ 2); % fft scaled to original amplitude, in case you want to plot it.
freq = -1*fMax:fMin:fMax; % frequencies of the fft
xd_fft = xfft .* fftshift(freq) * 1i*2*pi; % note the extra 2 pi in here.
xd = ifft(xd_fft, 'symmetric') * (length(xd_fft)/ 2); % reverse the scaling and take the ifft.
xd2 = -4*pi*sin(2*pi*time);
plot(time, xd2, '.')

我这样做的方程式是：

x(t) = Xe^(iwt) 和

x'(t) = iwXe^(iwt)

回想一下 w = 2*pi*f。

如果我知道如何在 SO 上使用希腊符号，我会的。但我想你明白我在说什么。

【讨论】：

我只是尝试提供尽可能简单的示例。在您的代码中，您在 ifft 函数中使用了对称选项，这消除了答案中的复数。不确定答案是否正确，您可能只是通过高频。
如果我错过了 2*pi 我不会得到复杂的答案而不是真实的
我实际检查过：“对称”选项消除了最高频率。
但这是你真正需要做的，所以从技术上讲你的答案也是正确的。我发布了我的解释以防万一
@MikhailGenkin 你是对的，symmetric 选项清理了答案，因为它强制 Matlab 使用对称形式。但是，您的原始方程式缺少2 pi。否则，您的导数是错误的。要对其进行测试，只需更改我示例的第 10 行。如果您从该行中删除 2 pi，则答案是错误的。

【解决方案3】：

Luis Mendo 为我的问题提出了正确的解决方案。我还想出了如何解决我的方法：想法是，正弦分量在高频上始终为零，只有余弦谐波可见：

signal=sin(2.*(linspace(0,2*pi*(1-1/4),4)));
q=fftshift(fft(signal))./4

这里 q 为零。但如果使用 cos(2x) 信号：

signal=cos(2.*(linspace(0,2*pi*(1-1/4),4)));
q=fftshift(fft(signal))./4

这里 q(1)=1。因此，在我的方法中，只要正弦谐波不可见，我必须在乘以 ik 后将最高谐波的虚部设置为 0。正如马特建议的那样，可以简单地在 ifft 例程中使用“对称”选项

【讨论】：