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