Если бы найти грех (x ^ 2)?

Question

Я хочу использовать обратное БПФ для вычисления обратных преобразований Фурье. Я обнаружил, что могу легко сделать это с помощью интегрируемых с квадратом функций, но не с помощью распределений.

Сначала я устанавливаю вектор волновых чисел 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 ;

Затем я проверяю, что все работает на двух примерах: функция unit 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');

Вы заметите, что вычисленные значения не похожи на целевую функцию.

Я не знаю, почему это не сработает. Единственное, что я замечаю в Википедии , - это то, что 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 является временем, а другой - частью вашей частоты, которая теперь линейно растет со временем. В спектре (показывая положительную и отрицательную частоту, так как БПФ определяется как сложный FT), вы можете видеть, что частота синусоидальной волны стабильна, тогда как sin (x ** 2) охватывает диапазон частот.

Это всегда помогает составить график. x ** 2) вы можете легко нарушить теорему Найквиста, если аргумент возрастает слишком быстро, что приводит к недостаточной дискретизации во временном домене (слева) и наложению алиасов в частотной области (справа). Это не показано, просто попробуйте с более высокой базой частота.