Fortran 傅立叶变换中的`munmap_chunk(): 无效指针` 和 SIGABRT

Question

对于我硕士期间的一个项目，我必须编写一个程序来计算信号的傅里叶变换。
我打开一个 .dat 文件的数据，然后我必须使用 子例程 来 计算我的信号的傅里叶变换并 return 它在一个数组中。我计算傅里叶变换的实部和虚部。

"npts"是我的正弦波的点数。
"nbft"是我傅立叶变换的频率数
"temps" 在法语中的意思是 "time"。**
"F_reel"如果傅里叶变换的实部。
"F_im"为虚部
"delta_nu" 是时间分辨率。
"freq"是频率。
"signal"为正弦波值

我试着翻译了这个程序，这样你就能理解我想做什么。

一切似乎都很好，但我一直收到 SIGABRT 错误。
我认为问题出在我的子程序中，在第 2 步，但我找不到解决问题的方法。

这是我的代码：

program puissance_sinus1_nom
                    !!!!! MAIN PROGRAM !!!!!

! DECLARATION OF CONSTANTS AND VARIABLES 
implicit none
integer :: i, ios=0, npts, nbft
double precision, dimension(:), allocatable :: freq, puis
double precision, dimension(:), allocatable :: temps, signal

! INSTRUCTIONS
! 1. Open the file
open(25, file='sinus1_nom.dat', status='old', &
         action='read', iostat=ios)

! 2. Lines count
read(25,*)npts        ! The number of points (npts) is indicated
                      ! in the 1st line of the file I open
allocate(temps(npts))
allocate(signal(npts))
write(*,*)" "

do i = 1,npts
   read(25,*)temps(i),signal(i)
enddo


if (mod(npts,2) .eq. 0) then
   nbft = npts/2 +1
else
   nbft = (npts +1)/2
endif

write(*,*)"Number of frequencies :",nbft
write(*,*) " "
write(*,*) " "

! 3. Use of subroutine
allocate(freq(nbft))
allocate(puis(nbft))
call calcul_fourier(npts,signal,temps,nbft,freq,puis)


! Close everything
deallocate(freq)
deallocate(puis)
deallocate(temps)
deallocate(signal)

close(25)
end program puissance_sinus1_nom

                    !!!!! SUBROUTINE !!!!!

subroutine calcul_fourier(npts,signal,temps,nbft,freq,puis)
! DECLARATION OF VARIABLES
implicit none 
integer :: i, k
double precision :: delta_nu, pi=4*atan(1.)
double precision, dimension(nbft) :: F_reel, F_im

integer, intent(in) :: npts, nbft
double precision, dimension(npts), intent(in) :: temps, signal
double precision, dimension(nbft), intent(out) :: freq, puis

! INSTRUCTIONS
! 1. Frequency resolution calculation
delta_nu = (1./temps(npts))*1.d6
write(*,*)"delta_nu = ",delta_nu," micro Hz"

! 2. Calculation of real and imaginary parts of Fourier transform
do k=1, nbft
   freq(k) = (k-1)*delta_nu
enddo

do k=0, nbft+1
   F_reel(k) = 0
   F_im(k) = 0
   do i=1, npts
   F_reel(k) = F_reel(k) + (2./npts)*signal(i)*cos(2*pi*freq(k)*temps(i))
   F_im(k) = F_im(k) + (2./npts)*signal(i)*sin(2*pi*freq(k)*temps(i))
   enddo
   write(*,*)F_reel(k)    ! display to see if it works
enddo

write(*,*) " "
write(*,*) " "
write(*,*) " "

! 3. power calculation (magnitude^2)
do k=1,nbft
   puis(k) = F_reel(k)**2 + F_im(k)**2
enddo

end subroutine calcul_fourier

我使用 write(*,*)F_reel(k) 来显示傅里叶变换实部的值，以检查程序是否正确运行，我希望它 return 我计算出它的所有值。

但是，我不断收到此消息：

*** Error in `./puissance_sinus1_nom.x': munmap_chunk(): invalid pointer: 0x0000000000e07d10 ***

Program received signal SIGABRT: Process abort signal.

Backtrace for this error:
#0  0x7F4E9392D407
#1  0x7F4E9392DA1E
#2  0x7F4E92E4A0DF
#3  0x7F4E92E4A067
#4  0x7F4E92E4B447
#5  0x7F4E92E881B3
#6  0x7F4E92E8D98D
#7  0x4013A6 in calcul_fourier_
#8  0x401C0A in MAIN__ at puissance_sinus1_nom.f90:?
Abandon

显然我做错了什么，但我看不出是什么。
你有什么主意吗？

Answer 1

你的循环

do k=0, nbft+1
   F_reel(k) = 0
   F_im(k) = 0
   do i=1, npts
   F_reel(k) = F_reel(k) + (2./npts)*signal(i)*cos(2*pi*freq(k)*temps(i))
   F_im(k) = F_im(k) + (2./npts)*signal(i)*sin(2*pi*freq(k)*temps(i))
   enddo
   write(*,*)F_reel(k)    ! display to see if it works
enddo

访问索引 0 和 nbft+1，但是你的数组声明为从 1 到 nbft:

double precision, dimension(nbft) :: F_reel, F_im

这是一个错误，您不能越界访问数组。

这个错误可以被你的编译器发现。只需启用错误检查，对 gfortran 使用 -Wall -fcheck=all 或对 Intel 使用 -warn -check 等标志。请查阅您的编译器手册。

顺便说一句，根据定义计算傅里叶变换效率很低，你应该使用快速傅里叶变换（FFT）。您可以使用许多库。

Fortran 傅立叶变换中的`munmap_chunk(): 无效指针` 和 SIGABRT

`munmap_chunk(): invalid pointer` and SIGABRT in Fortran Fourier transform

fortran

fft

sigabrt