Fix the floating-point exception caused by a c2r FFT when not enough space allocated for data array.

This commit is contained in:
Joe Taylor 2021-04-12 15:21:38 -04:00
parent 47332aa9ac
commit 3f98270681
2 changed files with 8 additions and 7 deletions

View File

@ -19,9 +19,10 @@ subroutine four2a(a,nfft,ndim,isign,iform)
! This version of four2a makes calls to the FFTW library to do the ! This version of four2a makes calls to the FFTW library to do the
! actual computations. ! actual computations.
use fftw3
parameter (NPMAX=2100) !Max numberf of stored plans parameter (NPMAX=2100) !Max numberf of stored plans
parameter (NSMALL=16384) !Max size of "small" FFTs parameter (NSMALL=16384) !Max size of "small" FFTs
complex a(nfft) !Array to be transformed complex a(nfft+1) !Array to be transformed
complex aa(NSMALL) !Local copy of "small" a() complex aa(NSMALL) !Local copy of "small" a()
integer nn(NPMAX),ns(NPMAX),nf(NPMAX) !Params of stored plans integer nn(NPMAX),ns(NPMAX),nf(NPMAX) !Params of stored plans
integer*8 nl(NPMAX),nloc !More params of plans integer*8 nl(NPMAX),nloc !More params of plans
@ -29,7 +30,6 @@ subroutine four2a(a,nfft,ndim,isign,iform)
logical found_plan logical found_plan
data nplan/0/ !Number of stored plans data nplan/0/ !Number of stored plans
common/patience/npatience,nthreads !Patience and threads for FFTW plans common/patience/npatience,nthreads !Patience and threads for FFTW plans
include 'fftw3.f90' !FFTW definitions
save plan,nplan,nn,ns,nf,nl save plan,nplan,nn,ns,nf,nl
if(nfft.lt.0) go to 999 if(nfft.lt.0) go to 999

View File

@ -8,7 +8,7 @@ subroutine q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol,xpol, &
type(hdr) h !Header for the .wav file type(hdr) h !Header for the .wav file
integer*2 iwave(60*12000) integer*2 iwave(60*12000)
complex ca(MAXFFT1),cb(MAXFFT1) !FFTs of raw x,y data complex ca(MAXFFT1),cb(MAXFFT1) !FFTs of raw x,y data
complex cx(0:MAXFFT2-1),cy(0:MAXFFT2-1),cz(0:MAXFFT2-1) complex cx(0:MAXFFT2-1),cy(0:MAXFFT2-1),cz(0:MAXFFT2)
logical xpol,first logical xpol,first
real*8 fcenter real*8 fcenter
character*12 mycall0,hiscall0 character*12 mycall0,hiscall0
@ -87,10 +87,11 @@ subroutine q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol,xpol, &
npol=1 npol=1
if(xpol) npol=4 if(xpol) npol=4
do ipol=1,npol do ipol=1,npol
if(ipol.eq.1) cz=cx if(ipol.eq.1) cz(0:MAXFFT2-1)=cx
if(ipol.eq.2) cz=0.707*(cx+cy) if(ipol.eq.2) cz(0:MAXFFT2-1)=0.707*(cx+cy)
if(ipol.eq.3) cz=cy if(ipol.eq.3) cz(0:MAXFFT2-1)=cy
if(ipol.eq.4) cz=0.707*(cx-cy) if(ipol.eq.4) cz(0:MAXFFT2-1)=0.707*(cx-cy)
cz(MAXFFT2)=0.
call four2a(cz,2*nfft2,1,1,-1) !Transform to time domain (real), fsample=12000 Hz call four2a(cz,2*nfft2,1,1,-1) !Transform to time domain (real), fsample=12000 Hz
do i=0,nfft2-1 do i=0,nfft2-1
j=nfft2-1-i j=nfft2-1-i