フーリエ空間での微分による相違

Question

私は、Matlabを使用してフーリエ空間で微分方程式を解いています。しかし、私は問題があります。私の実際の信号を差別化した後、私は複雑な答えを得ます（これは間違っています）。

はx（フーリエ空間におけるik乗算）上diffirentiation有する例を考える：

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

出力は複雑です。私はなぜそれが起こるかを考えました：私たちの周波数スペクトルは-3,-2,-1,0,1,2です。したがって、最も高い周波数のペアはありません（-3がありますが、3はありません）。私はそれを修正する方法が不思議です。

私たちはそれについて考えると、最高頻度から技術的に非ゼロ寄与があります。

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

私は好奇心、どのように実装する：周波数のフーリエ振幅が-3分化した後、我々が得るように、現実には、我々は、周波数-3と3の振幅c0/2を持っていることを、意味c0であれば、正しい差別化。私の問題は2Dだから、私はfft2とifft2を使う。しかし、問題は同じ起源です。

おかげ

Answer 1

は、あなたのアカウントに3つのことを取る必要があります。

• がNが信号周期であり、fr = 0:N-1は周波数サンプルです1-exp(-1j*2*pi/N*fr)によってDFTを掛けるまでの時間に対応して差別化。これは、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 

Answer 2

私はあなたがそこにある2パイが欠けていると思います。ここで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を思い出してください。

私はギリシャのシンボルを使用する方法を知っていました。しかし、私はあなたが私が言っていることを得ると思います。

Answer 3

ルイス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に設定する必要があります。 Mattが提案したように、ifftルーチンで単に 'symmetric'オプションを使うことができます