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?
-gflag. Then, to get error checking and other warnings, use flags such as-Wall -fcheck=allfor gfortran or-warn -checkfor Intel. Consult your compiler's manual. - Vladimir Fsubroutine calcul_fourieris 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