Блог

Для ускорения вы определенно хотите повторно использовать планы FFTW, если это возможно. При повторном использовании плана в нескольких вызовах hilbert удалите fftw_destroy_plan(plan) строку внутри hilbert, иначе план не будет действителен для использования в следующем вызове. Вместо этого уничтожьте план в конце программы.

Редактировать: ранее я пропустил, что вы используете два плана: один для прямого преобразования, другой для обратного. В hilbert() , удалите все вызовы fftw_plan_dft_1d , fftw_destroy_plan , и fftw_cleanup ; единственный вызов FFTW должен быть fftw_execute . Вместо этого создайте оба плана заранее один раз в начале программы и уничтожьте их в конце программы.

После нескольких попыток вот успешное hilbert() улучшение FFTW. Переместил два fftw_plan и позволил им инициализироваться в main() , плюс, fftw_destroy_plan они тоже были перемещены.

 #include <iostream>
#include <fftw3.h>
#include <time.h>


using namespace std;

//macros for real and imaginary parts
#define REAL 0
#define IMAG 1
//length of complex array
#define N 8
fftw_complex y[N];



void hilbert(const double* in, fftw_complex* out, fftw_plan plan_forward, fftw_plan plan_backward)
{
    // copy the data to the complex array
    for (int i = 0; i < N;   i) {
        out[i][REAL] = in[i];
        out[i][IMAG] = 0;
    }
    // creat a DFT plan and execute it
    //fftw_plan plan = fftw_plan_dft_1d(N, out, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(plan_forward);
    // destroy a plan to prevent memory leak
    //fftw_destroy_plan(plan_forward);
    int hN = N >> 1; // half of the length (N/2)
    int numRem = hN; // the number of remaining elements
    // multiply the appropriate value by 2
    //(those should multiplied by 1 are left intact because they wouldn't change)
    for (int i = 1; i < hN;   i) {
        out[i][REAL] *= 2;
        out[i][IMAG] *= 2;
    }
    // if the length is even, the number of the remaining elements decrease by 1
    if (N % 2 == 0)
        numRem--;
    else if (N > 1) {
        out[hN][REAL] *= 2;
        out[hN][IMAG] *= 2;
    }
    // set the remaining value to 0
    // (multiplying by 0 gives 0, so we don't care about the multiplicands)
    memset(amp;out[hN   1][REAL], 0, numRem * sizeof(fftw_complex));
    // creat a IDFT plan and execute it
    //plan = fftw_plan_dft_1d(N, out, out, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute(plan_backward);
    // do some cleaning
    //fftw_destroy_plan(plan_backward);
    //fftw_cleanup();
    // scale the IDFT output
    for (int i = 0; i < N;   i) {
        out[i][REAL] /= N;
        out[i][IMAG] /= N;
    }








}
/* Displays complex numbers in the form a  /- bi. */
void displayComplex(fftw_complex* y)
{
    for (int i = 0; i < N; i  ) {
        if (y[i][IMAG] < 0)
            cout << y[i][REAL] << "-" << abs(y[i][IMAG]) << "i" << endl;
        else
            cout << y[i][REAL] << " " << y[i][IMAG] << "i" << endl;
    }
}



/* Test */
int main()
{
    // input array
    double x[N];
    // output array
    //fftw_complex y[N];
    fftw_plan plan_forward = fftw_plan_dft_1d(N, y, y, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan plan_backward = fftw_plan_dft_1d(N, y, y, FFTW_BACKWARD, FFTW_ESTIMATE);
    



   
    // fill the first of some numbers
    for (int i = 0; i < N;   i) {
        x[i] = i   1;  // i.e.{1 2 3 4 5 6 7 8}
    }
    // compute the FFT of x and store the result in y.
    hilbert(x, y, plan_forward, plan_backward);
    
    // display the result
    cout << "Hilbert =" << endl;
    
    displayComplex(y);
    fftw_destroy_plan(plan_forward);
    fftw_destroy_plan(plan_backward);
   
}
 

Предваряется моими лучшими комментариями…

Хорошо… После нескольких неудачных запусков…

В вашей второй версии были некоторые проблемы.

Примечательно, что вы не выделяли out .

Кроме того, на всякий случай я считаю, что вам нужны две fftw_plan структуры, одна для прямого и одна для обратного.

Они должны быть выделены / инициализированы один раз main .

И удалите все вызовы выделения / уничтожения hilbert . Таким образом, его можно вызывать несколько раз, просто изменяя размещенные данные out .

Код очистки может перейти к нижней части main .

Я создал очищенную и рабочую версию вашей второй версии. Я показал разницу между старым и новым кодом с cpp помощью условных обозначений:

 #if 0
// old code ...
#else
// new code ...
#endif
 

Итак, вот оно:

 #include <iostream>
#include <cstring>
#include <fftw3.h>
#include <time.h>

using namespace std;

//macros for real and imaginary parts
#define REAL 0
#define IMAG 1

//length of complex array
#define N 8

fftw_plan plan_fwd;
fftw_plan plan_bwd;
fftw_complex *out;

void
hilbert(const double *in, fftw_complex *out)
{
    // copy the data to the complex array
    for (int i = 0; i < N;   i) {
        out[i][REAL] = in[i];
        out[i][IMAG] = 0;
    }

    // creat a DFT plan and execute it
    // fftw_plan plan = fftw_plan_dft_1d(N, out, out, FFTW_FORWARD, FFTW_ESTIMATE); // This line has been moved to the main function

    fftw_execute(plan_fwd);

    // destroy a plan to prevent memory leak
#if 0
    fftw_destroy_plan(plan_fwd);
#endif
    int hN = N >> 1;                    // half of the length (N/2)
    int numRem = hN;                    // the number of remaining elements

    // multiply the appropriate value by 2
    // (those should multiplied by 1 are left intact because they wouldn't change)
    for (int i = 1; i < hN;   i) {
        out[i][REAL] *= 2;
        out[i][IMAG] *= 2;
    }

    // if the length is even, the number of the remaining elements decrease by 1
    if (N % 2 == 0)
        numRem--;
    else if (N > 1) {
        out[hN][REAL] *= 2;
        out[hN][IMAG] *= 2;
    }

    // set the remaining value to 0
    // (multiplying by 0 gives 0, so we don't care about the multiplicands)
    memset(amp;out[hN   1][REAL], 0, numRem * sizeof(fftw_complex));
    // creat a IDFT plan and execute it
#if 0
    plan_bwd = fftw_plan_dft_1d(N, out, out, FFTW_BACKWARD, FFTW_ESTIMATE);
#endif
    fftw_execute(plan_bwd);
    // do some cleaning
#if 0
    fftw_destroy_plan(plan_bwd);
    fftw_cleanup();
#endif

    // scale the IDFT output
    for (int i = 0; i < N;   i) {
        out[i][REAL] /= N;
        out[i][IMAG] /= N;
    }
}

/* Displays complex numbers in the form a  /- bi. */
void
displayComplex(fftw_complex * y)
{
    for (int i = 0; i < N; i  ) {
        if (y[i][IMAG] < 0)
            cout << y[i][REAL] << "-" << abs(y[i][IMAG]) << "i" << endl;
        else
            cout << y[i][REAL] << " " << y[i][IMAG] << "i" << endl;
    }
}

/* Test */
int
main()
{

#if 1
    out = (fftw_complex *) malloc(sizeof(fftw_complex) * N);
#endif

    plan_fwd = fftw_plan_dft_1d(N, out, out, FFTW_FORWARD, FFTW_ESTIMATE);
    plan_bwd = fftw_plan_dft_1d(N, out, out, FFTW_BACKWARD, FFTW_ESTIMATE);

    // input array
    double x[N];

    // output array
    fftw_complex y[N];

    // fill the first of some numbers
    for (int i = 0; i < N;   i) {
        x[i] = i   1;                   // i.e.{1 2 3 4 5 6 7 8}
    }
    clock_t begin = clock();

    // compute the FFT of x and store the result in y.
    hilbert(x, y);
    // display the result
    cout << "Hilbert =" << endl;
    displayComplex(y);
    clock_t end = clock();
    double time_spent = (double) (end - begin) / CLOCKS_PER_SEC;

    printf("%f", time_spent);

    fftw_destroy_plan(plan_fwd);
    fftw_destroy_plan(plan_bwd);
    fftw_cleanup();
}
 

Вот вывод программы:

 Hilbert =
0.125 0i
0.5 0i
0.75 0i
1 0i
0.625 0i
0 0i
0 0i
0 0i
0.000171