FFTW в Фортране результат содержит только нули - PullRequest
Новая
       51

FFTW в Фортране результат содержит только нули

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

Я пытался написать простую программу для выполнения fft на входном массиве 1D с использованием fftw3. Здесь я использую сейсмограмму в качестве входных данных. Выходной массив, однако, будет содержать только нули.

Я знаю, что ввод правильный, так как я пытался сделать fft того же самого входного файла в MATLAB, что дает правильные результаты. Там нет ошибки компиляции. Я использую f95 для компиляции, однако gfortran также давал почти те же результаты. Вот код, который я написал: - 

program fft

    use functions
    implicit none
    include 'fftw3.f90'
    integer nl,row,col
    double precision, allocatable :: data(:,:),time(:),amplitude(:)
    double complex, allocatable :: out(:)
    integer*8 plan


    open(1,file='test-seismogram.xy')
    nl=nlines(1,'test-seismogram.xy')
    allocate(data(nl,2))
    allocate(time(nl))
    allocate(amplitude(nl))
    allocate(out(nl/2+1))
    do row = 1,nl
        read(1,*,end=101) data(row,1),data(row,2)
        amplitude(row)=data(row,2)
    end do
    101 close(1)


    call dfftw_plan_dft_r2c_1d(plan,nl,amplitude,out,FFTW_R2HC,FFTW_PATIENT)
    call dfftw_execute_dft_r2c(plan, amplitude, out)
    call dfftw_destroy_plan(plan)


    do row=1,(nl/2+1)
        print *,out(row)
    end do


    deallocate(data)
    deallocate(amplitude)
    deallocate(time)
    deallocate(out)
end program fft

Функция nlines() - это функция, которая используется для вычисления количества строк в файле, и она работает правильно. Это определено в модуле, называемом функциями.

Эта программа старается следовать примеру на http://www.fftw.org/fftw3_doc/Fortran-Examples.html

Возможно, я просто допускаю очень простую логическую ошибку, но я серьезно не могу понять, что здесь происходит не так. Любые указатели были бы очень полезны.

Это выглядит примерно так: - 

           .
           .
           .
           (0.0000000000000000,0.0000000000000000)
           (0.0000000000000000,0.0000000000000000)
           (0.0000000000000000,0.0000000000000000)
           (0.0000000000000000,0.0000000000000000)
           (0.0000000000000000,0.0000000000000000)
           .
           .
           .

Мои сомнения касаются непосредственно fftw, поскольку на SO есть тег для fftw, поэтому я надеюсь, что этот вопрос не по теме

1 Ответ

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

Как объяснили в комментариях сначала @roygvib и @Ross, подпрограммы плана перезаписывают входные массивы, потому что они пытаются преобразовать много раз с различными параметрами. Я добавлю некоторые практические соображения.

Вы утверждаете, что заботитесь о производительности. Тогда есть две возможности:

  1. Вы выполняете преобразование только один раз , как показано в вашем коде. Тогда нет смысла использовать FFTW_MEASURE. Подпрограмма планирования во много раз медленнее , чем подпрограмма фактического выполнения плана. Используйте FFTW_ESTIMATE, и это будет намного быстрее.

FFTW_MEASURE говорит FFTW найти оптимизированный план на самом деле вычисление нескольких БПФ и измерение времени их выполнения. в зависимости на вашей машине это может занять некоторое время (часто несколько секунд). FFTW_MEASURE - опция планирования по умолчанию.

FFTW_ESTIMATE указывает, что вместо фактических измерений различные алгоритмы, простая эвристика используется для выбора (вероятно, неоптимальный) план быстро . С этим флагом массив ввода / вывода не перезаписывается во время планирования.

http://www.fftw.org/fftw3_doc/Planner-Flags.html

  1. Вы выполняете одно и то же преобразование много раз для разных данных. Затем вы должны выполнить планирование только один раз перед первым преобразованием, а затем повторно использовать план. Просто сначала составьте план, и только потом вы заполняете массив первыми входными данными. Составление плана перед каждой перевозкой сделает программу чрезвычайно медленной.
...