Взяв вторую производную в FFTW3 - PullRequest
Новая
       39

Взяв вторую производную в FFTW3

1 голос
/ 23 июня 2011

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

Однако я хотел бы взять вторую производную от реальной функции.Для простоты я беру sin(2*pi*t) в качестве контрольного примера.Вот соответствующий код, который я использую (функции FFT в библиотеке): 

int main(void)
{
    int i;
    int nyh = (N/2) + 1;

    double result_array[nyh][2];

    double x_k[nyh][2];
    double x_r[N];

    FILE* psit;
    psit=fopen("psitest.txt","w");  

    init();

    fft(x, result_array);  //function in a library, this has been tested
    psi(result_array, x_k);  
    ifft(x_k, x_r);        //function in a library, this has been tested

    for(i=0;i<N;i++)
    {
        fprintf(psit, "%f\n", x_r[i]);
    }


    fclose(psit);

    return 0;
}

void psi(double array[nyh][2], double out[nyh][2])
{
    int i;

    for ( i = 0; i < N/2; i++ )
    {
        out[i][0] = -4.0*pi*pi*i*i*array[i][0];
        out[i][1] = -4.0*pi*pi*i*i*array[i][1];
    }   
    out[N/2][0]=0.0;
    out[N/2][1]=0.0;
}

void init()
{
    int i;

    for(i=0;i<N;i++) 
    {
        x[i] = sin(2.0*pi*i/N);
    }
}

Теперь возникает проблема: этот алгоритм отлично работает для любой функции вида sin( 2*pi*t*K), где K - целое число,но если я возьму в качестве тестовой функции sin(3*pi*t), алгоритм потерпит неудачу.Я не могу увидеть ошибку в моем коде.

Обратите внимание, что, поскольку функция действительна, мне нужно взять только половину значений k.Это не проблема.

Ответы [ 2 ]

1 голос
/ 23 июня 2011

Я предполагаю, что sin(3*pi*t) вводит разрыв, поскольку он не дает целого числа циклов в интервале выборки.Для большинства приложений, связанных с БПФ, вы бы применили оконную функцию для обработки таких разрывов, но очевидно, что это приведет к появлению термина ошибки в вашей производной, и я не уверен, сможете ли вы исправить это.

0 голосов
/ 03 апреля 2012

Я не знаю, исправили ли вы эту проблему ... Но я предполагаю, что основная проблема заключается в том, что sin (3 Pi t) не является периодическим в области [0,1] (sin (0)! = Sin (3 * Пи)).

БПФ не может работать должным образом ...

...