1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
|
! This source code file was last time modified by Igor UA3DJY on January 18th, 2018
! All changes are shown in the patch file coming together with the full JTDX source code.
subroutine downsam9s(npts8,newdat,fpk,c2)
!Downsample from dd() into c2() so as to yield nspsd samples per symbol,
!mixing from fpk down to zero frequency. The downsample factor is 432.
use, intrinsic :: iso_c_binding
use FFTW3
use timer_module, only: timer
use jt65_mod6
include 'constants.f90'
integer(C_SIZE_T) NMAX1
parameter (NMAX1=653184)
parameter (NFFT1=653184,NFFT2=1512)
type(C_PTR) :: plan !Pointers plan for big FFT
logical, intent(inout) :: newdat
real*4, pointer :: x1(:)
complex c1(0:NFFT1/2)
complex c2(0:NFFT2-1)
real s(5000)
logical first
data first/.true./
common/patience/npatience,nthreads
save plan,first,c1,s,x1
df1=12000.0/NFFT1
npts9=8*npts8 ! 597888 and 618624 for dec on 52nd second
if(npts9.gt.NFFT1) npts9=NFFT1 !### Fix! ###
if(first) then
nflags=FFTW_ESTIMATE
if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT
if(npatience.eq.2) nflags=FFTW_MEASURE
if(npatience.eq.3) nflags=FFTW_PATIENT
if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE
! Plan the FFTs just once
!$omp critical(fftw) ! serialize non thread-safe FFTW3 calls
plan=fftwf_alloc_real(NMAX1)
call c_f_pointer(plan,x1,[NMAX1])
x1(0:NMAX1-1) => x1 !remap bounds
call fftwf_plan_with_nthreads(nthreads)
plan=fftwf_plan_dft_r2c_1d(NFFT1,x1,c1,nflags)
!$omp end critical(fftw)
first=.false.
endif
if(newdat) then
x1(0:npts9-1)=dd(1:npts9)
x1(npts9:NFFT1-1)=0. !Zero the rest of x1
call timer('FFTbig9 ',0)
call fftwf_execute_dft_r2c(plan,x1,c1)
call timer('FFTbig9 ',1)
nadd=int(1.0/df1)
s=0.
do i=1,5000
j=int((i-1)/df1)
do n=1,nadd
j=j+1
s(i)=s(i)+real(c1(j))**2 + aimag(c1(j))**2
enddo
enddo
newdat=.false.
endif
nh2=NFFT2/2
nf=nint(fpk)
i0=int(fpk/df1)
nw=100
ia=max(1,nf-nw)
ib=min(5000,nf+nw)
call pctile(s(ia),ib-ia+1,40,avenoise)
fac=sqrt(1.0/avenoise)
do i=0,NFFT2-1
j=i0+i
if(i.gt.nh2) j=j-NFFT2
if(j.lt.0) j=0 ! protection of getting to negative index values
c2(i)=fac*c1(j)
enddo
call four2a(c2,NFFT2,1,1,1) !FFT back to time domain
return
end subroutine downsam9s
|
