Как исправить ошибки сдвига и масштабирования в программе дифференцирования на основе БПФ? - PullRequest
0 голосов
/ 07 сентября 2018

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

Я пытаюсь использовать как можно меньше кода, полностью используя функциональность numpy для поддержания чистоты кода при использовании его операций с массивами, которые, как мне сказали, быстрее, чем явные циклы python.

Из строки 106 таблицы Функциональные отношения, одномерные , в разделе о Важных преобразованиях Фурье статьи Википедии о преобразованиях Фурье, я понял, что преобразование Фурье обладает этим свойством, когда дифференцирование в реальном пространстве эквивалентно умножению в частотном пространстве, так что d^nf/dx^n = ifft[F*(i*w)^n], где F - это преобразование Фурье f.

В поисках более надежной ссылки для определения того, что будет w в d^nf/dx^n = ifft[F*(i*w)^n], я нашел документ PDF по математике дифференцирования с использованием ДПФ. В статье автор обсуждает, что существует бесконечно много так называемых тригонометрических интерполяций для данной функции, и показывает способ получения «менее волнистого» решения, которое является уникальным.

Итак, я использовал его определения и написал небольшой код на Python для проведения различий. Вот код, который я придумал (обратите внимание на последние две строки, так как они изменятся)

N = 51
a = 0
b = 2*np.pi
X  = np.linspace(a,b, N)

L = abs(b-a)
TwoPiOverL = 2.0*np.pi/L
freqs = np.fft.fftfreq(N,1./N)

# THESE COME FROM [1] #
D1 = TwoPiOverL*freqs*1j     
D2 = -(TwoPiOverL*freqs)**2
#=========================#

S = np.sin(X)
fftS=np.fft.fft(S)
C = np.cos(X)

# COMPUTING THE DERIVATIVES #
Y1 = np.fft.ifft(D1*fftS).real   
Y2 = np.fft.ifft(D2*fftS).real

На рисунках сгенерированных графиков запись ~d/dx sin означает производную функции синуса, вычисленную согласно умножению с преобразованием Фурье, вместо аналитического выражения.

Тест 1

Первый тест, используя только определения, найденные в PDF:

First Test, using only the definitions found in the PDF

Итак, первый тест был беспорядок. Производные хорошо интерполируются вокруг середины области, но конечные точки сумасшедшие. В самом PDF-документе в начале раздела 6 автор пишет:

Все вышеизложенное относится к спектральному дифференцированию периодических функций y (x). Если функция непериодическая, артефакты появятся в производных из-за неявных разрывов на конечных точках.

Звучит правдоподобно, и это означает, что я где-то напортачу, потому что, если вся статья написана с обсуждением периодических сигналов, и если мой сигнал попадает в непериодический класс, то есть место для улучшения, если я можно использовать периодический сигнал.

Но я использую DFT, что означает, что я применяю числовую процедуру к массиву. Дискретные данные, которые не несут никакой информации о периодичности. И тогда я подумал:

Ну, как я могу сказать своему дифференциатору Фурье, что сигнал периодический? Что это влечет за собой в коде? Где, в коде, я должен измениться, чтобы передать эту информацию в мою программу?

То, что казалось ответом, пришло из стека переполнения. В этом посте о численном дифференцировании с использованием БПФ (да, избыточно, но я туда доберусь), @dkato дает инструкции, которые я интерпретировал как изменение двух последних строк в моем коде с

Y1 = np.fft.ifft(D1*fftS).real   
Y2 = np.fft.ifft(D2*fftS).real

до

Y1 = np.fft.ifft(np.exp(D1)*fftS).real   
Y2 = np.fft.ifft(np.exp(D2)*fftS).real

Итак, я сделал это, а затем проверил, что я получу.

Тест 2

Второй тест с использованием поправок @dkato:

Second Test, using @dkato corrections

Как видно из рисунка выше, при использовании поправок, предоставленных @dkato, покачивание конечной точки исчезает, но оценка производной все еще довольно далека. Я использую 51 точку данных, поэтому я предполагаю, что это не связано с неточностями в расширении. Кажется, что способ вычисления производных неверен.

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

Y2 = np.fft.ifft(-np.exp(D2)*fftS).real

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

The Point, наконец

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

В этом качестве у меня есть два вопроса к сообществу, если вы будете так любезны, чтобы помочь мнес этой проблемой.

1) Если это правда, что я как-то сказал своему дифференциатору Фурье, что я использую периодическую функцию, как именно я это сделал?В чем причина кода, который работает, и кода, который не работает?

и

2) Что является источником этих странных несоответствий масштаба / сдвига междуаналитические и оценочные производные?Как я могу их исправить?

1 Ответ

0 голосов
/ 10 сентября 2018

Давайте сначала ответим на суть вопроса:

Каков источник этих странных несоответствий масштаба / сдвига между аналитическими и расчетными производными? Как я могу их исправить?

Проблема начинается со следующей строки:

X  = np.linspace(a,b, N)

По умолчанию numpy.linspace включает конечную точку, что вводит разрыв в вычисляемое расширение периодической функции, поскольку конечные точки повторяются.

Чтобы проиллюстрировать, как это может быть проблемой, вы можете взглянуть на вычисленную функцию sin(X) для N=4:

sin(0)
sin(2*np.pi/3)
sin(4*np.pi/3)
sin(2*np.pi)

, для которого фаза увеличивается на 2*np.pi/3 при каждом последующем значении. Повторение этих значений, чтобы получить периодическое продление с периодом N=4, даст

sin(0)
sin(2*np.pi/3)
sin(4*np.pi/3)
sin(2*np.pi)
sin(0)         # == sin(2*np.pi)
sin(2*np.pi/3) # == sin(2*np.pi + 2*np.pi/3)
sin(4*np.pi/3) # == sin(2*np.pi + 4*np.pi/3)
sin(2*np.pi)   # == sin(2*np.pi + 2*np.pi)

Где вы могли видеть, что фаза увеличивается на 2*np.pi/3 на каждом шаге, за исключением границы блока, где она остается неизменной. Этот небольшой разрыв значительно усложняет получение хорошего приближения на основе тригонометрической интерполяции. Результирующая интерполяция содержит значительные колебания, которые дополнительно усиливаются операцией деривации.

Чтобы избавиться от проблемы, вы должны исключить конечную точку с помощью (и в противном случае оставить исходное решение без изменений):

X  = np.linspace(a,b, N, endpoint=False)

Что касается решения в посте @ dkato , похоже, что оно реализует только операцию сдвига во времени. Для синусоидального сигнала, который может выглядеть не так уж плохо, но он не имеет ничего общего с дифференцированием.

В заключительной ноте, чтобы ответить на часть

Если это правда, что я как-то сказал своему дифференциатору Фурье, что я использую периодическую функцию, как именно я это сделал?

Да, это правда, и вы сделали это в тот момент, когда начали использовать дискретное преобразование Фурье (и его реализацию FFT).

...