1!------------------------------------------------------------------------------- 2! 3! This file is part of the WSPR application, Weak Signal Propogation Reporter 4! 5! File Name: fchisq.f90 6! Description: 7! 8! Copyright (C) 2001-2014 Joseph Taylor, K1JT 9! License: GPL-3+ 10! 11! This program is free software; you can redistribute it and/or modify it under 12! the terms of the GNU General Public License as published by the Free Software 13! Foundation; either version 3 of the License, or (at your option) any later 14! version. 15! 16! This program is distributed in the hope that it will be useful, but WITHOUT 17! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 18! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 19! details. 20! 21! You should have received a copy of the GNU General Public License along with 22! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 23! Street, Fifth Floor, Boston, MA 02110-1301, USA. 24! 25!------------------------------------------------------------------------------- 26real function fchisq(cx,npts,fsample,a,lag1,lag2,ccfmax,dtmax) 27 28 parameter (NMAX=120*375) 29 complex cx(npts) 30 real a(5) 31 complex*16 w1,ws1 32 complex*16 w2,ws2 33 complex*16 w3,ws3 34 complex*16 w4,ws4 35 complex*16 cs1(0:NMAX) 36 complex*16 cs2(0:NMAX) 37 complex*16 cs3(0:NMAX) 38 complex*16 cs4(0:NMAX) 39 complex z1,z2,z3,z4 40 real*8 twopi,baud,p2 41 real ss(5624) 42 save 43 44 twopi=8.d0*atan(1.d0) 45 baud=12000.d0/8192 46 47! Mix and integrate four channels 48 cs1(0)=0. 49 cs2(0)=0. 50 cs3(0)=0. 51 cs4(0)=0. 52 w1=1.0 53 w2=1.0 54 w3=1.0 55 w4=1.0 56 x0=0.5*(npts+1) !Middle sample 57 s=2.0/npts 58 dt=1.0/fsample 59 do i=1,npts 60 x=s*(i-x0) !x runs from -1 to +1 61 if(mod(i,100).eq.1) then 62 p2=1.5*x*x - 0.5 63! p3=2.5*(x**3) - 1.5*x 64! p4=4.375*(x**4) - 3.75*(x**2) + 0.375 65 dphi1=twopi*dt*(a(1) + x*a(2) + p2*a(3) + 1.5*baud) 66 dphi2=twopi*dt*(a(1) + x*a(2) + p2*a(3) + 0.5*baud) 67 dphi3=twopi*dt*(a(1) + x*a(2) + p2*a(3) - 0.5*baud) 68 dphi4=twopi*dt*(a(1) + x*a(2) + p2*a(3) - 1.5*baud) 69 ws1=cmplx(cos(dphi1),sin(dphi1)) 70 ws2=cmplx(cos(dphi2),sin(dphi2)) 71 ws3=cmplx(cos(dphi3),sin(dphi3)) 72 ws4=cmplx(cos(dphi4),sin(dphi4)) 73 endif 74 w1=w1*ws1 75 w2=w2*ws2 76 w3=w3*ws3 77 w4=w4*ws4 78 cs1(i)=cs1(i-1) + w1*cx(i) 79 cs2(i)=cs2(i-1) + w2*cx(i) 80 cs3(i)=cs3(i-1) + w3*cx(i) 81 cs4(i)=cs4(i-1) + w4*cx(i) 82 enddo 83 84! Compute full-symbol powers at 1/16-symbol steps. 85 nsps=nint(fsample/baud) !Samples per symbol 86 ndiv=16 !Steps per symbol 87 nout=ndiv*npts/nsps !Total steps 88 dtstep=1.0/(ndiv*baud) !Time per output step 89 fac=1.e-5 90 91 ss=0. 92 do i=1,nout 93 j=i*nsps/ndiv 94 k=j - nsps 95 ss(i)=0. 96 if(k.ge.1) then 97 z1=cs1(j)-cs1(k) 98 z2=cs2(j)-cs2(k) 99 z3=cs3(j)-cs3(k) 100 z4=cs4(j)-cs4(k) 101 102 p1=real(z1)**2 + aimag(z1)**2 103 p2=real(z2)**2 + aimag(z2)**2 104 p3=real(z3)**2 + aimag(z3)**2 105 p4=real(z4)**2 + aimag(z4)**2 106! ss(i)=fac*(max(p2,p4) - max(p1,p3)) 107 ss(i)=fac*((p2+p4) - (p1+p3)) 108 endif 109 enddo 110 111 ccfmax=0. 112 call ccf2(ss,nout,lag1,lag2,ccf,lagpk) 113 if(ccf.gt.ccfmax) then 114 ccfmax=ccf 115 dtmax=lagpk*dtstep 116 endif 117 118! Reverse sign (and offset!) because we will be minimizing fchisq 119 fchisq=-ccfmax + 100.0 120 121 return 122end function fchisq 123