реализация периодического gaussian

Question

Я пытаюсь реализовать периодический гауссовский язык в C, MATLAB или Python.

Что такое правильный способ оценки периодической функции Гаусса, как определено ниже

Я в настоящее время оценивают по формуле ниже, чтобы избежать суммирование по минус до плюс бесконечности:

Заранее спасибо.

Answer 1

Ну, вам не нужно оценивать бесконечную сумму, потому что как только вы доберетесь до (x-kL) >> 2sigma, вы достигнете пределов точности с плавающей запятой.

Итак, вы должны начать с нахождения минимума x - kL (т. Е. Просто установите x = x mod L и k = 0, законно делать, потому что это бесконечная сумма), а затем добавление членов в k = +/- 1, +/- 2, ... до тех пор, пока вы не достигнете пределов с плавающей запятой. Вот пример кода MATLAB, который иллюстрирует идею - я просто взбивал это, поэтому я не могу обещать, что это ошибка, но похоже, что он демонстрирует некоторые из основных ожидаемых действий.

function [result] = Periodic_Gaussian(x, L, sigma) 

    gaussian = @(y) 1/(2*pi*sigma)*exp(-y.^2 ./ 4 ./sigma^2); 

    x = mod(x, L); 

    oldresult = NaN; 
    newresult = gaussian(x); 
    k = 1; 
    while any(newresult ~= oldresult) 
     oldresult = newresult; 
     newresult = oldresult + gaussian(x-k*L) + gaussian(x+k*L); 
     k = k+1; 
    end 
    result = newresult; 

Надеюсь, это полезно!

EDIT: Отсутствует коэффициент 4 в знаменателе аргумента экспоненте и обновляется код, чтобы при необходимости взять вектор x.

Answer 2

function [result] = PeriodicGaussian(x, L, sigma)                                                



    gaussian = @(y, sigma) 1/(2*pi*sigma)*exp(-y.^2 ./ 2 ./sigma^2);                                           


    x0 = mod(x, L)                                                        
    x1 = mod(x, -1 * L)                                                       
    result = gaussian(x0, sigma) + gaussian(x1, sigma);                                               
    correctionIdx = (x0 == 0 & x1 == 0);                                                  
    result(correctionIdx) = 0.5 * result(correctionIdx);                                              

end