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

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

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

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?


python numpy scipy derivative differentiation
person TomCho    schedule 13.02.2017    source источник

Ответы (2)


arrow_upward
5
arrow_downward

Функция 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>

участок2

person Warren Weckesser    schedule 13.02.2017
comment
Я наивно думал, что period относится к периоду между одним измерением и другим (следовательно, 1 в моем случае). Но да, это имеет большой смысл, очевидно. Теперь, когда я это понял, я также могу воспроизвести это с помощью numpy, умножив fft на i*2*pi. - person TomCho; 13.02.2017

arrow_upward
0
arrow_downward

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

person f5r5e5d    schedule 13.02.2017