1
votes

For a project in my mater's degree I have to make a program to calculate the Fourier transform of a signal.
I open the data of a .dat file and then I have to use a subroutine to calculate the Fourier transform of my signal and return it in an array. I calculate the real part and the imaginary part of the Fourier transform.

  • "npts" is the number of points in my sine wave.
  • "nbft" is the number of frequency of my Fourier transform.
  • "temps" mean "time" in French.**
  • "F_reel" if the real part of the Fourier transform.
  • "F_im" is the imaginary part.
  • "delta_nu" is the temporal resolution.
  • "freq" is the frequency.
  • "signal" is the sine wave value.

I tried to translate the program so you can understand what I try to do.

Everything seems to be fine but I keep getting the error SIGABRT.
I think the problem is in my subroutine, at the step 2 but I am not able to find how to solve the problem.

Here is my code :

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

I use write(*,*)F_reel(k) to display the value of the real part of the Fourier transform to check if the program behave correctly, and I am expecting it to return me all the value it calculates.

However, I keep getting this message:

*** 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

Obviously I am doing something wrong, but I can't see what it is.
Do you have any idea?

1
Welcome. Be sure to take the Tour and read How to Ask. It is always important to compile your program for debugging. First, to get a meaningful backtrace, you should use the -g flag. Then, to get error checking and other warnings, use flags such as -Wall -fcheck=all for gfortran or -warn -check for Intel. Consult your compiler's manual. - Vladimir F
You subroutine subroutine calcul_fourier is external. It is not in a module (preferred!) nor it is internal to the program (learn contains). It is very easy to make errors this way. All subroutines should be placed into a module to avoid errors. - Vladimir F

1 Answers

2
votes

Your loop

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

accesses indexes 0 and nbft+1, but you arrays are declared from 1 to nbft:

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

That is an error, you cannot access arrays out of bounds.

This error can be found by your compiler. Just enable error checking, use flags such as -Wall -fcheck=all for gfortran or -warn -check for Intel. Consult your compiler's manual.

BTW, computing the Fourier transform from the definition is very inefficient, you should use the Fast Fourier Transform (FFT). There are many libraries you can use.