SVN r8748
This commit is contained in:
Jordan Sherer
@@ -1,107 +0,0 @@
|
||||
subroutine genwspr5(msg,msgsent,itone)
|
||||
|
||||
! Encode a WSPR-LF message, producing array itone().
|
||||
|
||||
use crc
|
||||
include 'wsprlf_params.f90'
|
||||
|
||||
character*22 msg,msgsent
|
||||
character*60 cbits
|
||||
integer*1,target :: idat(9)
|
||||
integer*1 msgbits(KK),codeword(ND)
|
||||
logical first
|
||||
integer icw(ND)
|
||||
integer id(NS+ND)
|
||||
integer jd(NS+ND)
|
||||
integer isync(48) !Long sync vector
|
||||
integer ib13(13) !Barker 13 code
|
||||
integer itone(NN)
|
||||
integer*8 n8
|
||||
data ib13/1,1,1,1,1,-1,-1,1,1,-1,1,-1,1/
|
||||
data first/.true./
|
||||
save first,isync
|
||||
|
||||
if(first) then
|
||||
n8=z'cbf089223a51'
|
||||
do i=1,48
|
||||
isync(i)=-1
|
||||
if(iand(n8,1).eq.1) isync(i)=1
|
||||
n8=n8/2
|
||||
enddo
|
||||
first=.false.
|
||||
endif
|
||||
|
||||
idat=0
|
||||
call wqencode(msg,ntype0,idat) !Source encoding
|
||||
id7=idat(7)
|
||||
if(id7.lt.0) id7=id7+256
|
||||
id7=id7/64
|
||||
icrc=crc10(c_loc(idat),9) !Compute the 10-bit CRC
|
||||
idat(8)=icrc/256 !Insert CRC into idat(8:9)
|
||||
idat(9)=iand(icrc,255)
|
||||
call wqdecode(idat,msgsent,itype)
|
||||
|
||||
write(cbits,1004) idat(1:6),id7,icrc
|
||||
1004 format(6b8.8,b2.2,b10.10)
|
||||
read(cbits,1006) msgbits
|
||||
1006 format(60i1)
|
||||
|
||||
! call chkcrc10(msgbits,nbadcrc)
|
||||
! print*,msgsent,itype,crc10_check(c_loc(idat),9),nbadcrc
|
||||
|
||||
call encode300(msgbits,codeword) !Encode the test message
|
||||
icw=2*codeword - 1 !NRZ codeword
|
||||
|
||||
! Message structure:
|
||||
! I channel: R1 48*(S1+D1) S13 48*(D1+S1) R1
|
||||
! Q channel: R1 D109 R1
|
||||
! Generate QPSK with no offset, then shift the y array to get OQPSK.
|
||||
|
||||
! I channel:
|
||||
n=0
|
||||
k=0
|
||||
do j=1,48 !Insert group of 48*(S1+D1)
|
||||
n=n+1
|
||||
id(n)=2*isync(j)
|
||||
k=k+1
|
||||
n=n+1
|
||||
id(n)=icw(k)
|
||||
enddo
|
||||
|
||||
do j=1,13 !Insert Barker 13 code
|
||||
n=n+1
|
||||
id(n)=2*ib13(j)
|
||||
enddo
|
||||
|
||||
do j=1,48 !Insert group of 48*(S1+D1)
|
||||
k=k+1
|
||||
n=n+1
|
||||
id(n)=icw(k)
|
||||
n=n+1
|
||||
id(n)=2*isync(j)
|
||||
enddo
|
||||
|
||||
! Q channel
|
||||
do j=1,204
|
||||
k=k+1
|
||||
n=n+1
|
||||
id(n)=icw(k)
|
||||
enddo
|
||||
|
||||
! Map I and Q to tones.
|
||||
n=0
|
||||
jz=(NS+ND+1)/2
|
||||
do j=1,jz-1
|
||||
jd(2*j-1)=id(j)/abs(id(j))
|
||||
jd(2*j)=id(j+jz)/abs(id(j+jz))
|
||||
enddo
|
||||
jd(NS+ND)=id(jz)/abs(id(jz))
|
||||
itone=0
|
||||
do j=1,jz-1
|
||||
itone(2*j+1)=(jd(2*j)*jd(2*j-1)+1)/2;
|
||||
itone(2*j+2)=-(jd(2*j)*jd(2*j+1)-1)/2;
|
||||
enddo
|
||||
itone(NS+ND+2)=jd(NS+ND) !### Is this correct ??? ###
|
||||
|
||||
return
|
||||
end subroutine genwspr5
|
||||
@@ -0,0 +1,40 @@
|
||||
subroutine extractmessage174(decoded,msgreceived,ncrcflag)
|
||||
use iso_c_binding, only: c_loc,c_size_t
|
||||
use crc
|
||||
use packjt
|
||||
|
||||
character*22 msgreceived
|
||||
character*87 cbits
|
||||
integer*1 decoded(87)
|
||||
integer*1, target:: i1Dec8BitBytes(11)
|
||||
integer*4 i4Dec6BitWords(12)
|
||||
|
||||
! Write decoded bits into cbits: 75-bit message plus 12-bit CRC
|
||||
write(cbits,1000) decoded
|
||||
1000 format(87i1)
|
||||
read(cbits,1001) i1Dec8BitBytes
|
||||
1001 format(11b8)
|
||||
read(cbits,1002) ncrc12 !Received CRC12
|
||||
1002 format(75x,b12)
|
||||
|
||||
i1Dec8BitBytes(10)=iand(i1Dec8BitBytes(10),128+64+32)
|
||||
i1Dec8BitBytes(11)=0
|
||||
icrc12=crc12(c_loc(i1Dec8BitBytes),11) !CRC12 computed from 75 msg bits
|
||||
|
||||
if(ncrc12.eq.icrc12 .or. sum(decoded(57:87)).eq.0) then !### Kludge ###
|
||||
! CRC12 checks out --- unpack 72-bit message
|
||||
do ibyte=1,12
|
||||
itmp=0
|
||||
do ibit=1,6
|
||||
itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*6+ibit))
|
||||
enddo
|
||||
i4Dec6BitWords(ibyte)=itmp
|
||||
enddo
|
||||
call unpackmsg(i4Dec6BitWords,msgreceived,.false.,' ')
|
||||
ncrcflag=1
|
||||
else
|
||||
msgreceived=' '
|
||||
ncrcflag=-1
|
||||
endif
|
||||
return
|
||||
end subroutine extractmessage174
|
||||
@@ -0,0 +1,164 @@
|
||||
program chkfft
|
||||
|
||||
! Tests and times one-dimensional FFTs computed by four2a().
|
||||
! An all-Fortran version of four2a() is available, but the preferred
|
||||
! version uses calls to the FFTW library.
|
||||
|
||||
parameter (NMAX=8*1024*1024) !Maximum FFT length
|
||||
complex a(NMAX),b(NMAX)
|
||||
real ar(NMAX),br(NMAX)
|
||||
real mflops
|
||||
character infile*12,arg*8
|
||||
logical list
|
||||
common/patience/npatience
|
||||
equivalence (a,ar),(b,br)
|
||||
|
||||
nargs=iargc()
|
||||
if(nargs.ne.5) then
|
||||
print*,'Usage: chkfft <nfft | infile> nr nw nc np'
|
||||
print*,' nfft: length of FFT'
|
||||
print*,' nfft=0: do lengths 2^n, n=2^4 to 2^23'
|
||||
print*,' infile: name of file with nfft values, one per line'
|
||||
print*,' nr: 0/1 to not read (or read) wisdom'
|
||||
print*,' nw: 0/1 to not write (or write) wisdom'
|
||||
print*,' nc: 0/1 for real or complex data'
|
||||
print*,' np: 0-4 patience for finding best algorithm'
|
||||
go to 999
|
||||
endif
|
||||
|
||||
list=.false.
|
||||
nfft=-1
|
||||
call getarg(1,infile)
|
||||
open(10,file=infile,status='old',err=1)
|
||||
list=.true. !A valid file name was provided
|
||||
go to 2
|
||||
1 read(infile,*) nfft !Takje first argument to be nfft
|
||||
2 call getarg(2,arg)
|
||||
read(arg,*) nr
|
||||
call getarg(3,arg)
|
||||
read(arg,*) nw
|
||||
call getarg(4,arg)
|
||||
read(arg,*) ncomplex
|
||||
call getarg(5,arg)
|
||||
read(arg,*) npatience
|
||||
|
||||
call sgran()
|
||||
|
||||
if(list) write(*,1000) infile,nr,nw,ncomplex,npatience
|
||||
1000 format(/'infile: ',a12,' nr:',i2,' nw',i2,' nc:',i2,' np:',i2/)
|
||||
if(.not.list) write(*,1002) nfft,nr,nw,ncomplex,npatience
|
||||
1002 format(/'nfft: ',i10,' nr:',i2,' nw',i2,' nc:',i2,' np:',i2/)
|
||||
|
||||
open(12,file='chkfft.out',status='unknown')
|
||||
open(13,file='fftwf_wisdom.dat',status='unknown')
|
||||
|
||||
if(nr.ne.0) then
|
||||
call import_wisdom_from_file(isuccess,13)
|
||||
if(isuccess.eq.0) then
|
||||
write(*,1010)
|
||||
1010 format('Failed to import FFTW wisdom.')
|
||||
go to 999
|
||||
endif
|
||||
endif
|
||||
|
||||
idum=-1 !Set random seed
|
||||
ndim=1 !One-dimensional transforms
|
||||
do i=1,NMAX !Set random data
|
||||
x=gran()
|
||||
y=gran()
|
||||
b(i)=cmplx(x,y) !Generate random data
|
||||
enddo
|
||||
|
||||
iters=1000000
|
||||
if(list .or. (nfft.gt.0)) then
|
||||
n1=1
|
||||
n2=1
|
||||
if(nfft.eq.-1) n2=999999
|
||||
write(*,1020)
|
||||
1020 format(' NFFT Time rms MHz MFlops iters', &
|
||||
' tplan'/61('-'))
|
||||
else
|
||||
n1=4
|
||||
n2=23
|
||||
write(*,1030)
|
||||
1030 format(' n N=2^n Time rms MHz MFlops iters', &
|
||||
' tplan'/63('-'))
|
||||
endif
|
||||
|
||||
do ii=n1,n2 !Test one or more FFT lengths
|
||||
if(list) then
|
||||
read(10,*,end=900) nfft !Read nfft from file
|
||||
else if(n2.gt.n1) then
|
||||
nfft=2**ii !Do powers of 2
|
||||
endif
|
||||
|
||||
iformf=1
|
||||
iformb=1
|
||||
if(ncomplex.eq.0) then
|
||||
iformf=0 !Real-to-complex transform
|
||||
iformb=-1 !Complex-to-real (inverse) transform
|
||||
endif
|
||||
|
||||
if(nfft.gt.NMAX) go to 900
|
||||
a(1:nfft)=b(1:nfft) !Copy test data into a()
|
||||
t0=second()
|
||||
call four2a(a,nfft,ndim,-1,iformf) !Get planning time for forward FFT
|
||||
call four2a(a,nfft,ndim,+1,iformb) !Get planning time for backward FFT
|
||||
t2=second()
|
||||
tplan=t2-t0 !Total planning time for this length
|
||||
|
||||
total=0.
|
||||
do iter=1,iters !Now do many iterations
|
||||
a(1:nfft)=b(1:nfft) !Copy test data into a()
|
||||
|
||||
t0=second()
|
||||
call four2a(a,nfft,ndim,-1,iformf) !Forward FFT
|
||||
call four2a(a,nfft,ndim,+1,iformb) !Backward FFT on same data
|
||||
t1=second()
|
||||
total=total+t1-t0
|
||||
if(total.ge.1.0) go to 40 !Cut iterations short if t>1 s
|
||||
enddo
|
||||
iter=iters
|
||||
|
||||
40 time=0.5*total/iter !Time for one FFT of current length
|
||||
tplan=0.5*tplan-time !Planning time for one FFT
|
||||
if(tplan.lt.0) tplan=0.
|
||||
a(1:nfft)=a(1:nfft)/nfft
|
||||
|
||||
! Compute RMS difference between original array and back-transformed array.
|
||||
sq=0.
|
||||
if(ncomplex.eq.1) then
|
||||
do i=1,nfft
|
||||
sq=sq + real(a(i)-b(i))**2 + imag(a(i)-b(i))**2
|
||||
enddo
|
||||
else
|
||||
do i=1,nfft
|
||||
sq=sq + (ar(i)-br(i))**2
|
||||
enddo
|
||||
endif
|
||||
rms=sqrt(sq/nfft)
|
||||
|
||||
freq=1.e-6*nfft/time
|
||||
mflops=5.0/(1.e6*time/(nfft*log(float(nfft))/log(2.0)))
|
||||
if(n2.eq.1 .or. n2.eq.999999) then
|
||||
write(*,1050) nfft,time,rms,freq,mflops,iter,tplan
|
||||
write(12,1050) nfft,time,rms,freq,mflops,iter,tplan
|
||||
1050 format(i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1)
|
||||
else
|
||||
write(*,1060) ii,nfft,time,rms,freq,mflops,iter,tplan
|
||||
write(12,1060) ii,nfft,time,rms,freq,mflops,iter,tplan
|
||||
1060 format(i2,i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1)
|
||||
endif
|
||||
if(mod(ii,50).eq.0) call four2a(0,-1,0,0,0)
|
||||
enddo
|
||||
|
||||
900 continue
|
||||
if(nw.eq.1) then
|
||||
rewind 13
|
||||
call export_wisdom_to_file(13)
|
||||
! write(*,1070)
|
||||
!1070 format(/'Exported FFTW wisdom')
|
||||
endif
|
||||
|
||||
999 call four2a(0,-1,0,0,0)
|
||||
end program chkfft
|
||||
Reference in New Issue
Block a user