傅立叶空间中的微分具有差异

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

考虑具有（在傅立叶空间中乘以ik）diffirentiation超过x一个例子：

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。但问题的起源相同。

谢谢

来源

2016-11-21 Mikhail Genkin

你需要采取三个方面考虑：

及时鉴别对应于由1-exp(-1j*2*pi/N*fr)，其中N是信号周期和fr = 0:N-1是频率样本的DFT倍增。这源于DFT的时移属性（参见例如here）。
这种差别应该是圆形（参见上面的链接），因为DFT本质上假定时间信号是周期性的。因此，在时域的第1微分样品将是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

来源

2016-11-21 20:03:38

谢谢！我总是在域[-N，N]中处理频率。这种方法可能有效。我只是不明白，为什么你必须乘以1-exp（i * 2pi * n/N）'。你有解释这个链接吗？ –

@Mikhail请参阅答案中的链接。它会告诉你信号的（循环）移位版本的DFT。要区分你需要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）=氙^（IWT）和

X'（T）= iwXe ^（ iwt）

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

如果我知道如何在SO上使用希腊符号，但我认为你明白我的意思。

来源

2016-11-21 20:03:35 Matt

我只是试图提供尽可能简单的例子。在你的代码中，你在ifft函数中使用了对称选项，消除了答案中的复数。不确定答案是否正确，你可能只是消除高频。 –

而如果我只是错过了2 * PI我不会得到复杂的答案，而不是真正的 –

我实际上检查：'对称'选项淘汰最高的频率。 –

路易斯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。正如马特所说，人们可以简单地使用“对称”选项作为例程

来源

2016-11-21 21:55:18

傅立叶空间中的微分具有差异

回答

相关问题