IFFT找到sin(x^2)？

Question

我想使用逆FFT来计算逆傅立叶变换。我发现我可以很容易地用平方可积函数来做，但不能用分布来做。

首先我建立一个波数矢量k和一个空间坐标x，

clear;

nx = 2^10;
L = 20;
dx = L/nx;
x = [0:nx-1]' * dx - L/2;

k = zeros(nx,1);
k(1:nx/2+1) = 2*[0:nx/2]/nx;
k(nx:-1:nx/2+2) = -k(2:nx/2);
k = k*pi/dx ;

接下来，我将验证两个示例是否一切正常:单元boxcar函数和sech：

% boxcar

Ghat = sin( k/2 )  ./ (k/2);
Ghat(1) = 1;
Gi = ifft(Ghat) / dx ;
Gi = ifftshift(Gi);
figure; plot(x,Gi);

% sech( x )

Ghat = pi*sech( pi*k /2 );
Gi = ifft(Ghat) / dx ;
Gi = ifftshift(Gi);
figure; plot(x,Gi,'o'); hold on;
analytical = sech(x);
plot(x,analytical,'-');

这两个看起来都很好。

这就是事情停止工作的地方：

% sin(x^2)

Ghat = -sqrt(pi)*sin( (k.^2-pi)/4 );
Gi = ifft(Ghat) / dx ;
Gi = ifftshift(Gi);
analytical = sin(x.^2);
figure;
plot(x,analytical,'-'); hold on;
plot(x,Gi,'o');

您会注意到，计算值看起来与预期的函数完全不同。

我不知道为什么这不管用。我在Wikipedia上唯一注意到的是sin(x^2)是一个分布，因此不是平方可积的。这是我的问题的根源吗？有解决方案吗？

Answer 1

我将给出一个python的例子，应该很容易翻译成Matlab。

import numpy as np
import matplotlib.pyplot as p
%matplotlib inline

dt=0.01
time = np.arange(0,10,dt) # 10 secs, sampled at 10 ms, so nyquist is 100 Hz
f=3 # hz
dat= np.sin(2*np.pi* f*time)

p.figure(figsize=(10,5))
p.subplot(221)
p.plot(time,dat)
p.subplot(222)
freqs = np.fft.fftfreq(len(time),d=dt)
#print(freqs)
spec=np.fft.fft(dat)
p.plot( freqs,np.abs(spec)  );


f=0.5
dat2= np.sin(2*np.pi* f*time**2)
p.subplot(223)
p.plot(time,dat2)
p.subplot(224)
p.plot( freqs,np.abs(np.fft.fft(dat2))  );

左边是正弦和带平方变元的正弦的时间行为(有效地，你有一个线性频率变化，一个upchirp。你可以把一个x想象成时间，另一个结合到你的频率中，它现在随着时间线性上升。在频谱中(显示正负频率，因为FFT被定义为复FT)，您可以看到正弦波的频率是稳定的，而正弦(x**2)扫描的频率范围。

描绘事物总是有帮助的。对于sin(x**2)，如果参数上升太快，很容易违反奈奎斯特定理，导致时间-domain (左)中的欠采样和频域(右)的混叠。这没有显示，只需尝试使用更高的基频。

​

​