傅里叶微分的差异 space

Question

我正在使用 Matlab 求解傅立叶 space 中的微分方程。然而，我遇到了一个问题：在对我的真实信号进行微分后，我得到了复杂的答案（这是不正确的）。

考虑一个在 x 上微分的例子（在傅里叶 space 中乘以 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。但问题同源。

谢谢

Answer 1

我想你只是少了一个 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 上使用希腊符号，我会的。但我想你明白我在说什么。

Answer 2

您需要考虑三件事：

时间微分对应于 将 DFT 乘以 1-exp(-1j*2*pi/N*fr)，其中 N 是信号周期，fr = 0:N-1 是频率样本.这源于 DFT 的时移属性（参见示例 here）。
这个微分应该是循环（见上文link），因为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

Answer 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 例程

中使用 'symmetric' 选项

傅里叶微分的差异 space

Discrepancy with differentiation in Fourier space

matlab

signal-processing

fft

dft