傅立葉空間中的微分具有差異

我使用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

傅立葉空間中的微分具有差異

回答

相關問題