`numpy.diff` и` scipy.fftpack.diff`, дающие разные результаты при дифференцировании

Question

Я пытаюсь вычислить производную некоторых данных, и я пытался сравнить выход конечной разности и выход спектрального метода. Но результаты очень разные, и я не могу понять, почему именно.

Рассмотрим пример кода ниже

import numpy as np 
from scipy import fftpack as sp 
from matplotlib import pyplot as plt 
x = np.arange(-100,100,1) 
y = np.sin(x) 

plt.plot(np.diff(y)/np.diff(x)) 
plt.plot(sp.diff(y)) 

plt.show() 

Это выводит следующий результат

Оранжевый выход, являющуюся fftpack выход. Не обращайте внимания на тонкости, это только ради примера.

Итак, почему они такие разные? Разве они не должны быть (примерно) одинаковыми?

Я уверен, что различные амплитуды могут быть скорректированы с ключевым словом fftpack.diff, но я не могу понять, какой правильный период (я думал, что это должно быть period=1, но это не работает).

Кроме того, как я могу использовать собственное спектральное дифференцирование с помощью numpy?

Answer 1

Функция scipy.fftpack.diff вычисляет производную, но предполагает, что вход является периодическим. Аргумент period дает период (т. Е. Общую длину интервала x) входной последовательности.

В вашем случае это len(x)*dx, где dx = x[1] - x[0].

Вот некоторый код, который описывает простую (центрированную) конечную разницу (в синем) и результат diff с использованием аргумента period (красным). Переменные x и y такие же, как те, которые используются в вашем коде:

In [115]: plt.plot(0.5*(x[1:]+x[:-1]), np.diff(y)/np.diff(x), 'b') 
Out[115]: [<matplotlib.lines.Line2D at 0x1188d01d0>] 

In [116]: plt.plot(x, sp.diff(y, period=len(x)*(x[1]-x[0])), 'r') 
Out[116]: [<matplotlib.lines.Line2D at 0x1188fc9d0>] 

In [117]: plt.xlabel('x') 
Out[117]: <matplotlib.text.Text at 0x1157425d0> 

Обратите внимание, что если ваш вклад фактически не периодично, производная вычисляется diff будет неточным вблизи концах интервал.

Вот еще один пример, используя более короткую последовательность, которая содержит только один полный период синусоидальной функции в интервале [0, 1]:

In [149]: x = np.linspace(0, 1, 20, endpoint=False) 

In [150]: y = np.sin(2*np.pi*x) 

In [151]: plt.plot(0.5*(x[1:]+x[:-1]), np.diff(y)/np.diff(x), 'b') 
Out[151]: [<matplotlib.lines.Line2D at 0x119872d90>] 

In [152]: plt.plot(x, sp.diff(y, period=len(x)*(x[1]-x[0])), 'r') 
Out[152]: [<matplotlib.lines.Line2D at 0x119c49090>] 

In [153]: plt.xlabel('x') 
Out[153]: <matplotlib.text.Text at 0x1197823d0> 

Answer 2

1 рад довольно грубый шаг для разностной аппроксимации, и вы должны хотеть целое число периодов в наборе данных

x = np.arange(-200,200,1) 
y = np.sin(np.pi/50*x) 

plt.plot(np.diff(y)/np.diff(x)) 
plt.plot(sp.diff(y,order=1, period=400)) 

матчи довольно хорошо - но я не знаю точно рациональны для период/нормализация в подпрограмме fft