mirror
/
morse-wip
kopia lustrzana https://github.com/ag1le/morse-wip
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

Tested with SNR 8 dB file x.txt. Got the following 

QWEHROWNOOUMP*VTH*MZCOT123H67EQUICSRJWNOTKMP0<84>TH*AZROT123H/7*<A0>UIYUROWNO8UMPMVTH*AZ*O0113<88>567*9QUICRROFO

The actual text is QUICK BROWN FOX JUMPED OVER THE LAZY FOX 0123456789
master
AG1LE Mauri 2013-09-02 16:01:42 -04:00
rodzic de3bb04f54
commit 891e5ff92d
122 zmienionych plików z 72061 dodań i 36382 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

BIN
LTRHAT_output.ods
Wyświetl plik

Plik binarny nie jest wyświetlany.

3
TEXT.txt
Wyświetl plik

@ -1,3 +0,0 @@
CWA6 DE LAB IAW THE QUICK GREY FOX JUMPED OVER THE LAZY BROWN DOG ON A SUNNY SUMMER AFTERNOON.
THIS IS A TEST. VVV JVXI JGBA GBEY IQNH OPRP CIPU URUC RHIC MUJX SKYQ

BIN
a.out
Wyświetl plik

Plik binarny nie jest wyświetlany.

23
autocr.f
Wyświetl plik

@ -1,23 +0,0 @@
SUBROUTINE AUTOCR (S5, RS)
DIMENSION S5(512), RS(512), S(1000), RS1(500)
DATA S/1000*0./, XN/0./
XN=XN+1
DO 100 I=1, 500
S(I) = S5(I)
RS1(I)=0.
100 CONTINUE
DO 200 I=1, 500
DO 300 K=1, 500
RS1(I)= RS1(I)+S(K+I-1)*S(K)
300 CONTINUE
200 CONTINUE
DO 400 I=1, 500
RS(I)=(RS(I)*(XN-1.)+RS1(I))/XN
400 CONTINUE
RETURN
END

51
bpfdet.f
Wyświetl plik

@ -1,51 +0,0 @@
SUBROUTINE BPFDET(ZIN, Z)
C
C THIS SUBROUTINE IMPLEMENTS THE BANDPASS FILTER AND
C ENVELOPE DETECTOR FUNCTIONS. THE BPF IS A SIMPLE CASCADE
C OF TWO 2-POLE DIGITAL RESONATORS AT A CENTER FREQ OF
C 1000 HZ. THE LPF OF THE ENVELOPE DETECTOR IS A
C THREE - POLE CHEBYSCHEV 100 HZ LPF.
C
C
DIMENSION A ( 4 )
DATA A/ .000030051, 2.9507982, 2.90396345, -0.953135172/
DATA CK1/1.37158/ , CK2/ 0.9409/, CG/ 0.0150/
DATA C1 / 1.2726/ ,C2/ 0.8100/ , C / 0.1900/
C OFF IS TWO 2 - POLE RESONATORS :
W3 = W2
W2 = W1
W1=C1*W2-C2*W3+C*ZIN
X3=X2
X2=X1
X1=CK1*X2-CK2*X3+CG*W1
ZBPF=X1
C ENVELOPE DETECTOR (SQUARE-LAW):
C sQUARE-
XDET=SQRT(ZBPF**2)
C LOW-PASS FILTER-
Y3=Y2
Y2=Y1
Y1=Y0
Y0=XDET*A(1)
Z3=Z2
Z2=Z1
Z1=Z
Z=Y0+3.*(Y1+Y2)+Y3
Z=Z+A(2)*Z1-A(3)*Z2-A(4)*23
200 CONTINUE
RETURN
END

13470
data.txt
Wyświetl plik

Plik diff jest za duży Load Diff

15
displa.f.backup
Wyświetl plik

@ -1,15 +0,0 @@
SUBROUTINE DISPLA (S1,S2,S3,S4)
DIMENSION S1 (512), S2(512), S3(512), S4(512)
CALL ERASE
CALL PLOTR (S1,512,0,XM,400)
CALL PLOTR (S2,512,0,XM,275)
CALL PLOTR (S3,512, 1 ,1.,150)
CALL PLOTR (S4,512,0,XM,40)
CALL VIEW ('1')
READ 100, WAIT
100 FORMAT(A5)
RETURN
END

1426
errors.txt
Wyświetl plik

Plik diff jest za duży Load Diff

0
fort.20
Wyświetl plik

BIN
imax.ods
Wyświetl plik

Plik binarny nie jest wyświetlany.

108
initl.f
Wyświetl plik

@ -1,108 +0,0 @@
SUBROUTINE INITL
DIMENSION IELMST(400),ILAMI(16),ILAMX(6)
DIMENSION ELEMTR(16,6),RTRANS(5,2),ISX(6)
DIMENSION MEMFCN(400,6),LTRMAP(400)
DIMENSION MEMDEL(6,6),MEMPR(6,6),IBLANK(400)
DIMENSION IARRAY(8)
DIMENSION ITEXT(200)
CHARACTER*70 IALPH
COMMON /BLKLAM/ IELMST,ILAMI,ILAMX
COMMON /BLKRAT/ MEMDEL
COMMON /BLKELM/ ELEMTR
COMMON /BLKSPD/ RTRANS,MEMPR
COMMON /BLKMEM/ MEMFCN
COMMON /BLKS/ ISX
COMMON /BLKTRN/LTRMAP,IBLANK, IALPH
COMMON /BLKTXT/ ITEXT
DATA ISX/1,1,0,0,0,0/
C MEMORY FUNCTION
C MEMFCN(400,6) - ACCESSED BY (LAMBDA, K) - LETTER STATE TABLE
C 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
C .^ .~ .w .p -^ -~ -w -p ^. ^- ~. ~- w. w- p. p-
C K=1 DIT, K=2 DAH, K=3 E-SPC, K=4 CHR-SPC, K=5 WRD-SPC, K=6 PAUSE
C . - ^ ~ w p
DATA MEMFCN /
Clmbd:1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
&9, 11,13,15, 9,11,13,15, 9, 0,11, 0,13, 0,15, 0, 384*0,
&10,12,14,16,10,12,14,16, 0,10, 0,12, 0,14, 0,16, 384*0,
&1, 0, 0, 0, 5, 0, 0, 0, 1, 5, 1, 5, 1, 5, 1, 5, 384*0,
&0, 2, 0, 0, 0, 6, 0, 0, 2, 6, 2, 6, 2, 6, 2, 6, 384*0,
&0, 0, 3, 0, 0, 0, 7, 0, 3, 7, 3, 7, 3, 7, 3, 7, 384*0,
&0, 0, 0, 4, 0, 0, 0, 8, 4, 8, 4, 8, 4, 8, 4, 8, 384*0/
DATA IELMST/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,384*0/
DATA ILAMI/ 3,4,5,6,3,4,5,6,1,2,1,2,1,2,1,2/
DATA ILAMX/ 1,1,0,0,0,0/
DATA LTRMAP/3,4,5,6,3,4,5,6,1,2,1,2,1,2,1,2,384*0/
DATA IALPH /'ABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890;:%&'/
C /'A','B','C','D','E','F','G','H','I','J','K','L','M'/
C 2 'N','O','P','Q','R','S','T','U','V','W','X','Y','Z',
C 2 '1','2','3','4','5','6','7'/
C 2 '8','9','0',';',':','%','&',0,0,'^',',','AS','SN',
C 2 0,0,0,0,'NR','NO','GA','OK','AR','SK',0,0,0,0,
C 2 'IMI',0,0,0,0,'BT',0,0,0,'EEE'/
DATA IBLANK/400*0/
DATA ELEMTR/.55,.5,.5,.5,.55,.5,.5,.5,8*0.,
&.45,.5,.5,.5,.45,.5,.5,.5,8*0.,
&8*0.,.581,.54,.923,.923,.923,.923,.95,.95,
&8*0.,.335,.376,.062,.062,.062,.062,.04,.04,
&8*0.,.069,.069,.012,.012,.012,.012,.009,.009,
&8*0.,.015,.015,.003,.003,.003,.003,.001,.001/
DATA RTRANS /.1,.2,.4,.2,.1,.15,.2,.3,.2,.15/
DATA MEMDEL /0,0,2,2,5,10,0,0,2,2,5,10,
&2,2,0,0,0,0,2,2,0,0,0,0,2,2,0,0,0,0,
&2,2,0,0,0,0/
DATA MEMPR /0,0,1,2,1,2,0,0,1,2,1,2,1,1,0,0,0,0,
&1,1,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0/
C PRINT 5, IALPH(1:36)
C5 FORMAT(1X,A36)
OPEN (UNIT=20, FILE='MORSEM')
DO 10 I=1,300
READ (20,30) (IARRAY(K),K=1,8)
30 FORMAT(8I3)
C
C PROCESS ONE ROW OF MORSEM FILE - IARRAY(8)
C LTRMAP(1..16) = 3,4,5,6,3,4,5,6,1,2,1,2,1,2,1,2,0,...
C IELMST(1..16) = LAST TWO ELEMENTS
C 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
C .^ .~ .w .p -^ -~ -w -p ^. ^- ~. ~- w. w- p. p-
DO 11 K=1,6
11 MEMFCN(I,K) = IARRAY(K+2)
LTRMAP(I) = IARRAY(1)
IELMST(I) = IARRAY(2)
IF ((IELMST(I).EQ.7).OR.(IELMST(I).EQ.3)) IBLANK(I) = 1
IF ((IELMST(I).EQ.8).OR.(IELMST(I).EQ.4)) IBLANK(I) = 2
10 CONTINUE
ENDFILE 20
OPEN(UNIT=20,FILE='OUTPUT')
DO 50 I=1,300
WRITE(20,40) (MEMFCN(I,K),K=1,6)
40 FORMAT(10X,6(I3,2X))
50 CONTINUE
ENDFILE 20
OPEN(UNIT=20,FILE='TEXT')
DO 60 I=1,108
READ(20,70) ITEXT(I)
70 FORMAT(I2)
60 CONTINUE
ENDFILE 20
RETURN
END

70
inputl.f
Wyświetl plik

@ -1,70 +0,0 @@
SUBROUTINE INPUTL
DIMENSION ESEP(6),EDEV(6)
COMMON /BLK1/TAU
COMMON /BLK6/ DMEAN,XDUR,ESEP,EDEV
COMMON /BLK2/ WC, WCHIRP, ASIGMA, BSIGMA, PHISGM, RSIGM, TCHIRP,
&GAMM
DATA TAU/.000125/
DATA ESEP/1,3,1,3,7,14/
DATA EDEV/6*0./
DATA XDUR/0./
GO TO 600
PRINT 100
100 FORMAT(1X,'INPUT KEYING PARMS: RATE,MEAN ELEM DURATIONS')
C READ 200,RATE,(ESEP(K), K=1,6)
PRINT 150
150 FORMAT(1X,'INPUT ELEM DURATION STD OEVIATIONS')
C READ 200,(EDEV(K),K=1,6)
200 FORMAT(F7.1)
PRINT 300
300 FORMAT(1X,'INPUT SIG PARMS- AVAR,BVAR,FCHIRP,TCHIRP,PHIVAR')
C READ 210,AVAR,BVAR,FCHIRP,TCHIRP,PHIVAR
210 FORMAT(F7.1,F7.1,F7.1,F7.1,F7.1)
PRINT 400
400 FORMAT(1X,'INPUT SIG PARMS: GAMMA,FREQ,NOISE')
C READ 220, GAMM,FC,RNOISE
220 FORMAT(F7.1,F7.1,F7.1)
600 CONTINUE
C ADDED BY MAURI
AVAR = 4.
BVAR = 4.
RATE = 25.
FC = 4000.
FCHIRP = 20.
RNOISE = 4.
PHIVAR = 4.
TCHIRP = 2.
GAMM = 1.
ASIGMA=SQRT(AVAR)
BSIGMA=SQRT(BVAR)
PHISGM=SQRT(PHIVAR)
RSIGM=SQRT(RNOISE)
DMEAN=1200./RATE
RC=6.28319*FC
WCHIRP = 6.2831*FCHIRP
C IF (ESEP(1).NE.0.) GO TO 500
ESEP(1)=1.
ESEP(2)=3.
ESEP(3)=1.
ESEP(4)=3.
ESEP(5)=7.
ESEP(6)=14.
500 RETURN
END

BIN
ipmax.ods
Wyświetl plik

Plik binarny nie jest wyświetlany.

96
kalfil.f
Wyświetl plik

@ -1,96 +0,0 @@
SUBROUTINE KALFIL (Z, IP, RN, ILX, IXS, KELEM, JNODE, ISRATE, DUR,
&ILRATE, PIN, LKHDJ)
C
C THIS SUBROUTINE COMPUTES THE ARRAY OF KALMAN FILTER
C RECURSIONS USED TO DETERMINE THE LIKELIHOODS.
C
C VARIABLES:
C Z - INPUT MEASUREMENT
C IP - INPUT PATH IDENTITY
C RN - INPUT NOISE POWER ESTIMATE
C ILX - INPUT SAVED KEYSTATE ON PATH IP
C IXS - INPUT KEYSTAT OF NEW NODE
C KELEM - INPUT ELEM STATE OF NEW NODE
C ISRATE INPUT SPEED STATE OF NEW NODE
C DUR - INPUT CURRENT DURATION OF ELEMENT ON IP
C ILRATE INPUT SPEED STATE ON PATH IP
C PIN - TRANS PROB FROM PATH IP TO NODE N
C LKHDJ - OUTPUT CALCULATED LIKELIHOOD VALUE
C
C SUBROUTINES USED
C MODEL - OBTAINS THE SIGNAL-STATE-DEPENDENT LINEAR
C MODEL FOR THE KALMAN FILTER RECURSIONS
C
REAL LKHDJ
REAL A
DIMENSION YKKIP(25), PKKIP(25)
DIMENSION YKKSV(750), PKKSV(750)
COMMON/BLKSV1/YKKIP, PKKIP, YKKSV, PKKSV
DATA YKKIP/25*.5/,PKKIP/25*.10/
DATA PINMIN/.00010/
C IF TRANSITION PROBABILITY IS VERY SMALL, DON'T
C BOTHER WITH LIKELIHOOD CALCULATION:
C
IF (PIN.GT.PINMIN) GO TO 100
LKHDJ = 0.
GO TO 400
C
C OBTAIN STATE-DEPENDENT MODEL PARAMETERS:
C
100 CALL MODEL(DUR, KELEM, ILRATE, ISRATE, IXS, PHI, QA, HZ)
C
C GET PREVIOUS ESTIMATES FOR PATH IP
C
YKK = YKKIP(IP)
PKK = PKKIP(IP)
C
C IMPLEMENT KALMAN FILTER FOR THIS TRANSITION
C
YPRED=PHI*YKK
PPRED=PHI*PKK*PHI+QA
PZ=HZ*PPRED+RN
PZINV=1./PZ
G=PPRED*HZ*PZINV
PEST=(1.-G*HZ)*PPRED
ZR = Z-HZ*YPRED
YKKSV(JNODE)=YPRED+G*ZR
PKKSV(JNODE)=PEST
IF(YKKSV(JNODE) .LE. .01) YKKSV(JNODE)=.01
A=0.5*PZINV*ZR**2
IF(A .LE. 1000.) GO TO 200
LKHDJ=0.
GO TO 400
200 EXPA = EXP(-A)
LKHDJ=(1./SQRT(PZ))*EXP(-A)
GO TO 400
PRINT 3, IP, Z,HZ,QA,PHI,PZ,ZR,G,PEST,YKK,YKKSV(JNODE),A,LKHDJ
3 FORMAT(1X,'KALFIL:', I3,2X,12(F9.3,1X))
400 RETURN
END

116
key.f
Wyświetl plik

@ -1,116 +0,0 @@
SUBROUTINE KEY(DUR,X)
DIMENSION ESEP (6),EDEV(6),MORSE(10,40)
DIMENSION ITEXT(200)
C DIMENSION IOUTP(500)
C CHARACTER(500) IOUTP
CHARACTER*500 IOUTP
C DIMENSION ISYMBL(6)
C CHARACTER*6 ISYMBL
CHARACTER(LEN=6), PARAMETER :: ISYMBL = '.-^CWP'
CHARACTER*1 CHR
COMMON/BLKEND/IEND
COMMON /BLK1/TAU
COMMON /BLK6/ DMEAN,XDUR,ESEP,EDEV
COMMON /BLKTXT/ITEXT
DATA IEND/0/
DATA IK/0001000000000/
DATA LTR/1/,NELM/0/, N/0/, NLTR/1/
DATA MORSE/1,3,2,0,0,0,0,0,0,0,
&2,3,1,3,1,3,1,0,0,0,2,3,1,3,2,3,1,0,0,0,
&2,3,1,3,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,
&1,3,1,3,2,3,1,0,0,0,2,3,2,3,1,0,0,0,0,0,
&1,3,1,3,1,3,1,0,0,0,1,3,1,0,0,0,0,0,0,0,
&1,3,2,3,2,3,2,0,0,0,2,3,1,3,2,0,0,0,0,0,
&1,3,2,3,1,3,1,0,0,0,2,3,2,0,0,0,0,0,0,0,
&2,3,1,0,0,0,0,0,0,0,2,3,2,3,2,0,0,0,0,0,
&1,3,2,3,2,3,1,0,0,0,2,3,2,3,1,3,2,0,0,0,
&1,3,2,3,1,0,0,0,0,0,1,3,1,3,1,0,0,0,0,0,
&2,0,0,0,0,0,0,0,0,0,1,3,1,3,2,0,0,0,0,0,
&1,3,1,3,1,3,2,0,0,0,1,3,2,3,2,0,0,0,0,0,
&2,3,1,3,1,3,2,0,0,0,2,3,1,3,2,3,2,0,0,0,
&1,3,1,3,2,3,2,3,2,0,1,3,1,3,1,3,2,3,2,0,
&1,3,1,3,2,3,2,3,2,0,1,3,1,3,1,3,2,3,2,0,
&1,3,1,3,1,3,1,3,2,0,1,3,1,3,1,3,1,3,1,0,
&2,3,1,3,1,3,1,3,1,0,2,3,2,3,1,3,1,3,1,0,
&2,3,2,3,2,3,1,3,1,0,2,3,2,3,2,3,2,3,1,0,
&2,3,2,3,2,3,2,3,2,0,40*0/
DATA IELM/1/
C REPLACE WITH THIS: CHARACTER(LEN=6), PARAMETER :: ISYMBL = '.-^CWP'
C DATA ISYMBL/'.','-','^','C','W','P'/
BETA=1000.*TAU*DUR
IF(BETA .LT. XDUR) GO TO 200
NELM=NELM+1
IELM=MORSE(NELM,LTR)
IF(IELM.NE. 0) GO TO 100
NELM=0
C Y=RAND(IK)
IELM=4
IF(Y .GT. .9) IELM=5
IF((Y.LE. .9).AND.(Y.GT. .3)) IELM=6
C Y=RAND(IK)
C Y=35*(Y-0.001)+1.
C IY=Y
C LTR=IY+1
LTR = LTR +1
IF (LTR.GE.36) LTR = 1
GO TO 100
NLTR=NLTR+1
LTR=ITEXT(NLTR)
IF(LTR .EQ. 0) IELM=4
IF(LTR .EQ. 37)IELM=5
IF(LTR .EQ. 38)IELM=6
NLTR=NLTR+1
LTR=ITEXT(NLTR)
100 N=N+1
C TO MAKE COMPATIBLE WITH F2C - REPLACE BELOW WITH THIS
IOUTP(N:N+1) = ISYMBL(IELM:IELM+1)
C IOUTP(N) = ISYMBL(IELM)
IF(N .LT. 300) GO TO 110
N=0
NLTR=0
IEND=1
PRINT 90
90 FORMAT(1X,'END OF RUN; INPUT DATA WAS:')
DO 10 K=1,10
K1=(K-1)*50+1
K2=K*50
C PRINT 105,(IOUTP(L), L=K1,K2)
PRINT 105,IOUTP(K1:K2)
105 FORMAT(1X,A50)
10 CONTINUE
READ 105,WAIT
110 XM=ESEP(IELM)*DMEAN
XSIGM=EDEV(IELM)*DMEAN
C Y=RAND(IK)
C Y=2.*(Y-.5)
C XDUR=XM+Y+XSIGM
XDUR = XM + XSIGM
IF(XDUR .LT. 20.) XDUR=20.
C MARK x=1 SPACE x=0
X=1.
IF(IELM .GE. 3) X=0
200 CONTINUE
C PRINT 300, NELM, LTR, IELM, XM, BETA, XDUR
C300 FORMAT('KEY: ',3(I3,2X),3(F8.2,2X))
RETURN
END

62
likhd.f
Wyświetl plik

@ -1,62 +0,0 @@
SUBROUTINE LIKHD(Z,RN,IP,LAMBDA,DUR,ILRATE,P,LKHD)
C
C THIS SUBROUTINE CALCULATES,FOR EACH PATH
C EXTENSION TO STATE N, THE LIKELIHOOD OF THAT
C TRANSITION GIVEN THE MEASUREMENTZ. IT USES
C AN ARRAY OF LINEAR (KALMAN) FILTERS TO DO SO.
C
C VARIABLES:
C Z- INPUT MEASUREMENT
C RN- INPUT NOISE POWER ESTIMATE
C IP- INPUT SAVED PATH IDENTITY
C LAMBDA- INPUT SAVED LTR STATE IDENTITY
C DUR- INPUT SAVED DURATION OF ELEMENT ON PATH IP
C ILRATE- INPUT SAVED DATA RATE (SPEED)
C P- INPUT TRANSITION PROBABILITIES
C LKHD- OUTPUT COMPUTED LIKELIHOODS FOR EACH TRANS
C
C SUBROUTINES USED:
C KALFIL-KALMAN FILTER FOR EACH NEW PATH
C
REAL LKHD,LKHDJ
DIMENSION P(25,30),LKHD(750)
DIMENSION IELMST(400), ILAMI(16),ILAMX(6)
DIMENSION ISX(6)
COMMON/BLKLAM/IELMST,ILAMI,ILAMX
COMMON/BLKS/ISX
C OBTAIN SAVED KEYSTATE:
IF (LAMBDA .EQ. 0) GO TO 200
KELEM=ILAMI(IELMST(LAMBDA))
ILX=ILAMX(KELEM)
C FOR EACH STATE:
DO 100 K=1,6
DO 100 I=1,5
C OBTAIN KEYSTATE, RATE STATE, STATE N, NEW NODE:
IXS = ISX(K)
ISRATE = I
N = (I-1)*6+K
J = (IP-1)*30+N
PIN =P(IP,N)
C COMPUTE AND STORE LIKELIHOOD:
CALL KALFIL(Z,IP,RN,ILX,IXS,KELEM,J,ISRATE,DUR,ILRATE,PIN,LKHDJ)
LKHD(J) = LKHDJ
GO TO 100
IF (PIN.LE. 1.E-06) GO TO 100
PRINT 110,IP,N, Z,LAMBDA,K,I, ILRATE,DUR,PIN,LKHDJ,RN
110 FORMAT(1X,'LIKHD:',I2,1X,I2, 1X, F5.3,2X,I3,3(2X,I2),4(3X,F8.3))
C 2 F8.4,2X,F8.4)
100 CONTINUE
200 RETURN
END

50
model.f
Wyświetl plik

@ -1,50 +0,0 @@
SUBROUTINE MODEL(DUR,IELM,ILR,ISR,IXS,PHI,QA,HZ)
C
C THIS SUBROUTINE COMPUTES THE PARAMETERS OF THE
C OBSERVATION STATE TRANSITION MATRIX PHI, THE
C MEASUREMENT MATRIX, AND THE COVARIANCES.
C
C VARIABLES:
C DUR- INPUT ELEMENT DURATION
C IELM- INPUT ELEMENT TYPE
C ILR- INPUT SAVED RATE
C ISR- INPUT RATE OF NEW STATE
C IXS- INPUT KEYSTATE OF NEW STATE
C
C PHI- OUTPUT STATE TRANSITION MATRIX ENTRY FOR
C SIGNAL AMPLITUDE STATE
C QA- OUTPUT COVARIANCE FOR AMPLITUDE STATE
C HZ- OUTPUT MEASUREMENT MATRIX VALUE
C
C COMPUTE MEASUREMENT COEFFICIENT:
HZ=IXS
C COMPUTE PHI AND AMPLITUDE STATE VARIANCE (Q):
R1=1200./ILR
BAUDS=DUR/R1
IF (BAUDS.GE.14.) BAUDS=14.
IF(IELM.GE.3) GO TO 100
QA=.0001
PHI=1.
GO TO 300
100 IF(IXS.EQ.0) GO TO 200
PHI=1.
QA=0.15*EXP(0.6*(BAUDS-14.))
QA=QA+.01*BAUDS*EXP(.2*(1.-BAUDS))
GO TO 300
200 XSAMP=22.4*R1
PHI=10.**(-2/XSAMP)
IF(BAUDS.GE.14.) PHI=1.
QA=0.
300 RETURN
END

41
morse.f
Wyświetl plik

@ -1,41 +0,0 @@
PROGRAM MORSE
INTEGER ELMHAT, XHAT, IPMAX
REAL X, ZOUT, RN,SPDHAT
DIMENSION S1(512), S2(512), S3(512), S4(512)
DATA RN/.1/
DATA NP/0/
C CALL SRAND(86456)
CALL INITL
CALL INPUTL
PRINT 4
4 FORMAT('MORSE:, X, PMAX, IPMAX, ELMHAT, ZSIG,RN,SPDHAT,ZDET')
1 DO 2 N1=1,512
DO 3 N2=1,18
CALL SIMSGI(X,ZSIG)
CALL RCVR(ZSIG,ZRCV)
CALL BPFDET(ZRCV,ZDET)
NP=NP+1
C DECIMATE 4 kHz by 40 down to 100Hz - 5 ms sample time for PROCES
IF(NP.LT.40) GO TO 3
NP=0
CALL NOISE(ZDET,RN,ZOUT)
C RN = RAND()
RN=0.01
CALL PROCES(X,RN,XHAT,PX,ELMHAT,LTRHAT, SPDHAT, IMAX,PMAX)
3 CONTINUE
N=N1
C PRINT 5, X, PMAX, IMAX, ELMHAT,ZSIG, RN, SPDHAT, ZDET
C5 FORMAT('MORSE:',2(F10.3,2X),2(I7,1X),4(F10.3,2X) )
CALL STATS(ZDET,Z,PX,XHAT,S1,S2,S3,S4,N)
2 CONTINUE
C CALL DISPLA(S1,S2,S3,S4)
GO TO 1
STOP
END

56
noise.f
Wyświetl plik

@ -1,56 +0,0 @@
SUBROUTINE NOISE(ZIN,RN,Z)
C
C THIS SUBROUTINE ESTIMATES THE NOISE POWER IN THE
C ENVELOPE DETECTED OUTPUT FOR USE BY THE KALMAN
C FILTERS. AN ESTIMATE OF THE NOISE POWER IS
C ALSO SUBTRACTED FROM THE ENVELOPE DETECTED SIGNAL
C IN ORDER TO PRODUCE A ZERO-MEAN NOISE PROCESS
C AT THE INPUT TO THE MORSE PROCESSOR.
C
DIMENSION YLONG(200),YSHORT(50)
DATA YLONG/200*1./
DATA YSHORT/50*1./
DATA KL/1/,KKL/1/,KS/1/,KKS/1/
DATA YMIN1/1./,YMIN2/1./,YMAVG/.05/
KL=KL+1
IF(KL .EQ. 201) KL=1
FS=KS+1
IF(KS.EQ.51)KS=1
KKL=KKL+1
IF(KKL.GE.200) KKL=200
KKS=KKS+1
IF (KKS.GE.50) KKS=50
IF (KKS.LE.2) GO TO 10
YLONG (KL)=ZIN
YSHORT (KS)=ZIN
YMIN1=ZIN
YMIN2=ZIN
10 DO 100 I=1,KKL
IF (YLONG (I).GT.YMIN1) GO TO 100
YMIN1=YLONG(I)
100 CONTINUE
DO 200 I=1,KKS
IF (YSHORT(I).GT.YMIN2) GO TO 200
YMIN2=YSHORT(I)
200 CONTINUE
YMIN=YMIN1
IF (YMIN2.LT.YMIN1) YMIN=YMIN2
YMAVG=YMAVG+.004* (YMIN-YMAVG)
RN=0.30*YMAVG
IF (RN.LT.0.005) RN=0.005
Z=1.1* (ZIN-2.4*YMAVG-.05)
RETURN
END

4971
output_c.txt
Wyświetl plik

Plik diff jest za duży Load Diff

4576
output_f.txt
Wyświetl plik

Plik diff jest za duży Load Diff

82
path.f
Wyświetl plik

@ -1,82 +0,0 @@
SUBROUTINE PATH(IP,LAMBDA,DUR,ILRATE,LAMSAV,DURSAV,ILRSAV)
C
C PATH COMPUTES THE LTR STATE, DURATION, AND DATA RATE OF
C EACH NEW PATH EXTENDED TO STATE N
C
C VARIABLES:
C IP- SAVED PATH IDENTITY
C LAMBDA- LTR STATE OF SAVED PATH
C DUR- DURATION OF ELEMENT ON SAVED PATH
C ILRATE- DATA RATE OF ELEMENT ON SAVED PATH
C LAMSAV- NEW LTR STATES FOR EACH PATH EXTENSION
C DURSAV- NEW ELEM DURATIONS FOR EACH PATH EXTENSION
C ILRSAV- NEW DATA RATES FOR EACH PATH EXTENSION
C J- NEW PATH IDENTITY
C THE LETTER TRANSITION TABLE, MEMFCN, IS STORED IN COMMON.
C
DIMENSION LAMSAV(750),DURSAV(750),ILRSAV( 750)
DIMENSION MEMFCN(400,6),IELMST(400),ILAMI(16)
DIMENSION ILAMX(6),ISX(6),MEMDEL(6,6)
COMMON/BLKLAM/IELMST,ILAMI,ILAMX
COMMON/BLKMEM/MEMFCN
COMMON/BLKS/ISX
COMMON/BLKRAT/MEMDEL
C FOR EACH ELEM STATE K, AND EACH SPEED I, COMPUTE:
C
DO 100 K=1,6
DO 100 I=1,5
C STATE IDENTITY N:
C
N=(I-1)*6+K
C NEW PATH IDENTITY:
C
J=(IP-1)*30+N
C NEW LTR STATE:
C
IF(LAMBDA .NE. 0) GO TO 50
LAMSAV(J)=0
GO TO 100
50 LAMSAV(J)=MEMFCN(LAMBDA,K)
IF(LAMSAV(J).EQ. 0) GO TO 100
C NEW DURATION:
C OBTAIN KEYSTATE OF SAVED PATH AND NEW STATE:
ILELM=ILAMI(IELMST(LAMBDA))
IXL=ILAMX(ILELM)
IXS=ISX(K)
C CALCULATE DURATION - ADD SAMPLE DURATION 5 ms FOR EACH VALID PATH
DURSAV(J)=DUR*(1-IXS-IXL+2*IXS*IXL)+5.
C NEW DATA RATE:
ILRSAV(J)=ILRATE+(I-3)*MEMDEL(ILELM,K)
GOTO 100
PRINT 75, J, ILRSAV(J),LAMBDA, DUR, DURSAV(J)
75 FORMAT('PATH:',3(I3,2X), 2(F8.3,2X))
100 CONTINUE
200 RETURN
END

53
probp.f
Wyświetl plik

@ -1,53 +0,0 @@
SUBROUTINE PROBP(P,PIN,ISAVE,LKHD)
C
C PROBP COMPUTES THE POSTERIOR PROBABILITY OF EACH NEW PATH
C
C VARIABLES:
C P- INPUT: SAVED PROBS OF PRIOR PATHS
C OUTPUT: COMPUTED POSTERIOR PROBS OF NEW PATHS
C PIN- INPUT TRANSISTION PROBABILITIES
C LKHD- INPUT LIKELIHOODS OF EACH TRANSTION
C PSUM- NORMALIZING CONSTANT (SUM OF P(J))
C
REAL LKHD
DIMENSION P(750), PIN(25,30), LKHD(750)
DIMENSION PSAV(750)
PMAX=0.
PSUM=0.
C FOR EACH SAVED PATH, EACH TRANSITION:
DO 100 I=1, ISAVE
DO 100 N=1,30
C COMPUTE IDENTITY OF NEW PATH:
J=(I-1)*30+N
C PRODUCT OF PROBS, ADD TO PSUM
PSAV(J) = P(I)*PIN(I,N)*LKHD(J)
PSUM = PSUM+PSAV(J)
IF (PSAV(J).LE. PMAX) GO TO 100
PMAX = PSAV(J)
100 CONTINUE
C NORMALIZE TO GET PROBABILITIES; SAVE:
NI=30*ISAVE
DO 200 J=1,NI
P(J)=PSAV(J)/PSUM
C PRINT 50, J, P(J), PSAV(J), PMAX
C50 FORMAT('PROBP:',I3,3(2x,F8.3))
200 CONTINUE
RETURN
END

193
proces.f
Wyświetl plik

@ -1,193 +0,0 @@
SUBROUTINE VISARL(NAM,K)
DIMENSION K(25)
REAL K
CHARACTER *3 NAM
PRINT 10,NAM,(J,J=1,25)
10 FORMAT(A3,25(I6,1x))
PRINT 20,(K(J),J=1,25)
20 FORMAT(' :',25(F6.1,1x))
RETURN
END
SUBROUTINE VISARR(NAM,K,Z,ISAVE,ELMHAT)
DIMENSION K(25)
INTEGER K, ELMHAT
CHARACTER *3 NAM
C PRINT 10,NAM,(J,J=1,25)
C10 FORMAT(A3,25(I6,1x))
PRINT 20,Z,ELMHAT,(K(J),J=1,ISAVE)
20 FORMAT(F5.1,I3,' :',25(I6,1x))
RETURN
END
SUBROUTINE VIS(NAM,P)
DIMENSION P(25,30)
CHARACTER *3 NAM
PRINT 10,NAM,(J,J=1,30)
10 FORMAT(A3,30(I7,1x))
DO 30 I=1,25
PRINT 20,I,(P(I,J),J=1,30)
20 FORMAT(I2,':',30(F7.2,1x))
30 CONTINUE
RETURN
END
SUBROUTINE VISI(NAM, H)
INTEGER H
DIMENSION H(25,30)
CHARACTER *3 NAM
PRINT 10,NAM,(J,J=1,30)
10 FORMAT(A3,30(I3,1x))
DO 30 I=1,25
PRINT 20,I,(H(I,J),J=1,30)
20 FORMAT(I2,':',30(I3,1x))
30 CONTINUE
RETURN
END
SUBROUTINE PROCES(Z,RN,XHAT,PX,ELMHAT,LTRHAT, SPDHAT, IMAX,PMAX)
C THIS SUBROUTINE IMPLEMENTS THE PROCESSING ALGORITHM
C FOR JOINT DEMODULATION, DECODING, AND TRANSLATION OF
C THE RECEIVED MORSE PROCESS. IT TAKES IN A NEW MEASURE-
C MENT, Z, OF THE DETECTED SIGNAL EVERY 5 MSEC AND PRO-
C DUCES AN ESTIMATE OF THE CURRENT KEYSTATE, ELEMENT
C STATE, AND LETTER OF THE RECEIVED SIGNAL.
C
C VIS
C DEFINITIONS OF VARIABLE NAMES:
C Z- INPUT SAMPLE OF DETECTED SIGNAL
C RN- INPUT NOISE POWER ESTIMATE
C XHAT- OUTPUT ESTIMATE OF KEYSTATE
C ELMHAT- OUTPUT ESTIMATE OF ELEMENT STATE
C LTRHAT- OUTPUT ESTIMATE OF LETTER STATE
C
C
C ISAVE- NO. OF PREVIOUS PATHS SAVED
C IPATH- IDENTITY OF SAVED PATH
C LAMBDA (I) - IDENTITY OF LTR STATE OF SAVED PATH I
C DUR (I)- DURATION OF ELEMENT ON PATH I
C ILRATE(I)- IDENTITY OF DATA RATE ON PATH I
C PIN (I,N)- COMPUTED TRANS PROB FROM PATH I TO STATE N
C LAMSAV(J)- IDENTITY OF LTR STATE AT NEW NODE J
C ILRSAV(J)- IDENTITY OF DATA RATE AT NEW NODE J
C LKHD(J)- LIKELIHOOD VALUE FOR NODE J
C P(J)- COMPUTED POSTERIOR PROB OF PATH ENDING AT NEW NODE J
C PSELEM(K)-COMPUTED POSTERIOR PROB OF ELEM K
C SPDHAT -COND MEAN ESTIMATE OF INSTANT DATA RATE
C PX- POSTERIOR PROB THAT KEYSTATE EQUALS 1
C
C THE FOLLOWING SUBROUTINES ARE UTILIZED:
C TRPROB- COMPUTES TRANSITION PROBABILITIES
C PATH- COMPUTES IDENTITY OF NEW PATHS
C LIKHD- COMPUTES THE LIKELIHOOD OF EACH PATH EXTENSION
C PROBP- COMPUTES POSTERIOR PROBS OF EACH NEW PATH
C SPROB- COMPUTES POSTERIOR PROBS OF EACH STATE
C SAVE- SAVES THE HIGHEST PROB PATHS
C TRELIS- FORMS A TRELIS OF SAVED PATHS
C TRANSL- TRANSLATES THE LETTER ESTIMATE
C
C ALL TABLES OF CONSTANTS ARE STORED IN COMMON.
REAL LKHD
INTEGER ELMHAT,XHAT,PATHSV,SORT, IMAX
DIMENSION LAMBDA(25),DUR(25),ILRATE(25),PIN(25,30)
DIMENSION LAMSAV(750),DURSAV(750), ILRSAV(750)
DIMENSION LKHD(750),P(750),PSELEM(5)
DIMENSION PATHSV(25),SORT(25)
DATA ISAVE/25/
DATA LAMBDA/25*5/
DATA ILRATE/5*10,5*20,5*30,5*40,5*50/
DATA P/750*1./
DATA LAMSAV /750*5/
DATA DUR /25*1000./
DATA DURSAV/750*0./
DATA ILRSAV/750*20/, PATHSV/25*5/
C FOR EACH SAVED PATH, COMPUTE:
C TRANSITION PROBABILITY TO NEW STATE (TRPROB);
C IDENTITY OF EACH NEW PATH EXTENDED (PATH);
C LIKELIHOOD OF EACH STATE EXTENSION (LIKHD);
DO 110 I=1,750
C P(I) = 1.0
110 CONTINUE
DO 100 I=1, ISAVE
IPATH=I
CALL TRPROB(IPATH, LAMBDA(I), DUR(I), ILRATE(I),PIN)
C CALL VISARL('DR1',DUR)
CALL PATH(IPATH,LAMBDA(I), DUR(I), ILRATE(I), LAMSAV, DURSAV,
&ILRSAV)
C CALL VIS('DRS', DURSAV)
C CALL VISARR('LBD',LAMBDA)
C CALL VIS('PIN', PIN)
CALL LIKHD(Z, RN, IPATH, LAMBDA(I), DUR(I), ILRATE(I), PIN, LKHD)
C CALL VIS('LKH',LKHD)
100 CONTINUE
C HAVING OBTAINED ALL NEW PATHS, COMPUTE:
C POSTERIOR PROBABILITY OF EACH NEW PATH(PROBP);
C POSTERIOR PROBABILITY OF KEYSTATE, ELEM STATE,
C CONDITIONAL MEAN ESTIMATE OF SPEED(SPROB);
CALL PROBP(P, PIN, ISAVE, LKHD)
CALL SPROB(P, ISAVE, ILRSAV, PELM, ELMHAT, SPDHAT, PX)
C CALL VISI('ILR', ILRSAV)
C CALL VIS('P1:',P)
C PRINT 111, Z, PX, SPDHAT, ELMHAT
C111 FORMAT('PROCES:',3(3X,F15.9), 3X,I3)
XHAT=0
IF(PX.GT. 0.5) XHAT=1
C SAVE THE PATHS WITH HIGHEST PROBABILITY, AND
C STORE THE VALUES CORRESPONDING TO THESE PATHS:
C CALL VISI('ILS',ILRSAV)
CALL SAVEP(P, PATHSV, ISAVE, IMAX, LAMSAV, DURSAV, ILRSAV,
&LAMBDA, DUR, ILRATE, SORT,PMAX)
C CALL VISARR('LBD',LAMBDA,Z,ISAVE,ELMHAT)
C CALL VIS('P2:',P)
C CALL VISARL('DR2',DUR)
C CALL VISARR('ILR',ILRATE)
C GOTO 900
C
C
C DO 900 INL=1, ISAVE
C PRINT 112, INL, P(INL), PATHSV(INL), LAMBDA(INL), DUR(INL),
C &ILRATE(INL), LKHD(SORT(INL))
C112 FORMAT(1X, I3, 2X, F10.7, 2X, I3, 2X, I3, 2X,F6.1, 2X, I3, 2X,
C &F10.5)
C 2 F10.7)
900 CONTINUE
C UPDATE TRELLIS WITH NEW SAVED NODES, AND
C OBTAIN LETTER STATE ESTIMATE:
CALL TRELIS(ISAVE, PATHSV, LAMBDA, IMAX,IPMAX)
200 RETURN
END

52
ptrans.f
Wyświetl plik

@ -1,52 +0,0 @@
SUBROUTINE PTRANS (KELEM, IRATE, LAMBDA, ILRATE, PTRX, PSUM,PIN,N)
C THIS FUNCTION SUBROUTINE RETURNS THE PATH CONDITIONAL TRANSITION
C PROBABILITIES TO EACH ALLOWABLE STATE N.
C VARIABLES:
C KELEM- INPUT CURRENT ELEMENT STATE
C IRATE- INPUT CURRENT DATA RATE STATE
C LAMBDA- INPUT IDENTITY OF CURRENT LTR STATE
C PTRX- INPUT KEYSTATE TRANSITION PROBABILITY
C ELEMTR- ELEMENT TRANSITION PROBABILITY MATRIX
C
C FUNCTION SUBROTINE USED:
C SPDTR- RETURNS DATA RATE TRANSITION PROBS,CONDITIONED ON CURRENT SPACE TYPE.
C
DIMENSION IELMST(400),ILAMI(16),ELEMTR(16,6)
DIMENSION ILAMX(6),PIN(30)
COMMON /BLKLAM/ IELMST,ILAMI,ILAMX
COMMON /BLKELM/ ELEMTR
C IF THE SAVED ELEMENT AND THE ELEMENT OF THE STATE
C N TO WHICH THE PATH IS BEING EXTENDED ARE THE
C SAME, THEN THE STATE TRANS PROB IS SIMPLY
C KEYSTATE TRANS PROB:
C
IF(KELEM .NE. ILAMI(IELMST(LAMBDA))) GO TO 100
PIN(N)=PTRX
C HOWEVER, IF CURRENT DATA RATE STATE = 3, THEN TRANS PROB = 0 ... WHY ?
IF(IRATE .NE. 3) PIN(N)=0.
GO TO 200
C OTHERWISE:
C OBTAIN ELEM TRANS PROBS TABLE:
C
100 PELEM=ELEMTR(IELMST(LAMBDA),KELEM)
C NEXT COMPUTE ELEM-CONDITIONAL SPEED TRANS PROB:
PRATE=SPDTR(IRATE,ILRATE,KELEM,ILAMI(IELMST(LAMBDA)))
C TRANS IS THE PRODUCT:
PIN(N)=(1.-PTRX)*PELEM*PRATE
200 PSUM=PSUM+PIN(N)
C PRINT 50, N, KELEM, IRATE, ILAMI(IELMST(LAMBDA)), PIN(N), PTRX, PELEM, PRATE, PSUM
C50 FORMAT(1X,'PTRANS:',4(I3,2X),5(F8.3,2X))
RETURN
END

0
puppu.txt
Wyświetl plik

4154
puppu2.txt
Wyświetl plik

Plik diff jest za duży Load Diff

26
rcvr.f
Wyświetl plik

@ -1,26 +0,0 @@
SUBROUTINE RCVR (ZIN, ZOUT)
C
C THIS SUBROUTINE CONVERTS THE INPUT SIGNAL AT
C RADIAN FREQ WC TO 1000 Hz.
C
C
COMMON/BLK1/TAU/BLK2/WC
DATA THETA/0./,THETLO/0./
THETA = THETA + WC*TAU
THETA = AMOD(THETA,6.28319)
ZI = ZIN*COS(THETA)
ZQ = ZIN*SIN(THETA)
ZILP = ZILP+ .070* (ZI-ZILP)
ZQLP = ZQLP+ .070* (ZQ-ZQLP)
THETLO = THETLO+6283.2* TAU
THETLO = AMOD(THETLO, 6.28319)
ZOUT = ZILP*COS(THETLO) + ZQLP* SIN(THETLO)
RETURN
END

178
savep.f
Wyświetl plik

@ -1,178 +0,0 @@
SUBROUTINE SAVEP(P,PATHSV,ISAVE,IMAX,LAMSAV,DURSAV,ILRSAV,LAMBDA,
& DUR,ILRATE,SORT,PMAX)
C
C THIS SUBROUTINE PERFORMS THE ALGORITM TO SAVE
C THE PATHS WITH HIGHEST POSTERIOR PROBABILITY.
C IT WILL SAVE A MINIMUM OF 7 PATHS (ONE FOR EACH
C ELEMENT STATE AND ONE ADDITIONAL NODE), AND
C A MAXIMUM OF 25 PATHS. WITHIN THESE LIMITS, IT
C SAVED ONLY ENOUGH TO MAKE THE TOTAL SAVED PROBABILITY
C EQUAL TO POPT.
C
C
C ADDITIONALLY, IT RESORTS THE LAMBDA,DUR,AND ILRATE
C ARRAYS TO CORRESPOND TO THE SAVED NODES.
C
C
C VARIABLES:
C P- INPUT PROBABILITY ARRAY OF NEW NODES
C PATHSV- OUTPUT ARRAY OF THE PREVIOUS NODES TO
C WHICH THE SAVED NODES ARE CONNECTED.
C ISAVE- INPUT: NO. OF PREVIOUS NODES SAVED
C
C OUPUT: NO. OF NODES SAVED AT CURRENT STAGE
C IMAX- INDEX OF HIGHEST PROBABILITY NODE
C LAMSAV- INPUT ARRAY OF LTR STATES AT EACH NEW NODE
C DURSAV- INPUT ARRAY OF SAVED DURATIONS
C ILRSAV- INPUT ARRAY OF SAVED RATES
C LAMBDA- OUTPUT ARRAY OF SAVED LTR STATES, SORTED
C
C ACCORDING TO PROBABILITY
C DUR- OUTPUT ARRAY OF SORTED DURATIONS
C ILRATE- OUTPUT ARRAY OF SORTED RATES
C
C
INTEGER PATHSV,SORT
DIMENSION P(750),PATHSV(25),PSAV(25),SORT(25)
DIMENSION LAMSAV(750),DURSAV(750),ILRSAV(750)
DIMENSION LAMBDA (25), DUR (25), ILRATE(25)
DIMENSION YKKIP(25), PKKIP(25)
DIMENSION YKKSV(750), PKKSV(750)
DIMENSION ICONV(25)
COMMON/BLKSV1/ YKKIP,PKKIP,YKKSV,PKKSV
DATA POPT/.90/
NSAV=0
PSUM=0.
C SELECT SIX HIGHEST PROB ELEMENT STATE NODES:
DO 200 K=1,6
PMAX=0.
DO 100 IP=1, ISAVE
DO 100 I=1,5
J=(IP-1)*30+(I-1)*6+K
IF (P(J) .LT. PMAX) GO TO 100
PMAX=P(J)
JSAV=J
IPSAV=IP
100 CONTINUE
C IF FOLLOWING TWO LINES COMMENTED OUT WE GET NSAV = 7
IF (PMAX .GE. 0.000001) GO TO 150
GO TO 200
150 NSAV=NSAV+1
PSUM=PSUM+PMAX
PSAV(NSAV)=PMAX
PATHSV(NSAV)=IPSAV
SORT(NSAV)= JSAV
200 CONTINUE
C SELECT ENOUGH ADDITIONAL NODES TO MAKE TOTAL
C PROBABILITY SAVED EQUAL TO POPT, OR A MAX OF 25:
520 PMAX=0.
DO 500 IP=1, ISAVE
DO 500 N=1, 30
J=(IP-1)*30+N
DO 510 I=1, NSAV
IF (J .EQ. SORT(I)) GO TO 500
510 CONTINUE
IF (P(J).LE.PMAX) GO TO 500
PMAX = P(J)
JSAV=J
IPSAV=IP
500 CONTINUE
PSUM=PSUM+PMAX
NSAV=NSAV+1
PSAV(NSAV)=PMAX
PATHSV(NSAV)=IPSAV
SORT(NSAV)=JSAV
IF(PSUM.GE.POPT) GO TO 600
IF(NSAV.GE.25) GO TO 600
GO TO 520
C NEW ISAVE EQUALS NO. OF NODES SAVED:
600 ISAVE=NSAV
C SORT THE SAVED ARRAYS TO OBTAIN THE ARRAYS
C TO BE USED FOR THE NEXT ITERATION:
DO 700 I=1,ISAVE
P(I)=PSAV(I)/PSUM
LAMBDA(I)=LAMSAV(SORT(I))
DUR(I)=DURSAV(SORT(I))
ILRATE(I)=ILRSAV(SORT(I))
YKKIP(I)=YKKSV(SORT(I))
PKKIP(I)=PKKSV(SORT(I))
C PRINT 650, I, ILRATE(I), SORT(I), LAMBDA(I), P(I), DUR(I),
C &YKKIP(I), PKKIP(I)
C650 FORMAT('SAVEP:',4(I4,1X), 4(F10.4,1X))
700 CONTINUE
DO 790 I=1,ISAVE
ICONV(I)=1
790 CONTINUE
ISAVM1=ISAVE-1
DO 800 N=1,ISAVM1
IF(ICONV(N).EQ.0) GO TO 800
NPLUS1=N+1
DO 810 K=NPLUS1,ISAVE
IF(ICONV(K).EQ.0) GO TO 810
IF(ILRATE(K).NE.ILRATE(N)) GO TO 810
IF(DUR(K).NE.DUR(N)) GO TO 810
IF(LAMBDA(K).NE.LAMBDA(N)) GO TO 810
ICONV(K)=0
810 CONTINUE
800 CONTINUE
PSUM=0.
N=1
DO 900 I=2,ISAVE
IF (ICONV(I).EQ.0) GO TO 900
N=N+1
LAMBDA(N)=LAMBDA(I)
DUR(N)=DUR(I)
ILRATE(N)=ILRATE(I)
YKKIP(N)=YKKIP(I)
PKKIP(N)=PKKIP(I)
PATHSV(N)=PATHSV(I)
SORT(N)=SORT(I)
P(N)=P(I)
PSUM=PSUM+P(N)
900 CONTINUE
C ALSO OBTAIN HIGHEST PROBABILITY NODE:
ISAVE=N
PMAX=0.
DO 950 I=1,ISAVE
P(I)=P(I)/PSUM
IF(P(I).LE.PMAX) GO TO 950
PMAX=P(I)
IMAX=I
950 CONTINUE
RETURN
END

47
simsgi.f
Wyświetl plik

@ -1,47 +0,0 @@
SUBROUTINE SIMSGI(X,SIG)
COMMON/BLK1/TAU
COMMON/BLK2/WC,WCHIRP,ASIGMA,BSIGMA,PHISMG, RSIGM,TCHIRP,GAMM
DATA XLAST/1./,BETA/1./
DATA AMP/1./,BFADE/0./,THETA/1./,PHI/1./
DUR=BETA
CALL KEY(DUR,X)
BETA=BETA*(1.-X-XLAST+2.*X*XLAST)+1.
TK=X*(1.-XLAST)
XLAST=X
C CALL RANDN(W,1,0.,ASIGMA)
AMP=AMP+TK*W
IF(AMP .LT. 0.01) AMP=.01
C CALL RANDN(W,1,0.,BSIGMA)
BFADE=GAMM*BFADE+W
AMPB=AMP+BFADE
IF(AMPB .LT. 0.001) BFADE=0.001-AMP
AMPB=AMP+BFADE
TDUR=1000.*TAU*BETA
WCHRP=X*WCHIRP*EXP(-TDUR/TCHIRP)
THETA=THETA+(NC+WCHRP)*TAU
THETA=AMOD(THETA,6.28319)
C CALL RANDN(W,1,0.,PHISGM)
PHI=PHI+TK*w
PHI=AMOD(PHI,6.28319)
SIG=X*AMPB*SIN(THETA+PHI)
C CALL RANDN(ZN,1,0.,RSIGM)
C PRINT 100,DUR, X
C100 FORMAT('SIMSG:', 2(2X,F10.3))
RETURN
END

45
spdtr.f
Wyświetl plik

@ -1,45 +0,0 @@
FUNCTION SPDTR(ISRT,ILRT,ISELM,ILELM)
C
C THIS FUNCTION RETURNS THE DATA RATE (SPEED) TRANSITION
C PROBABILITY BASED ON THE CURRENT ELEM TYPE. THE ALLOWABLE
C TRANSITION PROBS ARE STORED IN THE TABLE RTRANS.
C
C VARIABLES:
C ISRT - DATA RATE IDENTITY FOR STATE TO WHICH PATH
C IS BEING EXTENDED
C ILRT - DATA RATE ON CURRENT PATH
C ISELM - ELEM TYPE FOR NEXT STATE
C ILELM - ELEM TYPE ON CURRENT PATH
C
DIMENSION RTRANS(5,2), MEMPR(6,6), MEMDEL(6,6)
COMMON /BLKSPD/RTRANS,MEMPR
COMMON /BLKRAT/MEMDEL
C
C IF SAVED ELEMENT AND NEW ELEMENT ARE THE SAME
C THEN THERE CAN BE NO SPEED CHANGE:
C
IF (ILELM .NE. ISELM) GO TO 100
SPDTR = 1.
IF (ISRT .NE. 3) SPDTR = 0.
GO TO 300
C
C OTHERWISE, OBTAIN SPEED TRANSITION PROB
C
100 IDEL=MEMDEL(ILELM,ISELM)
IND1=MEMPR(ILELM,ISELM)
IF (IND1 .NE. 0) GOTO 200
SPDTR = 0.
GO TO 300
200 IDELSP = (ISRT-3)*IDEL
SPDTR = RTRANS(ISRT,IND1)
ISRATE=ILRT+IDELSP
IF ( ISRATE .GT. 60) SPDTR = 0.
IF ( ISRATE .LT. 10) SPDTR = 0.
300 RETURN
END

57
sprob.f
Wyświetl plik

@ -1,57 +0,0 @@
SUBROUTINE SPROB (P,ISAVE,ILRSAV,PELM,KHAT,SPDHAT,PX)
C
C SPROB COMPUTES THE POSTERIOR PROBS OF THE ELEMENT
C STATES, DATA RATE STATES, AND KEYSTATES BY SUMMING
C OVER THE APPROPRIETE PATHS.
C
C VARIABLE:
C P- INPUT PATH PROBABILITIES
C ISAVE- NUMBER OF PATHS SAVED
C PSELEM- OUTPUT ELEMENT PROB
C SPDHAT- OUTPUT SPEED ESTIMATE (DATA RATE WPM)
C PX- OUTPUT KEYSTATE PROBABILITY
C
DIMENSION P(750),PSELEM(6),ILRSAV(750)
C INITIALIZE:
SPDHAT=0.
PX=0.
C FOR EACH STATE EXTENSION OF PATH M:
C OBTAIN ELEMENT STATE PROBS,KEYSTATE PROBS,SPEED EST:
DO 100 K=1,6
PSELEM(K)=0.
DO 100 I=1,5
N=(I-1)*6+K
DO 100 M=1,ISAVE
J=(M-1)*30+N
PSELEM(K)=PSELEM(K)+P(J)
SPDHAT=SPDHAT+ILRSAV(J)*P(J)
IF(K.GT.2) GO TO 100
PX=PX+P(J)
100 CONTINUE
PELM=0.
DO 200 K=1,6
IF(PSELEM(K).LT.PELM) GO TO 200
PELM=PSELEM(K)
KHAT=K
200 CONTINUE
RETURN
END

0
MORSEM.backup → src/MORSEM
Wyświetl plik

300
src/OUTPUT 100644
Wyświetl plik

@ -0,0 +1,300 @@
9 10 1 0 0 0
11 12 0 2 0 0
13 14 0 0 3 0
15 16 0 0 0 4
9 10 5 0 0 0
11 12 0 6 0 0
13 14 0 0 7 0
15 16 0 0 0 8
9 0 1 2 3 4
0 10 5 6 7 8
11 0 1 2 3 4
0 12 5 6 7 8
13 0 1 2 3 4
0 14 5 6 7 8
15 0 1 2 3 4
0 16 5 6 7 8
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0

108
src/TEXT 100644
Wyświetl plik

@ -0,0 +1,108 @@
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

0
autocr.c → src/autocr.c
Wyświetl plik

0
bpfdet.c → src/bpfdet.c
Wyświetl plik

0
common.h → src/common.h
Wyświetl plik

0
gnuplot_i.h → src/gnuplot_i.h
Wyświetl plik

0
initl.c → src/initl.c
Wyświetl plik

0
inputl.c → src/inputl.c
Wyświetl plik

0
kalfil.c → src/kalfil.c
Wyświetl plik

0
key.c → src/key.c
Wyświetl plik

0
likhd.c → src/likhd.c
Wyświetl plik

0
model.c → src/model.c
Wyświetl plik

32
morse.c → src/morse.c
Wyświetl plik

@ -12,6 +12,8 @@
#include "f2c.h"
#include <stdio.h>
#include <math.h>
#include "morse.h"
/* Main program */ int MAIN__(void)
{
@ -20,7 +22,7 @@
static real rn = .1f;
static integer np = 0;
/* Subroutine */ int s_stop(char *, ftnlen);
int s_stop(char *, ftnlen);
/* Local variables */
static integer n;
@ -29,34 +31,37 @@
static real s1[512], s2[512], s3[512], s4[512], px;
static integer imax, xhat;
static real pmax, zdet, zsig;
extern /* Subroutine */ int rcvr_(real *, real *);
static real zrcv, zout;
extern /* Subroutine */ int noise_(real *, real *, real *), initl_(void),
stats_(real *, real *, real *, integer *, real *, real *, real *,
real *, integer *), bpfdet_(real *, real *);
static integer elmhat;
static real spdhat;
extern /* Subroutine */ int simsgi_(real *, real *), proces_(real *, real
*, integer *, real *, integer *, integer *, real *, integer *,
real *);
static integer ltrhat;
extern /* Subroutine */ int inputl_(void);
int res;
FILE *fp;
/* CALL SRAND(86456) */
/* INITIALIZE DATA STRUCTURES */
initl_();
inputl_();
fp = fopen("x.txt","r");
L1:
for (n1 = 1; n1 <= 512; ++n1) {
for (n2 = 1; n2 <= 18; ++n2) {
simsgi_(&x, &zsig);
res = fscanf(fp,"%f",&x);
if ( res != 1) {
fclose(fp);
fp = fopen("x.txt","r");
x = 0.0;
}
rcvr_(&zsig, &zrcv);
bpfdet_(&zrcv, &zdet);
++np;
/* DECIMATE 4 kHz by 40 down to 100Hz - 5 ms sample time for PROCES */
if (np < 40) {
if (np < 20) {
goto L3;
}
np = 0;
@ -64,6 +69,9 @@ L1:
/* RN = RAND() */
rn = .01f;
proces_(&x, &rn, &xhat, &px, &elmhat, &ltrhat, &spdhat, &imax, &pmax);
// printf("\n %f %d %f %d %f ",x,xhat,px, elmhat, spdhat);
L3:
;
}

11
src/morse.h 100644
Wyświetl plik

@ -0,0 +1,11 @@
extern int rcvr_(real *, real *);
extern int noise_(real *, real *, real *);
extern int initl_(void);
extern int stats_(real *, real *, real *, integer *, real *, real *, real *, real *, integer *);
extern int bpfdet_(real *, real *);
extern int simsgi_(real *, real *);
extern int proces_(real *, real *, integer *, real *, integer *, integer *, real *, integer *, real *);
extern int inputl_(void);

0
noise.c → src/noise.c
Wyświetl plik

0
path.c → src/path.c
Wyświetl plik

0
probp.c → src/probp.c
Wyświetl plik

0
proces.c → src/proces.c
Wyświetl plik

0
ptrans.c → src/ptrans.c
Wyświetl plik

0
rcvr.c → src/rcvr.c
Wyświetl plik

0
savep.c → src/savep.c
Wyświetl plik

0
simsgi.c → src/simsgi.c
Wyświetl plik

0
spdtr.c → src/spdtr.c
Wyświetl plik

0
sprob.c → src/sprob.c
Wyświetl plik

0
stats.c → src/stats.c
Wyświetl plik

0
transl.c → src/transl.c
Wyświetl plik

0
trelis.c → src/trelis.c
Wyświetl plik

0
trprob.c → src/trprob.c
Wyświetl plik

71622
src/x.txt 100644
Wyświetl plik

Plik diff jest za duży Load Diff

0
xtrans.c → src/xtrans.c
Wyświetl plik

8
stats.f
Wyświetl plik

@ -1,8 +0,0 @@
SUBROUTINE STATS (XIN1, XIN2, XIN3, XIN4, S1, S2, S3, S4, N)
DIMENSION S1 (512), S2(512), S3(512), S4(512)
S1(N) = XIN1
S2(N) = XIN2
S3(N) = XIN3
S4(N) = XIN4
RETURN
END

111
tmpcsave/autocr.c
Wyświetl plik

@ -1,111 +0,0 @@
/* autocr.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Subroutine */ int autocr_(real *s5, real *rs)
{
/* Initialized data */
static real s[1000] = { 0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f };
static real xn = 0.f;
static integer i__, k;
static real rs1[500];
/* Parameter adjustments */
--rs;
--s5;
/* Function Body */
xn += 1;
for (i__ = 1; i__ <= 500; ++i__) {
s[i__ - 1] = s5[i__];
rs1[i__ - 1] = 0.f;
/* L100: */
}
for (i__ = 1; i__ <= 500; ++i__) {
for (k = 1; k <= 500; ++k) {
rs1[i__ - 1] += s[k + i__ - 2] * s[k - 1];
/* L300: */
}
/* L200: */
}
for (i__ = 1; i__ <= 500; ++i__) {
rs[i__] = (rs[i__] * (xn - 1.f) + rs1[i__ - 1]) / xn;
/* L400: */
}
return 0;
} /* autocr_ */

78
tmpcsave/bpfdet.c
Wyświetl plik

@ -1,78 +0,0 @@
/* bpfdet.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include <math.h>
/* Subroutine */ int bpfdet_(real *zin, real *z)
{
/* Initialized data */
static real a[4] = { 3.0051e-5f,2.9507982f,2.90396345f,-.953135172f };
static real ck1 = 1.37158f;
static real ck2 = .9409f;
static real cg = .015f;
static real c1 = 1.2726f;
static real c2 = .81f;
static real c = .19f;
/* System generated locals */
real r;
/* Builtin functions */
// double sqrt(doublereal);
extern double Cmovavg(double a);
/* Local variables */
static real w1, w2, w3, x3, x2, x1, y3, y2, y1, y0, z3, z2, z1, zbpf, xdet;
/* THIS SUBROUTINE IMPLEMENTS THE BANDPASS FILTER AND */
/* ENVELOPE DETECTOR FUNCTIONS. THE BPF IS A SIMPLE CASCADE */
/* OF TWO 2-POLE DIGITAL RESONATORS AT A CENTER FREQ OF */
/* 1000 HZ. THE LPF OF THE ENVELOPE DETECTOR IS A */
/* THREE - POLE CHEBYSCHEV 100 HZ LPF. */
/* OFF IS TWO 2 - POLE RESONATORS : */
w3 = w2;
w2 = w1;
w1 = c1 * w2 - c2 * w3 + c * *zin;
x3 = x2;
x2 = x1;
x1 = ck1 * x2 - ck2 * x3 + cg * w1;
zbpf = x1;
/* ENVELOPE DETECTOR (SQUARE-LAW): */
/* sQUARE- */
/* Computing 2nd power */
// r = zbpf;
// xdet = sqrt(r * r);
r = *zin;
xdet = r * r; //sqrt(r * r);
xdet = (real) Cmovavg((double)xdet);// added filter here
*z = xdet;
/* LOW-PASS FILTER-
y3 = y2;
y2 = y1;
y1 = y0;
y0 = xdet * a[0];
z3 = z2;
z2 = z1;
z1 = *z;
*z = y0 + (y1 + y2) * 3.f + y3;
*z = *z + a[1] * z1 - a[2] * z2 - a[3] * 23;
return 0; */
} /* bpfdet_ */

87
tmpcsave/common.h
Wyświetl plik

@ -1,87 +0,0 @@
/* Common Block Declarations */
struct BLKTRN {
integer ltrmap[400], iblank[400];
} blktrn;
struct Alpha {
char ialph[70];
} alpha;
struct BLK1 {
real tau;
} blk1;
struct BLK6 {
real dmean;
real xdur;
real esep[6];
real edev[6];
} blk6;
struct BLK2 {
real wc, wchirp, asigma, bsigma, phisgm, rsigm, tchirp, gamma;
} blk2;
struct BLKLAM {
integer ielmst[400], ilami[16], ilamx[6];
} blklam;
struct BLKELM {
real elemtr[6][16] /* was [16][6] */;
} blkelm;
struct BLKSPD {
real rtrans[2][5] /* was [5][2] */;
integer mempr[6][6] /* was [6][6] */;
} blkspd;
struct BLKRAT {
integer memdel[6][6] /* was [6][6] */;
} blkrat;
struct BLKMEM {
integer memfcn[6][400] /* was [400][6] */;
} blkmem;
struct BLKS {
integer isx[6];
} blks;
struct BLKSV {
real ykkip[25], pkkip[25], ykksv[750], pkksv[750];
} blksv;
struct BLKEND {
integer iend;
} blkend;
struct BLKTXT {
integer itext[200];
} blktxt;
extern int path_(integer ip, integer *lambda, real *dur, integer *ilrate, integer lamsav[][25], real dursav[][25], integer ilrsav[][25]);
extern int likhd_(real *z, real *rn, integer ip, integer *lambda, real *dur, integer *ilrate, real p[][25], real lkhd[][25]);
extern int savep_(real p[][25], integer *pathsv, integer *isave, integer *imax, integer lamsav[][25], real dursav[][25], integer ilrsav[][25], integer *lambda, real *dur, integer *ilrate, integer *sort);
extern int probp_(real p[][25], real pin[][25], integer *isave, real lkhd[][25]);
extern int sprob_(real p[][25], integer *isave, integer ilrsav[][25], real *pelm, integer *khat, real *spdhat, real *px);
extern int trelis_(integer *, integer *, integer *, integer *);
extern int trprob_(integer , integer *, real *, integer *, real [][25]);

50
tmpcsave/displa.c
Wyświetl plik

@ -1,50 +0,0 @@
/* displa.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "gnuplot_i.h"
#include <stdio.h>
/* Subroutine */ int displa_(real *s1, real *s2, real *s3, real *s4)
{
/* Local variables */
static real xm, wait;
int i;
int len = 512;
static gnuplot_ctrl *h=0;
if (h == 0) {
h = gnuplot_init();
}
gnuplot_setstyle(h,"linespoints");
gnuplot_set_xlabel(h,"time (ms)");
gnuplot_plot_x(h,(double *)s1,len,"x");
// gnuplot_plot_x(h,(double *)s2,len,"z");
// gnuplot_plot_x(h,(double *)s3,len,"zsig");
// gnuplot_plot_x(h,(double *)s4,len,"zrcv");
printf("\nPress 0 and return");
scanf("%d",&i);
/* Parameter adjustments */
--s4;
--s3;
--s2;
--s1;
return 0;
} /* displa_ */

31
tmpcsave/filter.c
Wyświetl plik

@ -1,31 +0,0 @@
#include "f2c.h"
#include "malloc.h"
#include <stdio.h>
#define LEN 5
static double in[LEN];
static int empty = 1;
static int pint = 0;
static double out;
double Cmovavg(double a)
{
int i;
if (empty==1) {
empty = 0;
for (i = 0; i < LEN; i++) {
in[i] = a;
}
out = a * (double)LEN;
pint = 0;
return a;
}
out = out - in[pint] + a;
in[pint] = a;
if (++pint >= LEN) pint = 0;
return out / (double)LEN;
}

716
tmpcsave/gnuplot_i.c
Wyświetl plik

@ -1,716 +0,0 @@
/*-------------------------------------------------------------------------*/
/**
@file gnuplot_i.c
@author N. Devillard
@date Sep 1998
@version $Revision: 2.10 $
@brief C interface to gnuplot.
gnuplot is a freely available, command-driven graphical display tool for
Unix. It compiles and works quite well on a number of Unix flavours as
well as other operating systems. The following module enables sending
display requests to gnuplot through simple C calls.
*/
/*--------------------------------------------------------------------------*/
/*
$Id: gnuplot_i.c,v 2.10 2003/01/27 08:58:04 ndevilla Exp $
$Author: ndevilla $
$Date: 2003/01/27 08:58:04 $
$Revision: 2.10 $
*/
/*---------------------------------------------------------------------------
Includes
---------------------------------------------------------------------------*/
#include "gnuplot_i.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdarg.h>
#include <assert.h>
#ifdef WIN32
#include <io.h>
#endif // #ifdef WIN32
/*---------------------------------------------------------------------------
Defines
---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------
Prototype Functions
---------------------------------------------------------------------------*/
/**
* Creates a temporary file name for writing
*
* @author Peter (12/9/2011)
*
* @param handle
*
* @return char const * Pointer to file name string.
*/
char const * gnuplot_tmpfile(gnuplot_ctrl * handle);
/**
* Plot a temporary file.
*
* @author Peter (12/9/2011)
*
* @param handle
* @param tmp_filename
* @param title
*/
void gnuplot_plot_atmpfile(gnuplot_ctrl * handle, char const* tmp_filename, char const* title);
/*---------------------------------------------------------------------------
Function codes
---------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------*/
/**
@brief Opens up a gnuplot session, ready to receive commands.
@return Newly allocated gnuplot control structure.
This opens up a new gnuplot session, ready for input. The struct
controlling a gnuplot session should remain opaque and only be
accessed through the provided functions.
The session must be closed using gnuplot_close().
*/
/*--------------------------------------------------------------------------*/
gnuplot_ctrl * gnuplot_init(void)
{
gnuplot_ctrl * handle ;
int i;
#ifndef WIN32
if (getenv("DISPLAY") == NULL) {
fprintf(stderr, "cannot find DISPLAY variable: is it set?\n") ;
}
#endif // #ifndef WIN32
/*
* Structure initialization:
*/
handle = (gnuplot_ctrl*)malloc(sizeof(gnuplot_ctrl)) ;
handle->nplots = 0 ;
gnuplot_setstyle(handle, "points") ;
handle->ntmp = 0 ;
handle->gnucmd = popen("gnuplot", "w") ;
if (handle->gnucmd == NULL) {
fprintf(stderr, "error starting gnuplot, is gnuplot or gnuplot.exe in your path?\n") ;
free(handle) ;
return NULL ;
}
for (i=0;i<GP_MAX_TMP_FILES; i++)
{
handle->tmp_filename_tbl[i] = NULL;
}
return handle;
}
/*-------------------------------------------------------------------------*/
/**
@brief Closes a gnuplot session previously opened by gnuplot_init()
@param handle Gnuplot session control handle.
@return void
Kills the child PID and deletes all opened temporary files.
It is mandatory to call this function to close the handle, otherwise
temporary files are not cleaned and child process might survive.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_close(gnuplot_ctrl * handle)
{
int i ;
if (pclose(handle->gnucmd) == -1) {
fprintf(stderr, "problem closing communication to gnuplot\n") ;
return ;
}
if (handle->ntmp) {
for (i=0 ; i<handle->ntmp ; i++) {
remove(handle->tmp_filename_tbl[i]) ;
free(handle->tmp_filename_tbl[i]);
handle->tmp_filename_tbl[i] = NULL;
}
}
free(handle) ;
return ;
}
/*-------------------------------------------------------------------------*/
/**
@brief Sends a command to an active gnuplot session.
@param handle Gnuplot session control handle
@param cmd Command to send, same as a printf statement.
This sends a string to an active gnuplot session, to be executed.
There is strictly no way to know if the command has been
successfully executed or not.
The command syntax is the same as printf.
Examples:
@code
gnuplot_cmd(g, "plot %d*x", 23.0);
gnuplot_cmd(g, "plot %.18e * cos(%.18e * x)", 32.0, -3.0);
@endcode
Since the communication to the gnuplot process is run through
a standard Unix pipe, it is only unidirectional. This means that
it is not possible for this interface to query an error status
back from gnuplot.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_cmd(gnuplot_ctrl * handle, char const * cmd, ...)
{
va_list ap ;
va_start(ap, cmd);
vfprintf(handle->gnucmd, cmd, ap);
va_end(ap);
fputs("\n", handle->gnucmd) ;
fflush(handle->gnucmd) ;
return ;
}
/*-------------------------------------------------------------------------*/
/**
@brief Change the plotting style of a gnuplot session.
@param h Gnuplot session control handle
@param plot_style Plotting-style to use (character string)
@return void
The provided plotting style is a character string. It must be one of
the following:
- lines
- points
- linespoints
- impulses
- dots
- steps
- errorbars
- boxes
- boxeserrorbars
*/
/*--------------------------------------------------------------------------*/
void gnuplot_setstyle(gnuplot_ctrl * h, char * plot_style)
{
if (strcmp(plot_style, "lines") &&
strcmp(plot_style, "points") &&
strcmp(plot_style, "linespoints") &&
strcmp(plot_style, "impulses") &&
strcmp(plot_style, "dots") &&
strcmp(plot_style, "steps") &&
strcmp(plot_style, "errorbars") &&
strcmp(plot_style, "boxes") &&
strcmp(plot_style, "boxerrorbars")) {
fprintf(stderr, "warning: unknown requested style: using points\n") ;
strcpy(h->pstyle, "points") ;
} else {
strcpy(h->pstyle, plot_style) ;
}
return ;
}
/*-------------------------------------------------------------------------*/
/**
@brief Sets the x label of a gnuplot session.
@param h Gnuplot session control handle.
@param label Character string to use for X label.
@return void
Sets the x label for a gnuplot session.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_set_xlabel(gnuplot_ctrl * h, char * label)
{
gnuplot_cmd(h, "set xlabel \"%s\"", label) ;
}
/*-------------------------------------------------------------------------*/
/**
@brief Sets the y label of a gnuplot session.
@param h Gnuplot session control handle.
@param label Character string to use for Y label.
@return void
Sets the y label for a gnuplot session.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_set_ylabel(gnuplot_ctrl * h, char * label)
{
gnuplot_cmd(h, "set ylabel \"%s\"", label) ;
}
/*-------------------------------------------------------------------------*/
/**
@brief Resets a gnuplot session (next plot will erase previous ones).
@param h Gnuplot session control handle.
@return void
Resets a gnuplot session, i.e. the next plot will erase all previous
ones.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_resetplot(gnuplot_ctrl * h)
{
int i ;
if (h->ntmp) {
for (i=0 ; i<h->ntmp ; i++) {
remove(h->tmp_filename_tbl[i]) ;
free(h->tmp_filename_tbl[i]);
h->tmp_filename_tbl[i] = NULL;
}
}
h->ntmp = 0 ;
h->nplots = 0 ;
return ;
}
/*-------------------------------------------------------------------------*/
/**
@brief Plots a 2d graph from a list of doubles.
@param handle Gnuplot session control handle.
@param d Array of doubles.
@param n Number of values in the passed array.
@param title Title of the plot.
@return void
Plots out a 2d graph from a list of doubles. The x-coordinate is the
index of the double in the list, the y coordinate is the double in
the list.
Example:
@code
gnuplot_ctrl *h ;
double d[50] ;
int i ;
h = gnuplot_init() ;
for (i=0 ; i<50 ; i++) {
d[i] = (double)(i*i) ;
}
gnuplot_plot_x(h, d, 50, "parabola") ;
sleep(2) ;
gnuplot_close(h) ;
@endcode
*/
/*--------------------------------------------------------------------------*/
void gnuplot_plot_x(
gnuplot_ctrl * handle,
double * d,
int n,
char * title
)
{
int i ;
FILE* tmpfd ;
char const * tmpfname;
if (handle==NULL || d==NULL || (n<1)) return ;
/* Open temporary file for output */
tmpfname = gnuplot_tmpfile(handle);
tmpfd = fopen(tmpfname, "w");
if (tmpfd == NULL) {
fprintf(stderr,"cannot create temporary file: exiting plot") ;
return ;
}
/* Write data to this file */
for (i=0 ; i<n ; i++) {
fprintf(tmpfd, "%.18e\n", d[i]);
}
fclose(tmpfd) ;
gnuplot_plot_atmpfile(handle,tmpfname,title);
return ;
}
/*-------------------------------------------------------------------------*/
/**
@brief Plot a 2d graph from a list of points.
@param handle Gnuplot session control handle.
@param x Pointer to a list of x coordinates.
@param y Pointer to a list of y coordinates.
@param n Number of doubles in x (assumed the same as in y).
@param title Title of the plot.
@return void
Plots out a 2d graph from a list of points. Provide points through a list
of x and a list of y coordinates. Both provided arrays are assumed to
contain the same number of values.
@code
gnuplot_ctrl *h ;
double x[50] ;
double y[50] ;
int i ;
h = gnuplot_init() ;
for (i=0 ; i<50 ; i++) {
x[i] = (double)(i)/10.0 ;
y[i] = x[i] * x[i] ;
}
gnuplot_plot_xy(h, x, y, 50, "parabola") ;
sleep(2) ;
gnuplot_close(h) ;
@endcode
*/
/*--------------------------------------------------------------------------*/
void gnuplot_plot_xy(
gnuplot_ctrl * handle,
double * x,
double * y,
int n,
char * title
)
{
int i ;
FILE* tmpfd ;
char const * tmpfname;
if (handle==NULL || x==NULL || y==NULL || (n<1)) return ;
/* Open temporary file for output */
tmpfname = gnuplot_tmpfile(handle);
tmpfd = fopen(tmpfname, "w");
if (tmpfd == NULL) {
fprintf(stderr,"cannot create temporary file: exiting plot") ;
return ;
}
/* Write data to this file */
for (i=0 ; i<n; i++) {
fprintf(tmpfd, "%.18e %.18e\n", x[i], y[i]) ;
}
fclose(tmpfd) ;
gnuplot_plot_atmpfile(handle,tmpfname,title);
return ;
}
/*-------------------------------------------------------------------------*/
/**
@brief Open a new session, plot a signal, close the session.
@param title Plot title
@param style Plot style
@param label_x Label for X
@param label_y Label for Y
@param x Array of X coordinates
@param y Array of Y coordinates (can be NULL)
@param n Number of values in x and y.
@return
This function opens a new gnuplot session, plots the provided
signal as an X or XY signal depending on a provided y, waits for
a carriage return on stdin and closes the session.
It is Ok to provide an empty title, empty style, or empty labels for
X and Y. Defaults are provided in this case.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_plot_once(
char * title,
char * style,
char * label_x,
char * label_y,
double * x,
double * y,
int n
)
{
gnuplot_ctrl * handle ;
if (x==NULL || n<1) return ;
if ((handle = gnuplot_init()) == NULL) return ;
if (style!=NULL) {
gnuplot_setstyle(handle, style);
} else {
gnuplot_setstyle(handle, "lines");
}
if (label_x!=NULL) {
gnuplot_set_xlabel(handle, label_x);
} else {
gnuplot_set_xlabel(handle, "X");
}
if (label_y!=NULL) {
gnuplot_set_ylabel(handle, label_y);
} else {
gnuplot_set_ylabel(handle, "Y");
}
if (y==NULL) {
gnuplot_plot_x(handle, x, n, title);
} else {
gnuplot_plot_xy(handle, x, y, n, title);
}
printf("press ENTER to continue\n");
while (getchar()!='\n') {}
gnuplot_close(handle);
return ;
}
void gnuplot_plot_slope(
gnuplot_ctrl * handle,
double a,
double b,
char * title
)
{
char const * cmd = (handle->nplots > 0) ? "replot" : "plot";
title = (title == NULL) ? "(none)" : title;
gnuplot_cmd(handle, "%s %.18e * x + %.18e title \"%s\" with %s",
cmd, a, b, title, handle->pstyle) ;
handle->nplots++ ;
return ;
}
void gnuplot_plot_equation(
gnuplot_ctrl * h,
char * equation,
char * title
)
{
char const * cmd = (h->nplots > 0) ? "replot" : "plot";
title = (title == NULL) ? "(none)" : title;
gnuplot_cmd(h, "%s %s title \"%s\" with %s",
cmd, equation, title, h->pstyle) ;
h->nplots++ ;
return ;
}
int gnuplot_write_x_csv(
char const * fileName,
double const * d,
int n,
char const * title)
{
int i;
FILE* fileHandle;
if (fileName==NULL || d==NULL || (n<1))
{
return -1;
}
fileHandle = fopen(fileName, "w");
if (fileHandle == NULL)
{
return -1;
}
// Write Comment.
if (title != NULL)
{
fprintf(fileHandle, "# %s\n", title) ;
}
/* Write data to this file */
for (i=0 ; i<n; i++)
{
fprintf(fileHandle, "%d, %.18e\n", i, d[i]) ;
}
fclose(fileHandle) ;
return 0;
}
int gnuplot_write_xy_csv(
char const * fileName,
double const * x,
double const * y,
int n,
char const * title)
{
int i ;
FILE* fileHandle;
if (fileName==NULL || x==NULL || y==NULL || (n<1))
{
return -1;
}
fileHandle = fopen(fileName, "w");
if (fileHandle == NULL)
{
return -1;
}
// Write Comment.
if (title != NULL)
{
fprintf(fileHandle, "# %s\n", title) ;
}
/* Write data to this file */
for (i=0 ; i<n; i++)
{
fprintf(fileHandle, "%.18e, %.18e\n", x[i], y[i]) ;
}
fclose(fileHandle) ;
return 0;
}
int gnuplot_write_multi_csv(
char const * fileName,
double const ** xListPtr,
int n,
int numColumns,
char const * title)
{
int i;
int j;
FILE* fileHandle;
if (fileName==NULL || xListPtr==NULL || (n<1) || numColumns <1)
{
return -1;
}
for (j=0;j<numColumns;j++)
{
if (xListPtr[j] == NULL)
{
return -1;
}
}
fileHandle = fopen(fileName, "w");
if (fileHandle == NULL)
{
return -1;
}
// Write Comment.
if (title != NULL)
{
fprintf(fileHandle, "# %s\n", title) ;
}
/* Write data to this file */
for (i=0 ; i<n; i++)
{
fprintf(fileHandle, "%d, %.18e", i, xListPtr[0][i]) ;
for (j=1;j<numColumns;j++)
{
fprintf(fileHandle, ", %.18e", xListPtr[j][i]) ;
}
fprintf(fileHandle, "\n");
}
fclose(fileHandle) ;
return 0;
}
char const * gnuplot_tmpfile(gnuplot_ctrl * handle)
{
static char const * tmp_filename_template = "gnuplot_tmpdatafile_XXXXXX";
char * tmp_filename = NULL;
int tmp_filelen = strlen(tmp_filename_template);
#ifndef WIN32
int unx_fd;
#endif // #ifndef WIN32
assert(handle->tmp_filename_tbl[handle->ntmp] == NULL);
/* Open one more temporary file? */
if (handle->ntmp == GP_MAX_TMP_FILES - 1) {
fprintf(stderr,
"maximum # of temporary files reached (%d): cannot open more",
GP_MAX_TMP_FILES) ;
return NULL;
}
tmp_filename = (char*) malloc(tmp_filelen+1);
if (tmp_filename == NULL)
{
return NULL;
}
strcpy(tmp_filename, tmp_filename_template);
#ifdef WIN32
if (_mktemp(tmp_filename) == NULL)
{
return NULL;
}
#else // #ifdef WIN32
unx_fd = mkstemp(tmp_filename);
if (unx_fd == -1)
{
return NULL;
}
close(unx_fd);
#endif // #ifdef WIN32
handle->tmp_filename_tbl[handle->ntmp] = tmp_filename;
handle->ntmp ++;
return tmp_filename;
}
void gnuplot_plot_atmpfile(gnuplot_ctrl * handle, char const* tmp_filename, char const* title)
{
char const * cmd = (handle->nplots > 0) ? "replot" : "plot";
title = (title == NULL) ? "(none)" : title;
gnuplot_cmd(handle, "%s \"%s\" title \"%s\" with %s", cmd, tmp_filename,
title, handle->pstyle) ;
handle->nplots++ ;
return ;
}
/* vim: set ts=4 et sw=4 tw=75 */

418
tmpcsave/gnuplot_i.h
Wyświetl plik

@ -1,418 +0,0 @@
/*-------------------------------------------------------------------------*/
/**
@file gnuplot_i.h
@author N. Devillard
@date Sep 1998
@version $Revision: 1.11 $
@brief C interface to gnuplot.
gnuplot is a freely available, command-driven graphical display tool for
Unix. It compiles and works quite well on a number of Unix flavours as
well as other operating systems. The following module enables sending
display requests to gnuplot through simple C calls.
*/
/*--------------------------------------------------------------------------*/
/*
$Id: gnuplot_i.h,v 1.11 2003/01/27 08:58:04 ndevilla Exp $
$Author: ndevilla $
$Date: 2003/01/27 08:58:04 $
$Revision: 1.11 $
*/
#ifndef _GNUPLOT_PIPES_H_
#define _GNUPLOT_PIPES_H_
/*---------------------------------------------------------------------------
Includes
---------------------------------------------------------------------------*/
#include <stdio.h>
/** Maximal number of simultaneous temporary files */
#define GP_MAX_TMP_FILES 64
/*---------------------------------------------------------------------------
New Types
---------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------*/
/**
@typedef gnuplot_ctrl
@brief gnuplot session handle (opaque type).
This structure holds all necessary information to talk to a gnuplot
session. It is built and returned by gnuplot_init() and later used
by all functions in this module to communicate with the session, then
meant to be closed by gnuplot_close().
This structure is meant to remain opaque, you normally do not need
to know what is contained in there.
*/
/*-------------------------------------------------------------------------*/
typedef struct _GNUPLOT_CTRL_ {
/** Pipe to gnuplot process */
FILE * gnucmd ;
/** Number of currently active plots */
int nplots ;
/** Current plotting style */
char pstyle[32] ;
/** Pointer to table of names of temporary files */
char* tmp_filename_tbl[GP_MAX_TMP_FILES] ;
/** Number of temporary files */
int ntmp ;
} gnuplot_ctrl ;
/*---------------------------------------------------------------------------
Function ANSI C prototypes
---------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------*/
/**
@brief Opens up a gnuplot session, ready to receive commands.
@return Newly allocated gnuplot control structure.
This opens up a new gnuplot session, ready for input. The struct
controlling a gnuplot session should remain opaque and only be
accessed through the provided functions.
The session must be closed using gnuplot_close().
*/
/*--------------------------------------------------------------------------*/
gnuplot_ctrl * gnuplot_init(void);
/*-------------------------------------------------------------------------*/
/**
@brief Closes a gnuplot session previously opened by gnuplot_init()
@param handle Gnuplot session control handle.
@return void
Kills the child PID and deletes all opened temporary files.
It is mandatory to call this function to close the handle, otherwise
temporary files are not cleaned and child process might survive.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_close(gnuplot_ctrl * handle);
/*-------------------------------------------------------------------------*/
/**
@brief Sends a command to an active gnuplot session.
@param handle Gnuplot session control handle
@param cmd Command to send, same as a printf statement.
This sends a string to an active gnuplot session, to be executed.
There is strictly no way to know if the command has been
successfully executed or not.
The command syntax is the same as printf.
Examples:
@code
gnuplot_cmd(g, "plot %d*x", 23.0);
gnuplot_cmd(g, "plot %g * cos(%g * x)", 32.0, -3.0);
@endcode
Since the communication to the gnuplot process is run through
a standard Unix pipe, it is only unidirectional. This means that
it is not possible for this interface to query an error status
back from gnuplot.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_cmd(gnuplot_ctrl * handle, char const * cmd, ...);
/*-------------------------------------------------------------------------*/
/**
@brief Change the plotting style of a gnuplot session.
@param h Gnuplot session control handle
@param plot_style Plotting-style to use (character string)
@return void
The provided plotting style is a character string. It must be one of
the following:
- lines
- points
- linespoints
- impulses
- dots
- steps
- errorbars
- boxes
- boxeserrorbars
*/
/*--------------------------------------------------------------------------*/
void gnuplot_setstyle(gnuplot_ctrl * h, char * plot_style);
/*-------------------------------------------------------------------------*/
/**
@brief Sets the x label of a gnuplot session.
@param h Gnuplot session control handle.
@param label Character string to use for X label.
@return void
Sets the x label for a gnuplot session.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_set_xlabel(gnuplot_ctrl * h, char * label);
/*-------------------------------------------------------------------------*/
/**
@brief Sets the y label of a gnuplot session.
@param h Gnuplot session control handle.
@param label Character string to use for Y label.
@return void
Sets the y label for a gnuplot session.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_set_ylabel(gnuplot_ctrl * h, char * label);
/*-------------------------------------------------------------------------*/
/**
@brief Resets a gnuplot session (next plot will erase previous ones).
@param h Gnuplot session control handle.
@return void
Resets a gnuplot session, i.e. the next plot will erase all previous
ones.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_resetplot(gnuplot_ctrl * h);
/*-------------------------------------------------------------------------*/
/**
@brief Plots a 2d graph from a list of doubles.
@param handle Gnuplot session control handle.
@param d Array of doubles.
@param n Number of values in the passed array.
@param title Title of the plot.
@return void
Plots out a 2d graph from a list of doubles. The x-coordinate is the
index of the double in the list, the y coordinate is the double in
the list.
Example:
@code
gnuplot_ctrl *h ;
double d[50] ;
int i ;
h = gnuplot_init() ;
for (i=0 ; i<50 ; i++) {
d[i] = (double)(i*i) ;
}
gnuplot_plot_x(h, d, 50, "parabola") ;
sleep(2) ;
gnuplot_close(h) ;
@endcode
*/
/*--------------------------------------------------------------------------*/
void gnuplot_plot_x(gnuplot_ctrl * handle, double * d, int n, char * title);
/*-------------------------------------------------------------------------*/
/**
@brief Plot a 2d graph from a list of points.
@param handle Gnuplot session control handle.
@param x Pointer to a list of x coordinates.
@param y Pointer to a list of y coordinates.
@param n Number of doubles in x (assumed the same as in y).
@param title Title of the plot.
@return void
Plots out a 2d graph from a list of points. Provide points through a list
of x and a list of y coordinates. Both provided arrays are assumed to
contain the same number of values.
@code
gnuplot_ctrl *h ;
double x[50] ;
double y[50] ;
int i ;
h = gnuplot_init() ;
for (i=0 ; i<50 ; i++) {
x[i] = (double)(i)/10.0 ;
y[i] = x[i] * x[i] ;
}
gnuplot_plot_xy(h, x, y, 50, "parabola") ;
sleep(2) ;
gnuplot_close(h) ;
@endcode
*/
/*--------------------------------------------------------------------------*/
void gnuplot_plot_xy(
gnuplot_ctrl * handle,
double * x,
double * y,
int n,
char * title
) ;
/*-------------------------------------------------------------------------*/
/**
@brief Open a new session, plot a signal, close the session.
@param title Plot title
@param style Plot style
@param label_x Label for X
@param label_y Label for Y
@param x Array of X coordinates
@param y Array of Y coordinates (can be NULL)
@param n Number of values in x and y.
@return
This function opens a new gnuplot session, plots the provided
signal as an X or XY signal depending on a provided y, waits for
a carriage return on stdin and closes the session.
It is Ok to provide an empty title, empty style, or empty labels for
X and Y. Defaults are provided in this case.
*/
/*--------------------------------------------------------------------------*/
void gnuplot_plot_once(
char * title,
char * style,
char * label_x,
char * label_y,
double * x,
double * y,
int n
);
/*-------------------------------------------------------------------------*/
/**
@brief Plot a slope on a gnuplot session.
@param handle Gnuplot session control handle.
@param a Slope.
@param b Intercept.
@param title Title of the plot.
@return void
Plot a slope on a gnuplot session. The provided slope has an
equation of the form y=ax+b
Example:
@code
gnuplot_ctrl * h ;
double a, b ;
h = gnuplot_init() ;
gnuplot_plot_slope(h, 1.0, 0.0, "unity slope") ;
sleep(2) ;
gnuplot_close(h) ;
@endcode
*/
/*--------------------------------------------------------------------------*/
void gnuplot_plot_slope(
gnuplot_ctrl * handle,
double a,
double b,
char * title
) ;
/*-------------------------------------------------------------------------*/
/**
@brief Plot a curve of given equation y=f(x).
@param h Gnuplot session control handle.
@param equation Equation to plot.
@param title Title of the plot.
@return void
Plots out a curve of given equation. The general form of the
equation is y=f(x), you only provide the f(x) side of the equation.
Example:
@code
gnuplot_ctrl *h ;
char eq[80] ;
h = gnuplot_init() ;
strcpy(eq, "sin(x) * cos(2*x)") ;
gnuplot_plot_equation(h, eq, "sine wave", normal) ;
gnuplot_close(h) ;
@endcode
*/
/*--------------------------------------------------------------------------*/
void gnuplot_plot_equation(gnuplot_ctrl * h, char * equation, char * title) ;
/**
* Writes a CSV file for use with gnuplot commands later. Allows files to also be saved for post
* analysis with excel for example. Arguments are similar to gnuplot_plot_x()
*
* Uses title as gnuplot "comment" on the first line.
*
* @author Peter H Maresh 11/8/2011
*
* @param fileName file name to write to. This should be the same name used in the gnuplot command
* @param d
* @param n
* @param title
* @return int <0 if error writing file.
*
*/
int gnuplot_write_x_csv(
char const * fileName,
double const * d,
int n,
char const * title);
/**
* Writes a CSV file for use with gnuplot commands later. Allows files to also be saved for post
* analysis with excel for example. Arguments are similar to gnuplot_plot_xy()
*
* Uses title as gnuplot "comment" on the first line.
*
* @author Peter H Maresh 11/8/2011
*
* @param fileName file name to write to. This should be the same name used in the gnuplot command
* @param x
* @param y
* @param n
* @param title
* @return int <0 if file wasn't written.
*/
int gnuplot_write_xy_csv(
char const * fileName,
double const * x,
double const * y,
int n,
char const * title);
/**
* Writes a multi column CSV file for use with gnuplot commands later. Allows files to also be
* saved for post analysis with excel for example. Note that when used with gnuplot, since there
* may be more than one column the whole "1:3" or whatever should be used.
*
* Uses title as gnuplot "comment" on the first line.
*
* @author Peter H Maresh 11/8/2011
*
* @param fileName file name to write to. This should be the same name used in the gnuplot
* command
* @param xListPtr A list of pointers to column buffers
* @param n
* @param numColumns Length of xListPtr list. Note that the CSV file will add an additional
* "index" column first.
* @param title Title to write for the first line of the .csv file, will be preceeded by
* "#"
* @return int <0 if file wasn't written.
*/
int gnuplot_write_multi_csv(
char const * fileName,
double const ** xListPtr,
int n,
int numColumns,
char const * title);
#endif

361
tmpcsave/initl.c
Wyświetl plik

@ -1,361 +0,0 @@
/* initl.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
/* Common Block Declarations */
/* Initialized data */
struct BLKLAM blklam = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
/*ilami*/ 3, 4, 5, 6, 3, 4, 5, 6, 1, 2, 1, 2, 1, 2, 1, 2,
/*ilamx*/ 1, 1, 0, 0, 0, 0 };
// Page 103
struct BLKRAT blkrat = {
//integer memdel[6][6]
0, 0, 2, 2, 5, 10,
0, 0, 2, 2, 5, 10,
2, 2, 0, 0, 0, 0,
2, 2, 0, 0, 0, 0,
2, 2, 0, 0, 0, 0,
2, 2, 0, 0, 0, 0 };
// TABLE XII Second Order Markov Symbol Transition Matrix - Page 105 Table XII
// elemtr[6][16]
struct BLKELM blkelm = {
/* .^ .~ .w .p -^ -~ -w -p ^. ^- ~. ~- w. w- p. p-
/* . */ .55, .5, .5, .5, .55, .5, .5, .5, 0., 0., 0., 0., 0., 0., 0., 0.,
/* - */ .45,.45, .5, .5, .5, .45, .5, .5, .5, 0., 0., 0., 0., 0., 0., 0.,
/* ^ */ 0., 0., 0., 0., 0., 0., 0., 0., .581, .54, .923, .923, .923, .923, .95, .95,
/* ~ */ 0., 0., 0., 0., 0., 0., 0., 0., .335,.376, .062, .062, .062, .062, .04, .04,
/* w */ 0., 0., 0., 0., 0., 0., 0., 0., .069,.069, .012, .012, .012, .012,.009,.009,
/* p */ 0., 0., 0., 0., 0., 0., 0., 0., .015,.015, .003, .003, .003, .003,.001,.001 };
struct BLKSPD blkspd = {
// rtrans[2][5] - symbol conditional speed transition probabilities - Page 104 - Table X
// used in spdtr.c
.1, .2, .4, .2, .1, // dot, dash, e-sp, w-sp
.15,.2, .3, .2, .15, // c-sp, pause
// mempr[6][6] - [
0, 0, 1, 2, 1, 2,
0, 0, 1, 2, 1, 2,
1, 1, 0, 0, 0, 0,
1, 1, 0, 0, 0, 0,
1, 1, 0, 0, 0, 0,
1, 1, 0, 0, 0, 0 };
struct BLKMEM blkmem = { //memfcn[6][400]
/*k=0*/ 9, 11, 13, 15, 9, 11, 13, 15, 9, 0, 11, 0, 13, 0, 15, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
/*k=1*/ 10, 12, 14, 16, 10, 12, 14, 16, 0, 10, 0, 12, 0, 14, 0, 16,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
/*k=2*/ 1, 0, 0, 0, 5, 0, 0, 0, 1, 5, 1, 5, 1, 5, 1, 5,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
/*k=3*/ 0, 2, 0, 0, 0, 6, 0, 0, 2, 6, 2, 6, 2, 6, 2, 6,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
/*k=4*/ 0, 0, 3, 0, 0, 0, 7, 0, 3, 7, 3, 7, 3, 7, 3, 7,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
/*k=5*/ 0, 0, 0, 4, 0, 0, 0, 8, 4, 8, 4, 8, 4, 8, 4, 8,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
struct BLKS blks = { 1, 1, 0, 0, 0, 0 };
struct BLKSV blksv = { .5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f,
.5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f, .5f,
.1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f,
.1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f, .1f };
struct BLKTRN blktrn = { 3, 4, 5, 6, 3, 4, 5, 6, 1, 2, 1, 2, 1, 2, 1, 2, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
struct Alpha alpha = {"ABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890;:%&##############################"};
/* Table of constant values */
static integer c__1 = 1;
/* Subroutine */ int initl_(void)
{
/* Format strings */
static char fmt_30[] = "(8i3)";
static char fmt_40[] = "(10x,6(i3,2x))";
static char fmt_70[] = "(i2)";
/* System generated locals */
olist o__1;
alist al__1;
/* Builtin functions */
integer f_open(olist *), s_rsfe(cilist *), do_fio(integer *, char *,
ftnlen), e_rsfe(void), f_end(alist *), s_wsfe(cilist *), e_wsfe(
void);
/* Local variables */
static integer i, k, iarray[8];
/* Fortran I/O blocks */
static cilist io___2 = { 0, 20, 0, fmt_30, 0 };
static cilist io___5 = { 0, 20, 0, fmt_40, 0 };
static cilist io___6 = { 0, 20, 0, fmt_70, 0 };
/* /'A','B','C','D','E','F','G','H','I','J','K','L','M'/ */
/* 2 'N','O','P','Q','R','S','T','U','V','W','X','Y','Z', */
/* 2 '1','2','3','4','5','6','7'/ */
/* 2 '8','9','0',';',':','%','&',0,0,'^',',','AS','SN', */
/* 2 0,0,0,0,'NR','NO','GA','OK','AR','SK',0,0,0,0, */
/* 2 'IMI',0,0,0,0,'BT',0,0,0,'EEE'/ */
o__1.oerr = 0;
o__1.ounit = 20;
o__1.ofnmlen = 6;
o__1.ofnm = "MORSEM";
o__1.orl = 0;
o__1.osta = 0;
o__1.oacc = 0;
o__1.ofm = 0;
o__1.oblnk = 0;
f_open(&o__1);
for (i = 0; i < 300; ++i) {
s_rsfe(&io___2);
for (k = 0; k < 8; ++k) {
do_fio(&c__1, (char *)&iarray[k], (ftnlen)sizeof(integer));
}
e_rsfe();
for (k = 0; k < 6; ++k) {
blkmem.memfcn[k][i] = iarray[k + 2];
}
blktrn.ltrmap[i] = iarray[0];
blklam.ielmst[i] = iarray[1];
if (blklam.ielmst[i] == 7 || blklam.ielmst[i] == 3) {
blktrn.iblank[i] = 1;
}
if (blklam.ielmst[i] == 8 || blklam.ielmst[i] == 4) {
blktrn.iblank[i] = 2;
}
/* L10: */
}
al__1.aerr = 0;
al__1.aunit = 20;
f_end(&al__1);
o__1.oerr = 0;
o__1.ounit = 20;
o__1.ofnmlen = 6;
o__1.ofnm = "OUTPUT";
o__1.orl = 0;
o__1.osta = 0;
o__1.oacc = 0;
o__1.ofm = 0;
o__1.oblnk = 0;
f_open(&o__1);
for (k = 0; k < 6; ++k) {
s_wsfe(&io___5);
for (i = 0; i < 300; ++i) {
do_fio(&c__1, (char *)&blkmem.memfcn[k][i], (
ftnlen)sizeof(integer));
}
e_wsfe();
/* L50: */
}
al__1.aerr = 0;
al__1.aunit = 20;
f_end(&al__1);
o__1.oerr = 0;
o__1.ounit = 20;
o__1.ofnmlen = 4;
o__1.ofnm = "TEXT";
o__1.orl = 0;
o__1.osta = 0;
o__1.oacc = 0;
o__1.ofm = 0;
o__1.oblnk = 0;
f_open(&o__1);
for (i = 1; i <= 108; ++i) {
s_rsfe(&io___6);
do_fio(&c__1, (char *)&blktxt.itext[i - 1], (ftnlen)sizeof(
integer));
e_rsfe();
/* L60: */
}
al__1.aerr = 0;
al__1.aunit = 20;
f_end(&al__1);
return 0;
} /* initl_ */

96
tmpcsave/inputl.c
Wyświetl plik

@ -1,96 +0,0 @@
/* inputl.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
/* Common Block Declarations */
/* Initialized data */
struct BLK1 blk1 = { 1.25e-4f };
struct BLK6 blk6 = { 0., 0.,
{1.f, 3.f, 1.f, 3.f, 7.f, 14.f},
{ 0.f, 0.f, 0.f, 0.f, 0.f, 0.f}};
/* Subroutine */ int inputl_(void)
{
/* Format strings */
static char fmt_100[] = "(1x,\002INPUT KEYING PARMS: RATE,MEAN ELEM DURA"
"TIONS\002)";
static char fmt_150[] = "(1x,\002INPUT ELEM DURATION STD OEVIATIONS\002)";
static char fmt_300[] = "(1x,\002INPUT SIG PARMS- AVAR,BVAR,FCHIRP,TCHIR"
"P,PHIVAR\002)";
static char fmt_400[] = "(1x,\002INPUT SIG PARMS: GAMMA,FREQ,NOISE\002)";
/* Builtin functions */
integer s_wsfe(cilist *), e_wsfe(void);
double sqrt(doublereal);
/* Local variables */
static real fc, rc, avar, bvar, rate, fchirp, phivar, rnoise;
/* Fortran I/O blocks */
static cilist io___1 = { 0, 6, 0, fmt_100, 0 };
static cilist io___2 = { 0, 6, 0, fmt_150, 0 };
static cilist io___3 = { 0, 6, 0, fmt_300, 0 };
static cilist io___4 = { 0, 6, 0, fmt_400, 0 };
// s_wsfe(&io___1);
// e_wsfe();
/* READ 200,RATE,(ESEP(K), K=1,6) */
// s_wsfe(&io___2);
// e_wsfe();
/* READ 200,(EDEV(K),K=1,6) */
/* L200: */
// s_wsfe(&io___3);
// e_wsfe();
/* READ 210,AVAR,BVAR,FCHIRP,TCHIRP,PHIVAR */
/* L210: */
// s_wsfe(&io___4);
// e_wsfe();
/* READ 220, GAMMA,FC,RNOISE */
/* L220: */
/* ADDED BY MAURI */
avar = 1.f;
bvar = 1.f;
rate = 20.f;
fc = 4e3f;
fchirp = 2.f;
rnoise = 1.f;
phivar = 4.f;
blk2.tchirp = 2.f;
blk2.gamma = 1.f;
blk2.asigma = sqrt(avar);
blk2.bsigma = sqrt(bvar);
blk2.phisgm = sqrt(phivar);
blk2.rsigm = sqrt(rnoise);
blk6.dmean = 1200.f / rate;
rc = fc * 6.28319f;
blk2.wchirp = fchirp * 6.2831f;
/* IF (ESEP(1).NE.0.) GO TO 500 */
blk6.esep[0] = 1.f;
blk6.esep[1] = 3.f;
blk6.esep[2] = 1.f;
blk6.esep[3] = 3.f;
blk6.esep[4] = 7.f;
blk6.esep[5] = 14.f;
/* L500: */
return 0;
} /* inputl_ */

108
tmpcsave/kalfil.c
Wyświetl plik

@ -1,108 +0,0 @@
/* kalfil.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
int kalfil_(real *z, integer ip, real *rn, integer *ilx, integer *ixs, integer *kelem, integer *jnode, integer *israte, real *dur, integer *ilrate, real *pin, real *lkhdj)
{
/* Initialized data */
static real pinmin = 1e-4f;
/* System generated locals */
real r;
/* Builtin functions */
double sqrt(doublereal), exp(doublereal);
/* Local variables */
static real a, g, qa, hz, pz, zr, phi, pkk, ykk, pest;
extern int model_(real *, integer *, integer *, integer *, integer *, real *, real *, real *);
static real ppred, ypred, pzinv;
/* THIS SUBROUTINE COMPUTES THE ARRAY OF KALMAN FILTER */
/* RECURSIONS USED TO DETERMINE THE LIKELIHOODS. */
/* VARIABLES: */
/* Z - INPUT MEASUREMENT */
/* IP - INPUT PATH IDENTITY */
/* RN - INPUT NOISE POWER ESTIMATE */
/* ILX - INPUT SAVED KEYSTATE ON PATH IP */
/* IXS - INPUT KEYSTAT OF NEW NODE */
/* KELEM - INPUT ELEM STATE OF NEW NODE */
/* ISRATE - INPUT SPEED STATE OF NEW NODE */
/* DUR - INPUT CURRENT DURATION OF ELEMENT ON IP */
/* ILRATE - INPUT SPEED STATE ON PATH IP */
/* PIN - TRANS PROB FROM PATH IP TO NODE N */
/* LKHDJ - OUTPUT CALCULATED LIKELIHOOD VALUE */
/* SUBROUTINES USED */
/* MODEL - OBTAINS THE SIGNAL-STATE-DEPENDENT LINEAR */
/* MODEL FOR THE KALMAN FILTER RECURSIONS */
/* IF TRANSITION PROBABILITY IS VERY SMALL, DON'T */
/* BOTHER WITH LIKELIHOOD CALCULATION: */
if (*pin <= pinmin) {
*lkhdj = 0.f;
return 0;
}
/* OBTAIN STATE-DEPENDENT MODEL PARAMETERS: */
L100:
model_(dur, kelem, ilrate, israte, ixs, &phi, &qa, &hz);
/* GET PREVIOUS ESTIMATES FOR PATH IP */
ykk = blksv.ykkip[ip];
pkk = blksv.pkkip[ip];
/* IMPLEMENT KALMAN FILTER FOR THIS TRANSITION */
ypred = phi * ykk;
ppred = phi * pkk * phi + qa;
pz = hz * ppred + *rn;
pzinv = 1.f / pz;
g = ppred * hz * pzinv;
pest = (1.f - g * hz) * ppred;
zr = *z - hz * ypred;
blksv.ykksv[*jnode ] = ypred + g * zr;
blksv.pkksv[*jnode ] = pest;
if (blksv.ykksv[*jnode ] <= .01f) {
blksv.ykksv[*jnode ] = .01f;
}
/* Computing 2nd power */
r = zr;
a = pzinv * .5f * (r * r);
if (a <= 1e3f) {
goto L200;
}
*lkhdj = 0.f;
goto L400;
L200:
*lkhdj = 1.f / sqrt(pz) * exp(-a);
L400:
#ifdef DEBUG
printf("\nKALFIL: lkhdj:%f pin:%f phi:%f qa:%f hz:%f a:%f jnode:%d", *lkhdj, *pin, phi,qa,hz,a, *jnode);
#endif
return 0;
} /* kalfil_ */

156
tmpcsave/key.c
Wyświetl plik

@ -1,156 +0,0 @@
/* key.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
#include <math.h>
//#define DEBUG 1
/* Common Block Declarations */
/* Table of constant values */
static integer c__1 = 1;
/* Subroutine */ int key_(real *dur, real *x)
{
/* Initialized data */
static integer ik = 0.5;
static integer row = 0;
static integer col = 0;
static integer n = 0;
static integer nltr = 1;
static integer morse[10][40] = {
1,3,2,0,0,0,0,0,0,0, 2,3,1,3,1,3,1,0,0,0, 2,3,1,3,2,3,1,0,0,0, 2,3,1,3,1,0,0,0,0,0, // ABCD
1,0,0,0,0,0,0,0,0,0, 1,3,1,3,2,3,1,0,0,0, 2,3,2,3,1,0,0,0,0,0, 1,3,1,3,1,3,1,0,0,0, // EFGH
1,3,1,0,0,0,0,0,0,0, 1,3,2,3,2,3,2,0,0,0, 2,3,1,3,2,0,0,0,0,0, 1,3,2,3,1,3,1,0,0,0, // IJKL
2,3,2,0,0,0,0,0,0,0, 2,3,1,0,0,0,0,0,0,0, 2,3,2,3,2,0,0,0,0,0, 1,3,2,3,2,3,1,0,0,0, // MNOP
2,3,2,3,1,3,2,0,0,0, 1,3,2,3,1,0,0,0,0,0, 1,3,1,3,1,0,0,0,0,0, 2,0,0,0,0,0,0,0,0,0, // QRST
1,3,1,3,2,0,0,0,0,0, 1,3,1,3,1,3,2,0,0,0, 1,3,2,3,2,0,0,0,0,0, 2,3,1,3,1,3,2,0,0,0, // UVWX
2,3,1,3,2,3,2,0,0,0, 1,3,1,3,2,3,2,3,2,0, 1,3,2,3,2,3,2,3,2,0, 1,3,1,3,2,3,2,3,2,0, // YZ12
1,3,1,3,1,3,2,3,2,0, 1,3,1,3,1,3,1,3,2,0, 1,3,1,3,1,3,1,3,1,0, 2,3,1,3,1,3,1,3,1,0, // 3456
2,3,2,3,1,3,1,3,1,0, 2,3,2,3,2,3,1,3,1,0, 2,3,2,3,2,3,2,3,1,0, 2,3,2,3,2,3,2,3,2,0, // 7890
0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0 };
// 6 morse symbols: 1=dot(.) 2=dash(-) 3=element space(^) 4=char space (~) 5=word space(W) 6=pause(P)
static char isymbl[6] = { ".-^~WP"};
/* Local variables */
static integer k, l;
static real y;
static integer k1, k2, iy;
static real xm, beta;
extern double randn(real mean, real sigma);
static integer ielm, kait;
static real tend, xsigm;
static char ioutp[500];
beta = 1e3f * blk1.tau * *dur;
if (beta < blk6.xdur) {
goto L200;
}
// Get Morse element from morse table [row][col]
++col;
ielm = morse[row][col];
if (ielm != 0) {
goto L100;
}
col = 0;
y = randn(ik,2);
ielm = 4;
if (y > .9f) {
ielm = 5;
}
if (y <= .9f && y > .3f) {
ielm = 6;
}
// implement random characters here
// y = rand(36);
// y = (y - .001f) * 35 + 1.f;
// iy = y;
// ltr = iy + 1;
row +=1;
if (row >= 10)
row =0;
goto L100;
// Fetch letter from pre-loaded text - not working now
++nltr;
row = blktxt.itext[nltr - 1];
if (row == 0) {
ielm = 4;
}
if (row == 37) { // %
ielm = 5;
}
if (row == 38) { // &
ielm = 6;
}
++nltr;
row = blktxt.itext[nltr - 1];
L100:
// collect morse symbols into buffer ioutp until 300 and then print
++n;
ioutp[n - 1] = isymbl[ielm - 1];
#ifdef DEBUG
printf("\nn:%d isymbl[%d]:%c col:%d row:%d",n,(ielm-1),isymbl[ielm-1],col,row);
#endif
if (n >= 300) {
printf("\nEND OF RUN; INPUT DATA WAS: %s", ioutp);
// for (n=0; n<301; n++)
// ioutp[n] = 0;
n = 0;
}
L110:
xm = blk6.esep[ielm - 1] * blk6.dmean;
xsigm = blk6.edev[ielm - 1] * blk6.dmean;
// implement random duration later
// y = randn(ik,2);
// y = (y - .5f) * 2.f;
// blk6.xdur = xm + y + xsigm;
blk6.xdur = xm + xsigm;
if (blk6.xdur < 20.f) {
blk6.xdur = 20.f;
}
// OUTPUT: mark = 1.0 space = 0.0
*x = 1.f;
if (ielm >= 3) {
*x = 0.f;
}
L200:
return 0;
} /* key_ */
#undef isymbl

89
tmpcsave/likhd.c
Wyświetl plik

@ -1,89 +0,0 @@
/* likhd.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
int likhd_(real *z, real *rn, integer ip, integer *lambda, real *dur, integer *ilrate, real p[][25], real lkhd[][25])
{
/* Local variables */
static integer i, j, k, n;
static real pin;
static integer ilx, ixs;
static real lkhdj;
static integer kelem;
extern int kalfil_(real *, integer , real *, integer *,integer *, integer *, integer *, integer *, real *, integer *, real *, real *);
static integer israte;
/* THIS SUBROUTINE CALCULATES,FOR EACH PATH */
/* EXTENSION TO STATE N, THE LIKELIHOOD OF THAT */
/* TRANSITION GIVEN THE MEASUREMENTZ. IT USES */
/* AN ARRAY OF LINEAR (KALMAN) FILTERS TO DO SO. */
/* VARIABLES: */
/* Z- INPUT MEASUREMENT */
/* BN-INPUT NOISE POWER ESTIMATE */
/* IP-INPUT SAVED PATH IDENTITY */
/* LAMBDA-INPUT SAVED LTR STATE IDENTITY */
/* DUR-INPUT SAVED DURATION OF ELEMENT ON PATH IP */
/* ILRATE-INPUT SAVED DATA RATE (SPEED) */
/* P-INPUT TRANSITION PROSABILITIES */
/* SBHD-OUTPUT COMPUTED LIKELIHOODS FOR EACH TRANS */
/* SUBRODTINES USED: */
/* KALFIL-KALMAN FILTER FOR EACH NEW PATH */
/* OBTAIN SAVED KEYSTATE: */
/* Function Body */
if (*lambda == 0) {
goto L200;
}
kelem = blklam.ilami[blklam.ielmst[*lambda - 1] - 1];
ilx = blklam.ilamx[kelem - 1];
/* FOR EACH STATE: */
for (k = 0; k < 6; ++k) {
for (i = 0; i < 5; ++i) {
/* OBTAIN KEYSTATE, RATE STATE, STATE N, NEW NODE: */
// n = (i - 1) * 6 + k;
// j = (ip ) * 30 + n;
ixs = blks.isx[k];
israte = i;
n = i * 6 + k;
j = ip * 30 + n;
pin = p[n][ip];
/* COMPUTE AND STORE LIKELIHOOD: */
kalfil_(z, ip, rn, &ilx, &ixs, &kelem, &j, &israte, dur, ilrate,&pin, &lkhdj);
lkhd[n][ip] = lkhdj;
#ifdef DEBUG
printf("\nLIKHD: p[%2d][%2d]:%f z:%f lambda:%2d k:%d ilrate:%d israte:%d dur:%f pin:%f lkhdj[%3d][%3d]:%f rn:%f ip:%d n:%d i:%d k:%d",n,ip,p[n][ip],*z,*lambda,k, *ilrate,israte,*dur,pin,n,ip, lkhd[n][ip],*rn, ip, n, i, k);
#endif
goto L100;
if (pin <= 1e-6f) {
goto L100;
}
L100:
;
}
}
L200:
return 0;
} /* likhd_ */

82
tmpcsave/model.c
Wyświetl plik

@ -1,82 +0,0 @@
/* model.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
/* Table of constant values */
static doublereal c_b5 = 10.;
/* Subroutine */ int model_(real *dur, integer *ielm, integer *ilr, integer *isr, integer *ixs, real *phi, real *qa, real *hz)
{
/* System generated locals */
doublereal d1;
/* Builtin functions */
double exp(doublereal), pow_dd(doublereal *, doublereal *);
/* Local variables */
static real r1;
static integer iflm;
static real bauds, xsamp;
/* THIS SUBROUTINE COMPUTES THE PARAMETERS OF THE */
/* OBSERVATION STATE TRANSITION MATRIX PHI, THE */
/* MEASUREMENT MATRIX, AND THE COVARIANCES. */
/* VARIABLES: */
/* DUR- INPUT ELEMENT DURATION */
/* IELM- INPUT ELEMENT TYPE */
/* ILR- INPUT SAVED RATE */
/* ISR- INPUT RATE OF NEW STATE */
/* IXS- INPUT KEYSTATE OF NEW STATE */
/* PHIA- OUTPUT STATE TRANSITION MATRIX ENTRY FOR SIGNAL AMPLITUDE STATE */
/* QA- OUTPUT COVARIANCE FOR AMPLITUDE STATE */
/* HZ- OUTPUT MEASUREMENT MATRIX VALUE */
/* COMPUTE MEASUREMENT COEFFICIENT: */
*hz = (real) (*ixs);
/* COMPUTE PHI AND AMPLITUDE STATE VARIANCE (Q): */
r1 = 1200.f / *ilr;
bauds = *dur / r1;
if (bauds >= 14.f) {
bauds = 14.f;
}
if (*ielm >= 2) {
goto L100;
}
*qa = 1e-4f;
*phi = 1.f;
goto L300;
L100:
if (*ixs == 0) {
goto L200;
}
*phi = 1.f;
*qa = exp((bauds - 14.f) * .6f) * .15f;
*qa += bauds * .01f * exp((1.f - bauds) * .2f);
goto L300;
L200:
xsamp = r1 * 22.4f;
d1 = (doublereal) (-2 / xsamp);
*phi = pow_dd(&c_b5, &d1);
if (bauds >= 14.f) {
*phi = 1.f;
}
*qa = 0.f;
L300:
return 0;
} /* model_ */

104
tmpcsave/morse.c
Wyświetl plik

@ -1,104 +0,0 @@
/* morse.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
/* Main program */ int MAIN__(void)
{
/* Initialized data */
static real rn = .1f;
static integer np = 0;
/* Builtin functions */
/* Subroutine */ int s_stop(char *, ftnlen);
/* Local variables */
static integer n;
static real x, z;
static integer n1, n2;
static real s1[512], s2[512], s3[512], s4[512], dp, px;
static integer xhat;
static real zdet, zsig;
extern /* Subroutine */ int rcvr_(real *, real *);
static real zrcv;
extern /* Subroutine */ int noise_(real *, real *, real *), initl_(void),
stats_(real *, real *, real *, real *, real *, real *, real *,
real *, integer *), bpfdet_(real *, real *);
static integer elmhat;
extern /* Subroutine */ int simsgi_(real *, real *), proces_(real *, real
*, integer *, real *, integer *, integer *);
static integer ltrhat;
extern /* Subroutine */ int inputl_(void);
static int sample = 0;
FILE *fp;
fp = fopen("data.txt","r");
if (fp == NULL) {
perror("Error opening file");
return (-1);
}
initl_();
inputl_();
#ifdef DEBUG
printf("\nsample\tx\tzsig\tzrcv\tzdet\tz\trn\txhat\tpx");
#endif
L1:
for (n1 = 1; n1 <= 512; ++n1) {
for (n2 = 1; n2 <= 18; ++n2) {
//
// if (fscanf(fp,"%f",&x) !=1) {
// return 0;
// }
sample++;
//printf("\nMORSE: n1:%d n2:%d sample:%d x:%f", n1,n2, sample,x);
simsgi_(&x, &zsig);
rcvr_(&zsig, &zrcv);
bpfdet_(&zrcv, &zdet);
++np;
if (np < 40) { // decimate by 40: 4000Hz -> 100Hz
goto L3;
}
np = 0;
noise_(&zdet, &rn, &z);
proces_(&x, &rn, &xhat, &px, &elmhat, &ltrhat);
#ifdef DEBUG
// printf("\nMORSE:n1:%d\tx:%f\tzsig:%f\tzrcv:%f\tzdet:%f\tz:%f\trn:%fxhat:%d\tpx:%f",n1,x,zsig,zrcv,zdet,z,rn,xhat,px);
printf("\n%d\t%f\t%f\t%f\t%f\t%f\t%f\t%d\t%f",sample,x,zsig,zrcv,zdet,z,rn,xhat,px);
#endif
L3:
;
}
n = n1;
stats_(&x, &z, &zsig, &zdet, s1, s2, s3, s4, &n);
/* L2: */
}
// displa_(s1, s2, s3, s4);
goto L1;
s_stop("", (ftnlen)0);
return 0;
} /* MAIN__ */
/* Main program alias */ int morse_ () { MAIN__ (); return 0; }

114
tmpcsave/noise.c
Wyświetl plik

@ -1,114 +0,0 @@
/* noise.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
int noise_(real *zin, real *rn, real *z__)
{
/* Initialized data */
static real ylong[200] = { 1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f };
static real yshort[50] = { 1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,
1.f,1.f,1.f,1.f,1.f,1.f,1.f };
static integer kl = 1;
static integer kkl = 1;
static integer ks = 1;
static integer kks = 1;
static real ymin1 = 1.f;
static real ymin2 = 1.f;
static real ymavg = .05f;
/* System generated locals */
integer i1;
/* Local variables */
static integer i;
static real fs, ymin;
/* THIS SUBROUTINE ESTIMATES THE NOISE POWER IN THE */
/* ENVELOPE DETECTED OUTPUT FOR USE BY THE KALMAN */
/* FILTERS. AN ESTIMATE OF THE NOISE POWER IS */
/* ALSO SUBTRACTED FROM THE ENVELOPE DETECTED SIGNAL */
/* IN ORDER TO PRODUCE A ZERO-MEAN NOISE PROCESS */
/* AT THE INPUT TO THE MORSE PROCESSOR. */
++kl;
if (kl == 201) {
kl = 1;
}
fs = (real) (ks + 1);
if (ks == 51) {
ks = 1;
}
++kkl;
if (kkl >= 200) {
kkl = 200;
}
++kks;
if (kks >= 50) {
kks = 50;
}
if (kks <= 2) {
goto L10;
}
ylong[kl - 1] = *zin;
yshort[ks - 1] = *zin;
ymin1 = *zin;
ymin2 = *zin;
L10:
i1 = kkl;
for (i = 1; i <= i1; ++i) {
if (ylong[i - 1] > ymin1) {
goto L100;
}
ymin1 = ylong[i - 1];
L100:
;
}
i1 = kks;
for (i = 1; i <= i1; ++i) {
if (yshort[i - 1] > ymin2) {
goto L200;
}
ymin2 = yshort[i - 1];
L200:
;
}
ymin = ymin1;
if (ymin2 < ymin1) {
ymin = ymin2;
}
ymavg += (ymin - ymavg) * .004f;
*rn = ymavg * .3f;
if (*rn < .005f) {
*rn = .005f;
}
*z__ = (*zin - ymavg * 2.4f - .05f) * 1.1f;
return 0;
} /* noise_ */

84
tmpcsave/path.c
Wyświetl plik

@ -1,84 +0,0 @@
/* path.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
/* Common Block Declarations */
int path_(integer ip, integer *lambda, real *dur, integer *ilrate, integer lamsav[][25], real dursav[][25], integer ilrsav[][25])
{
static integer i, j, k, n, ixl, ixs, ilelm;
/* PATH COMPUTES THE LTR STATE, DURATION, AND DATA RATE OF */
/* EACH NEW PATH EXTENDED TO STATE N */
/* VARIABLES: */
/* IP- SAVED PATH IDENTITY */
/* LAMBDA LTR STATE OF SAVED PATH */
/* DUR DURATION OF ELEMENT ON SAVED PATH */
/* ILRATE DATA RATE OF ELEMENT ON SAVED PATH */
/* LAMSAV NEW LTR STATES FOR EACH PATH EXTENSION */
/* DURSAV NEW ELEM DURATIONS FOR EACH PATH EXTENSION */
/* ILRSAV NEW DATA RATES FOR EACH PATH EXTENSION */
/* J NEW PATH IDENTITY */
/* THE LETTER TRANSITION TABLE, MEMFCN, IS STORED IN COMMON. */
/* FOR EACH ELEM STATE K, AND EACH SPEED I, COMPUTE: */
for (k = 0; k < 6; ++k) { // ELEM STATE K
for (i = 0; i < 5; ++i) { // SPEED I
/* STATE IDENTITY N: */
n = (i ) * 6 + k;
/* NEW PATH IDENTITY: */
j = (ip) * 30 + n;
/* NEW LTR STATE: */
if (*lambda == 0) {
lamsav[n][ip] = 0;
goto L100;
}
lamsav[ip][n] = blkmem.memfcn[k][*lambda-1];
if (lamsav[n][ip] == 0) {
goto L100;
}
/* NEW DURATION: */
/* OBTAIN KEYSTATE OF SAVED PATH AND NEW STATE: */
ilelm = blklam.ilami[blklam.ielmst[*lambda - 1] - 1];
ixl = blklam.ilamx[ilelm - 1];
ixs = blks.isx[k];
/* CALCULATE DURATION: Phi_k number of samples since last state transition
/ Page 90 Phi_k = Phi_k-1(1- x_k - x_k-1 + 2* x_k*x_k-1) +1 */
dursav[n][ip] = *dur * (1 - ixs - ixl + 2*ixs*ixl) * 5.f;
/* NEW DATA RATE: */
ilrsav[n][ip] = *ilrate + (i - 2) * blkrat.memdel[k][ilelm];
L100:
#ifdef DEBUG
printf("\nPATH: ilrsav[%2d][%2d]:%2d ip:%2d n:%2d k:%d i:%d ilrate:%2d ilelm:%d ixl:%d ixs:%d \t memdel[%d][%d]:%d dursav[%3d][%3d]:%f",
n,ip, ilrsav[n][ip], ip, n, k, i, *ilrate, ilelm, ixl, ixs, k, ilelm, blkrat.memdel[k][ilelm],n,ip,dursav[n][ip]);
#endif
;
}
}
return 0;
} /* path_ */

64
tmpcsave/probp.c
Wyświetl plik

@ -1,64 +0,0 @@
/* probp.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
int probp_(real p[][25], real pin[][25], integer *isave, real lkhd[][25])
{
/* Local variables */
static integer i, j, n;
static real pmax, psav[30][25], psum;
/* PROBP COMPUTES THE POSTERIOR PROBABILITY OF EACH NEW PATH */
/* VARIABLES: */
/* P- INPUT: SAVED PROBS OF PRIOR PATHS */
/* OUTPUT: COMPUTED POSTERIOR PROBS OF NEW PATHS */
/* PIN- INPUT TRANSISTION PROBABILITIES */
/* LKHD- INPUT LIKLIEHOODS OF EACH TRANSTION */
/* PSUM- NORMALIZING CONSTANT (SUM OF P(J)) */
/* Function Body */
pmax = 0.f;
psum = 0.f;
/* FOR EACH SAVED PATH, EACH TRANSITION: */
for (i = 0; i < *isave; ++i) {
for (n = 0; n < 30; ++n) {
/* COMPUTE IDENTITY OF NEW PATH: */
j = i*30 + n;
/* PRODUCT OF PROBS, ADD TO PSUM */
psav[n][i] = p[n][i] * pin[n][i] * lkhd[n][i];
psum += psav[n][i];
#ifdef DEBUG
printf("\nPROBP: psum:%f psav[%3d]:%f p[%d][%d]:%f pin[%2d][%2d]:%f lkhd[%3d]:%f", psum, j,psav[j],n,i,p[n][i],n,i,pin[n][i],j,lkhd[j]);
#endif
if (psav[n][i] > pmax) {
pmax = psav[n][i];
}
}
}
/* NORMALIZE TO GET PROBABILITIES; SAVE: */
for (j = 0; j < *isave; ++j) {
for (n = 0; n < 30; n++) {
p[n][j] = psav[n][j] / psum;
}
}
return 0;
} /* probp_ */

144
tmpcsave/proces.c
Wyświetl plik

@ -1,144 +0,0 @@
/* proces.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
#define DEBUG 1
int proces_(real *z, real *rn, integer *xhat, real *px, integer *elmhat, integer *ltrhat)
{
/* Initialized data */
static integer isave = 25;
static integer lambda[25] = { 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5 };
static integer ilrate[25] = { 10,10,10,10,10,20,20,20,20,20,30,30,30,30,
30,40,40,40,40,40,50,50,50,50,50 };
static real dur[25] = { 1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,
1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,1e3f,
1e3f,1e3f };
/* Local variables */
static integer i, k, inl;
static integer init = 0;
static real spdhat;
static real pelm;
static integer khat;
static real pin[30][25]; // N: 5 spd x 6 morse states I: 25 paths - transition probability from path I to state N
static real lkhd[30][25]; // 5 speeds x 6 morse states x 25 paths = 750 likelyhoods
static integer ilrsav[30][25];
static real dursav[30][25];
static integer lamsav[30][25];
static real p[30][25];
static integer imax, sort[25];
static integer pathsv[25];
static integer ipath;
if (init == 0) {
init = 1;
for (k= 0; k< 30; k++) {
for (i = 0; i<25; i++) {
ilrsav[k][i] = 0;
lamsav[k][i] = 5;
p[k][i] = 1.0;
}
}
}
/* THIS SUBROUTINE IMPLEMENTS THE PROCESSING ALGORITHM */
/* FOR JOINT DEMODULATION, DECODING, AND TRANSLATION OF */
/* THE RECEIVED MORSE PROCESS. IT TAKES IN A NEW MEASUREMENT, Z */
/* OF THE DETECTED SIGNAL EVERY 5 MSEC AND PRODUCES AN ESTIMATE */
/* OF THE CURRENT KEYSTATE, ELEMENT STATE, AND LETTER OF THE RECEIVED SIGNAL. */
/* DEFINITIONS OF VARIABLE NAMES: */
/* Z- INPUT SAMPLE OF DETECTED SIGNAL */
/* RN- INPUT NOISE POWER ESTIMATE */
/* XHAT- OUTPUT ESTIMATE OF KEYSTATE */
/* ELMHAT- OUTPUT ESTIMATE OF ELEMENT STATE */
/* LTRHAT- OUTPUT ESTIMATE OF LETTER STATE */
/* ISAVE- NO. OF PREVIOUS PATHS SAVED */
/* IPATH- IDENTITY OF SAVED PATH */
/* LAMBDA (I) IDENTITY OF LTR STATE OF SAVED PATH I */
/* DUR (I)- DURATION OF ELEMENT ON PATH I */
/* ILRATE(I)- IDENTITY OF DATA RATE ON PATH I */
/* PIN (I,N)- COMPUTED TRANS PROB FROM PATH I TO STATE N */
/* LAMSAV(J)- IDENTITY OF LTR STATE AT NEW NODE J */
/* ILRSAV(J)- IDENTITY OF DATA RATE AT NEW NODE J */
/* LKHD(J)- LIKELIHOOD VALUE FOR NODE J */
/* P(J)- COMPUTED POSTERIOR PROB OF PATH ENDING AT NEW NODE J */
/* PSELEM(K)- COMPUTED POSTERIOR PROB OF ELEM K */
/* SPDHAT - COND MEAN ESTIMATE OF INSTANT DATA RATE */
/* PX- POSTERIOR PROB THAT KEYSTATE EQUALS 1 */
/* THE FOLLOWING SUBROUTINES ARE UTILIZED: */
/* TRPROB- COMPUTES TRANSITION PROBABILITIES */
/* PATH- COMPUTES IDENTITY OF NEW PATHS */
/* LIKHD- COMPUTES THE LIKELIHOOD OF EACH PATH EXTENSION */
/* PROBP- COMPUTES POSTERIOR PROBS OF EACH NEW PATH */
/* SPROB- COMPUTES POSTERIOR PROBS OF EACH STATE */
/* SAVE- SAVES THE HIGHEST PROB PATHS */
/* TRELIS- FORMS A TRELIS OF SAVED PATHS */
/* TRANSL- TRANSLATES THE LETTER ESTIMATE */
/* ALL TABLES OF CONSTANTS ARE STORED IN COMMON. */
/* FOR EACH SAVED PATH, COMPUTE: */
/* TRANSITION PROBABILITY TO NEW STATE (TRPROB); */
/* IDENTITY OF EACH NEW PATH EXTENDED (PATH); */
/* LIKELIHOOD OF EACH STATE EXTENSION (LIKHD); */
for (k= 0; k< 30; k++) {
for (i = 0; i<25; i++) {
p[k][i] = 1.0;
}
}
for (i = 0; i < isave; ++i) {
trprob_(i, &lambda[i], &dur[i], &ilrate[i], pin);
path_(i, &lambda[i], &dur[i], &ilrate[i], lamsav, dursav, ilrsav);
likhd_(z, rn, i, &lambda[i], &dur[i], &ilrate[i], pin, lkhd);
}
/* HAVING OBTAINED ALL NEW PATHS, COMPUTE: */
/* POSTERIOR PROBABILITY OF EACH NEW PATH(PROBP); */
/* POSTERIOR PROBABILITY OF KEYSTATE, ELEM STATE, */
/* CONDITIONAL MEAN ESTIMATE OF SPEED(SPROB); */
probp_(p, pin, &isave, lkhd);
sprob_(p, &isave, ilrsav, &pelm, &khat, &spdhat, px);
#ifdef DEBUG
printf("\nPROCES: z:%f\tpx:%f\tspdhat:%3.2f\tisave:%d\tpelm:%f\tkhat:%d",*z, *px,spdhat, isave, pelm, khat);
#endif
*xhat = 0;
if (*px > .5f) {
*xhat = 1;
}
/* SAVE THE PATHS WITH HIGHEST PROBABILITY, AND STORE THE VALUES CORRESPONDING TO THESE PATHS: */
savep_(p, pathsv, &isave, &imax, lamsav, dursav, ilrsav, lambda, dur, ilrate, sort);
/* UPDATE TRELLIS WITH NEW SAVED NODES, AND OBTAIN LETTER STATE ESTIMATE: */
//trelis_(&isave, pathsv, lambda, &imax);
return 0;
} /* proces_ */

68
tmpcsave/ptrans.c
Wyświetl plik

@ -1,68 +0,0 @@
/* ptrans.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
int ptrans_(integer *kelem, integer *irate, integer *lambda, integer *ilrate, real *ptrx, real *psum, real *pin, integer *n)
{
static real pelem, prate;
extern doublereal spdtr_(integer *, integer *, integer *, integer *);
/* THIS FUNCTION SUBROUTINE RETURNS THE PATH CONDITIONAL TRANSITION */
/* PROBABILITIES TO EACH ALLOWABLE STATE N. */
/* VARIABLES: */
/* KELEM- INPUT CURRENT ELEMENT STATE */
/* IRATE- INPUT CURRENT DATA RATE STATE */
/* LAMBDA- INPUT IDENTITY OF CURRENT LTR STATE */
/* PTRX- INPUT KEYSTATE TRANSITION PROBABILITY */
/* ELEMTR- ELEMENT TRANSITION PROBABILITY MATRIX */
/* FUNCTION SUBROTINE USED: */
/* SPDTR- RETURNS DATA RATE TRANSITION PROGS,CONDITIONED ON CURRENT SPACE TYPE. */
/* IF THE SAVED ELEMENT AND THE ELEMENT OF THE STATE N */
/* TO WHICH THE PATH IS BEING EXTENDED ARE THE SAME, */
/* THEN THE STATE TRANS PROB IS SIMPLY KEYSTATE TRANS PROB: */
if (*kelem != blklam.ilami[blklam.ielmst[*lambda - 1] - 1]) {
goto L100;
}
pin[*n] = *ptrx;
if (*irate != 2) {
pin[*n] = 0.f;
}
goto L200;
/* OTHERWISE: */
/* OBTAIN ELEM TRANS PROBS TABLE: */
L100:
pelem = blkelm.elemtr[*kelem][blklam.ielmst[*lambda - 1]];
/* NEXT COMPUTE ELEM-CONDITIONAL SPEED TRANS PROB: */
// iselm= next state ilelm=current state
prate = spdtr_(irate, ilrate, kelem, &blklam.ilami[blklam.ielmst[*lambda - 1] - 1]);
/* TRANS IS THE PRODUCT: */
pin[*n] = (1.f - *ptrx) * pelem * prate;
L200:
*psum += pin[*n];
#ifdef DEBUG
printf("\nPTRANS: psum:%f pin[%2d]:%f pelem[%d][%d]:%f prate:%f ptrx:%f lambda:%d irate:%d", *psum, *n, pin[*n], *kelem, blklam.ielmst[*lambda - 1], pelem,prate,*ptrx, *lambda, *irate);
#endif
return 0;
} /* ptrans_ */

23
tmpcsave/randn.c
Wyświetl plik

@ -1,23 +0,0 @@
#include "f2c.h"
#include <stdio.h>
#include <math.h>
#define RAND_MAX 2147483647
double randn(real mean, real sigma)
{
double u = ((double) rand() / (RAND_MAX)) * 2 - 1;
double v = ((double) rand() / (RAND_MAX)) * 2 - 1;
double r = u * u + v * v;
if (r == 0 || r > 1) {
// printf(".");
return randn(mean,sigma);
}
double c = sqrt(-2 * log(r) / r);
return (mean + sigma * u * c);
}

52
tmpcsave/rcvr.c
Wyświetl plik

@ -1,52 +0,0 @@
/* rcvr.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
/* Common Block Declarations */
/* Table of constant values */
static real c_b2 = 6.28319f;
/* Subroutine */ int rcvr_(real *zin, real *zout)
{
/* Initialized data */
static real theta = 0.f;
static real thetlo = 0.f;
/* Builtin functions */
double r_mod(real *, real *), cos(doublereal), sin(doublereal);
/* Local variables */
static real z1, zq, zglp, zilp, zqlp;
/* THIS SUBROUTINE CONVERTS THE INPUT SIGNAL AT */
/* RADIAN FREQ WC TO 1000 Hz. */
theta += blk2.wc * blk1.tau;
theta = r_mod(&theta, &c_b2);
z1 = *zin * cos(theta);
zq = *zin * sin(theta);
zilp += (z1 - zilp) * .07f;
zqlp += (zq - zqlp) * .07f;
thetlo += blk1.tau * 6283.2f;
thetlo = r_mod(&thetlo, &c_b2);
*zout = zilp * cos(thetlo) + zglp * sin(thetlo);
return 0;
} /* rcvr_ */

214
tmpcsave/savep.c
Wyświetl plik

@ -1,214 +0,0 @@
/* savep.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
int savep_(real p[][25], integer *pathsv, integer *isave, integer *imax, integer lamsav[][25], real dursav[][25], integer ilrsav[][25], integer *lambda, real *dur, integer *ilrate, integer *sort)
{
/* Initialized data */
static real popt = .9f;
/* Local variables */
static integer i, j, k, n;
static integer ip, jsav;
static integer nsav,savn,nsave[25];
static integer iconv[25], ipsav, isavm1, nplus1;
static real pmax;
static real psav[25], psum;
/* THIS SUBROUTINE PERFORMS THE ALGORITM TO SAVE */
/* THE PATHS WITH HIGHEST POSTERIOR PROBABILITY. */
/* IT WILL SAVE A MINIMUM OF 7 PATHS (ONE FOR EACH */
/* ELEMENT STATE AND ONE ADDITIONAL NODE), AND */
/* A MAXIMUM OF 25 PATHS. WITHIN THESE LIMITS, IT */
/* SAVED ONLY ENOUGH TO MAKE THE TOTAL SAVED PROBABILITY */
/* EQUAL TO POPT. */
/* ADDITIONALLY, IT RESORTS THE LAMBDA,DUR,AND ILRATE */
/* ARRAYS TO CORRESPOND TO THE SAVED NODES. */
/* VARIABLES: */
/* P- INPUT PROBABILITY ARRAY OF NEW NODES */
/* PATHSV- OUTPUT ARRAY OF THE PREVIOUS NODES TO */
/* WHICH THE SAVED NODES ARE CONNECTED. */
/* ISAVE- INPUT: NO. OF PREVIOUS NODES SAVED */
/* OUPUT: NO. OF NODES SAVED AT CURRENT STAGE */
/* IMAX- INDEX OF HIGHEST PROBABILITY NODE */
/* LAMSAV- INPUT ARRAY OF LTR STATES AT EACH NEW NODE */
/* DURSAV- INPUT ARRAY OF SAVED DURATIONS */
/* ILRSAV- INPUT ARRAY OF SAVED RATES */
/* LAMBDA- OUTPUT ARRAY OF SAVED LTR STATES, SORTED ACCORDING TO PROBABILITY */
/* DUR- OUTPUT ARRAY OF SORTED DURATIONS */
/* ILRATE- OUTPUT ARRAY OF SORTED RATES */
/* Function Body */
nsav = 0;
psum = 0.f;
jsav = 0;
/* SELECT SIX HIGHEST PROB ELEMENT STATE NODES: */
for (k = 0; k < 6; ++k) {
pmax = 0.f;
for (i = 0; i < 5; ++i) {
n = 6*i+k;
for (ip = 0; ip < *isave; ++ip) {
j = (ip ) * 30 + (i) * 6 + k;
#ifdef DEBUG
//printf("\nSAVEP:p[%d][%d]:%f pmax:%f jsav:%d", n,ip,p[n][ip],pmax,jsav);
#endif
if (p[n][ip] >= pmax) {
pmax = p[n][ip];
jsav = j;
savn = n;
ipsav = ip;
}
}
}
if (pmax >= 0) { //1e-6f) {
psum += pmax;
psav[nsav] = pmax;
nsave[nsav] = savn;
pathsv[nsav] = ipsav;
sort[nsav] = jsav;
#ifdef DEBUG
printf("\nSAVEP: psav[%d]:%f nsave[%d]:%d pathsv:%d sort[%d]:%d pmax:%f k:%d i:%d j:%d ip:%d", nsav,psav[nsav],nsav,nsave[nsav], pathsv[nsav],nsav, sort[nsav],pmax,k,i,j,ip);
#endif
++nsav;
}
}
/* SELECT ENOUGH ADDITIONAL NODES TO MAKE TOTAL */
/* PROBABILITY SAVED EQUAL TO POPT, OR A MAX OF 25: */
do {
pmax = 0.f;
for (ip = 0; ip < *isave; ++ip) {
for (n = 0; n < 30; ++n) {
j = ip*30 + n;
for (i = 0; i < nsav; ++i) {
if (j == sort[i]) {
goto L500;
}
}
if (p[n][ip] > pmax) {
pmax = p[n][ip];
jsav = j;
savn = n;
ipsav = ip;
}
L500:
;
}
}
psum += pmax;
psav[nsav] = pmax;
nsave[nsav] = savn;
pathsv[nsav] = ipsav;
sort[nsav] = jsav;
#ifdef DEBUG
printf("\nSAVEP: psum:%f psav[%d]:%f nsave[%d]:%d pathsv[%d]:%d sort[%d]:%d", psum, nsav,psav[nsav], nsav,nsave[nsav],nsav, pathsv[nsav], nsav, sort[nsav]);
#endif
++nsav;
if (psum >= popt) {
break; //goto L600;
}
}
while (nsav < 25);
L600:
/* NEW ISAVE EQUALS NO. OF NODES SAVED: */
#ifdef DEBUG
//printf("\nNEW ISAVE EQUALS # OF NODES SAVED:%d", nsav);
#endif
// *isave = nsav;
/* SORT THE SAVED ARRAYS TO OBTAIN THE ARRAYS */
/* TO BE USED FOR THE NEXT ITERATION: */
/* ALSO OBTAIN HIGHEST PROBAILITY NODE: */
for (i = 0; i < *isave; ++i) {
// p[i] = psav[i] / psum;
lambda[i] = lamsav[nsave[i]][pathsv[i]]; //[sort[i]];
dur[i] = dursav[nsave[i]][pathsv[i]]; //[sort[i]];
ilrate[i] = ilrsav[nsave[i]][pathsv[i]]; //[sort[i]];
blksv.ykkip[i] = blksv.ykksv[sort[i]];
blksv.pkkip[i] = blksv.pkksv[sort[i]];
#ifdef DEBUG
//printf("\nSAVEP: p[%2d]:%f psav[%2d]:%f lambda[%2d]:%d dur[%2d]:%f ilrate[%2d]:%d ykkip[%2d]:%f isave:%d psum:%f", i, p[i], i ,psav[i], i, lambda[i],i,dur[i], i, ilrate[i], i, blksv.ykkip[i],*isave,psum);
#endif
}
for (i = 0; i < *isave; ++i) {
iconv[i] = 1;
}
// compare: if speedrate, duration and letter state has not changed then set iconv[k] to zero
for (n = 0; n < (*isave-1); ++n) {
if (iconv[n] == 0) goto L800;
nplus1 = n + 1;
for (k = nplus1; k < *isave; ++k) {
if (iconv[k] == 0) goto L810;
if (ilrate[k] != ilrate[n]) goto L810;
if (dur[k] != dur[n]) goto L810;
if (lambda[k] != lambda[n]) goto L810;
#ifdef DEBUG
//printf("\nSAVEP: ICONV[k:%d]=0 n:%d",k,n);
#endif
iconv[k] = 0;
L810:
;
}
L800:
;
}
/*
psum = 0.f;
n = 0;
for (i = 0; i < *isave; ++i) {
if (iconv[i] != 0) {
++n;
lambda[n] = lambda[i];
dur[n] = dur[i];
ilrate[n] = ilrate[i];
blksv.ykkip[n] = blksv.ykkip[i];
blksv.pkkip[n] = blksv.pkkip[i];
pathsv[n] = pathsv[i];
sort[n] = sort[i];
p[n] = p[i];
psum += p[n];
}
}
//printf("\nSAVEP: NEW ISAVE %d", n);
// *isave = n;
pmax = 0.f;
for (i = 1; i < *isave; ++i) {
p[i] /= psum;
if (p[i] > pmax) {
pmax = p[i];
*imax = i;
}
}
*/
return 0;
} /* savep_ */

78
tmpcsave/simsgi.c
Wyświetl plik

@ -1,78 +0,0 @@
/* simsgi.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
/* Table of constant values */
static integer c__1 = 1;
static real c_b3 = 0.f;
static real c_b6 = 6.28319f;
/* Subroutine */ int simsgi_(real *x, real *sig)
{
/* Initialized data */
static real xlast = 1.f;
static real beta = 1.f;
static real amp = 1.f;
static real bfade = 0.f;
static real theta = 1.f;
static real phi = 1.f;
/* Builtin functions */
double exp(doublereal), r_mod(real *, real *), sin(doublereal);
/* Local variables */
static real w;
static integer nc;
static real tk, zn;
extern /* Subroutine */ int key_(real *, real *);
static real dur, ampb, tdur;
extern double randn(real mean, real sigma);
static real wchrp, phisgm;
dur = beta;
key_(&dur, x);
beta = beta * (1.f - *x - xlast + *x * 2.f * xlast) + 1.f;
tk = *x * (1.f - xlast);
xlast = *x;
w = randn(c_b3, blk2.asigma);
amp += tk * w;
if (amp < .01f) {
amp = .01f;
}
w = randn(c_b3, blk2.bsigma);
bfade = blk2.gamma * bfade + w;
ampb = amp + bfade;
if (ampb < .001f) {
bfade = .001f - amp;
}
ampb = amp + bfade;
tdur = blk1.tau * 1e3f * beta;
wchrp = *x * blk2.wchirp * exp(-tdur / blk2.tchirp);
theta += (nc + wchrp) * blk1.tau;
theta = r_mod(&theta, &c_b6);
w = randn(c_b3, phisgm);
phi += tk * w;
phi = r_mod(&phi, &c_b6);
*sig = *x * ampb * sin(theta + phi);
zn = randn(c_b3, blk2.rsigm);
return 0;
} /* simsgi_ */

80
tmpcsave/spdtr.c
Wyświetl plik

@ -1,80 +0,0 @@
/* spdtr.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
/* Common Block Declarations */
doublereal spdtr_(integer *isrt, integer *ilrt, integer *iselm, integer *ilelm)
{
/* System generated locals */
real ret_val;
/* Local variables */
static integer ind1, idel, idelsp, israte;
/* THIS FUNCTION RETURNS THE DATA RATE (SPEED) TRANSITION */
/* PROBABILITY BASED ON THE CURRENT ELEM TYPE. THE ALLOWABLE */
/* TRANSITION PROBS ARE STORED IN THE TABLE RTRANS. */
/* VARIABLES: */
/* ISRT - DATA RATE IDENTITY FOR STATE TO WHICH PATH IS BEING EXTENDED */
/* ILRT - DATA RATE ON CURRENT PATH */
/* ISELM - ELEM TYPE FOR NEXT STATE */
/* ILELM - ELEM TYPE ON CURRENT PATH */
/* IF SAVED ELEMENT AND NEW ELEMENT ARE THE SAME THEN THERE CAN BE NO SPEED CHANGE: */
if (*ilelm != *iselm) {
goto L100;
}
ret_val = 1.f;
if (*isrt != 2) {
ret_val = 0.f;
}
goto L300;
/* OTHERWISE, OBTAIN SPEED TRANSITION PROB */
L100:
idel = blkrat.memdel[*iselm][*ilelm] ;
ind1 = blkspd.mempr[*iselm][*ilelm];
if (ind1 != 0) {
goto L200;
}
ret_val = 0.f;
goto L300;
L200:
idelsp = (*isrt - 2) * idel;
ret_val = blkspd.rtrans[ind1-1][*isrt];
israte = *ilrt + idelsp;
if (israte > 60) {
ret_val = 0.f;
}
if (israte < 10) {
ret_val = 0.f;
}
L300:
#ifdef DEBUG
printf("\nSPDTR: isrt:%d ilrt:%d idel:%d=memdel[%d][%d] ind1:%d idelsp:%d israte:%d retval:%f, rtrans[%d][%d]:%f", *isrt,*ilrt, idel, *iselm, *ilelm, ind1,idelsp, israte,ret_val,ind1-1,*isrt,blkspd.rtrans[ind1-1][*isrt]);
#endif
return ret_val;
} /* spdtr_ */

74
tmpcsave/sprob.c
Wyświetl plik

@ -1,74 +0,0 @@
/* sprob.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
int sprob_(real p[][25], integer *isave, integer ilrsav[][25], real *pelm, integer *khat, real *spdhat, real *px)
{
/* Local variables */
static integer i, j, k, m, n;
static real pselem[6];
/* SPROB COMPUTES THE POSTERIOR PROBS OF THE ELEMENT */
/* STATES, DATA RATE STATES, AND KEYSTATES BY SUMMING */
/* OVER THE APPROPRIATE PATHS. */
/* VARIABLE: */
/* P INPUT PATH PROBABILITIES */
/* ISAVE NUMBER OF PATHS SAVED */
/* PSELEM OUTPUT ELEMENT PROB */
/* SPDHAT OUTPUT SPEED ESTIMATE (DATA RATE WPM) */
/* PX OUTPUT KEYSTATE PROBABILITY */
/* Function Body */
*spdhat = 0.f;
*px = 0.f;
/* FOR EACH STATE EXTENSION OF PATH M: */
/* OBTAIN ELEMENT STATE PROBS,KEYSTATE PROBS,SPEED EST: */
for (k = 0; k < 6; ++k) {
pselem[k] = 0.f;
for (i = 0; i < 5; ++i) {
n = (i*6) + k;
for (m = 0; m < *isave; ++m) {
j = (m*30) + n;
pselem[k] += p[n][m];
*spdhat += ilrsav[n][m] * p[n][m];
#ifdef DEBUG
printf("\nSPROB: pselem[%d]:%f spdhat:%f ilrsav[%3d][%3d]:%d \t p[%3d][%3d]:%f \tpx:%f", k, pselem[k], *spdhat, n,m,ilrsav[n][m], n,m, p[n][m],*px);
#endif
if (k < 2) { // KEYSTATE DOT (0) OR DASH (1)
*px += p[n][m];
}
L100:
;
}
}
}
*pelm = 0.f;
for (k = 0; k < 6; ++k) {
if (pselem[k] >= *pelm) {
*pelm = pselem[k ];
*khat = k;
}
}
return 0;
} /* sprob_ */

31
tmpcsave/stats.c
Wyświetl plik

@ -1,31 +0,0 @@
/* stats.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Subroutine */ int stats_(real *xin1, real *xin2, real *xin3, real *xin4,
real *s1, real *s2, real *s3, real *s4, integer *n)
{
/* Parameter adjustments */
--s4;
--s3;
--s2;
--s1;
/* Function Body */
s1[*n] = *xin1;
s2[*n] = *xin2;
s3[*n] = *xin3;
s4[*n] = *xin4;
return 0;
} /* stats_ */

255
tmpcsave/transl.c
Wyświetl plik

@ -1,255 +0,0 @@
/* transl.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
/* Table of constant values */
static integer c__1 = 1;
/* Subroutine */ int transl_(integer *ltr)
{
/* Initialized data */
static char ispace[1] = " ";
static integer spflag = 0;
static integer nchar = 0;
/* Format strings */
static char fmt_5000[] = "(1x,i1,s)";
static char fmt_1500[] = "(2x,a1,/)";
static char fmt_1000[] = "(1x,a1,s)";
static char fmt_1100[] = "(1x,a2,$)";
static char fmt_1200[] = "(2x,a2,2x,$)";
static char fmt_1300[] = "(2x,a3,2x,$)";
static char fmt_1400[] = "(/,10x,a2,/,10x)";
static char fmt_1700[] = "(2x,a3,2x,$)";
static char fmt_1600[] = "(/,10x)";
/* System generated locals */
integer i__1;
/* Builtin functions */
integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
/* Local variables */
static integer i__;
static real o;
static integer ixl, nblank, elmhat, ltrhat, lstelm, ixlast, elmout;
static real txlast;
static integer lstltr, ltrout;
/* Fortran I/O blocks */
static cilist io___13 = { 0, 6, 0, fmt_5000, 0 };
static cilist io___15 = { 0, 6, 0, fmt_1500, 0 };
static cilist io___16 = { 0, 6, 0, fmt_1000, 0 };
static cilist io___18 = { 0, 6, 0, fmt_1000, 0 };
static cilist io___19 = { 0, 6, 0, fmt_1100, 0 };
static cilist io___20 = { 0, 6, 0, fmt_1000, 0 };
static cilist io___21 = { 0, 6, 0, fmt_1200, 0 };
static cilist io___22 = { 0, 6, 0, fmt_1000, 0 };
static cilist io___23 = { 0, 6, 0, fmt_1300, 0 };
static cilist io___24 = { 0, 6, 0, fmt_1400, 0 };
static cilist io___25 = { 0, 6, 0, fmt_1700, 0 };
static cilist io___26 = { 0, 6, 0, fmt_1000, 0 };
static cilist io___27 = { 0, 6, 0, fmt_1600, 0 };
/* THIS SUBROUTINE PRODUCES THE OUTPUT TEXT ON A CRT. */
/* IT USES THE SIMPLE FORMATTING RULES DESCRIBED IN THE */
/* TEXT. */
/* DETERMINE IF A CSP,WSP, OR PAUSE TO MARK TRASITION */
/* HAS OCCURED; IF SO LTR IS READY FOR OUTPUT: */
// elmhat = blklam.ilami[(1600 + (0 + (blklam.ielmst[(0 + (0 + (*ltr - 1<< 2))) / 4] - 1 << 2)) - 1600) / 4];
elmhat = blklam.ilami[blklam.ielmst[*ltr-1]-1];
ixl = blklam.ilamx[elmhat - 1];
#ifdef DEBUG
printf("\nTRANSL: elmhat:%d ltr:%d ixl:%d lstltr:%d", elmhat, *ltr, ixl, lstltr);
#endif
if (ixl == ixlast) {
goto L700;
}
if (ixl == 1 && lstelm >= 4) {
goto L10;
}
if (ixl == 0 && lstelm <= 2) {
goto L700;
}
goto L700;
L10:
ltrhat = lstltr;
ltrout = alpha.ialph[blktrn.ltrmap[ltrhat - 1] - 1];
nblank = blktrn.iblank[ltrhat - 1];
elmout = blklam.ilami[blklam.ielmst[ltrhat - 1] - 1];
goto L40;
s_wsfe(&io___13);
do_fio(&c__1, (char *)&elmout, (ftnlen)sizeof(integer));
e_wsfe();
++nchar;
L40:
if (blktrn.ltrmap[ltrhat - 1] == 43) {
goto L50;
}
if (blktrn.ltrmap[ltrhat - 1] <= 44) {
goto L100;
}
if (blktrn.ltrmap[ltrhat - 1] <= 46) {
goto L200;
}
if (blktrn.ltrmap[ltrhat - 1] <= 60) {
goto L300;
}
if (blktrn.ltrmap[ltrhat - 1] == 61) {
goto L400;
}
if (blktrn.ltrmap[ltrhat - 1] == 66) {
goto L500;
}
goto L550;
L50:
if ((real) spflag == o) {
goto L100;
}
nchar = 0;
s_wsfe(&io___15);
do_fio(&c__1, (char *)&ltrout, (ftnlen)sizeof(integer));
e_wsfe();
spflag = 1;
goto L600;
L100:
#ifdef DEBUG
printf("\nTRANSL: one char ltrhat:%d ltrmap:%d ltrout:%c", ltrhat, blktrn.ltrmap[ltrhat - 1],alpha.ialph[blktrn.ltrmap[ltrhat - 1] - 1]);
#endif
++nchar;
s_wsfe(&io___16);
do_fio(&c__1, (char *)&ltrout, (ftnlen)sizeof(integer));
e_wsfe();
spflag = 0;
if (nblank == 0) {
goto L110;
}
spflag = 1;
i__1 = nblank;
for (i__ = 1; i__ <= i__1; ++i__) {
++nchar;
s_wsfe(&io___18);
do_fio(&c__1, ispace, (ftnlen)1);
e_wsfe();
L110:
;
}
goto L600;
L200:
#ifdef DEBUG
printf("\nTRANSL: 2_ char:");
#endif
nchar += 2;
s_wsfe(&io___19);
do_fio(&c__1, (char *)&ltrout, (ftnlen)sizeof(integer));
e_wsfe();
spflag = 0;
if (nblank == 0) {
goto L210;
}
spflag = 1;
i__1 = nblank;
for (i__ = 1; i__ <= i__1; ++i__) {
++nchar;
s_wsfe(&io___20);
do_fio(&c__1, ispace, (ftnlen)1);
e_wsfe();
L210:
;
}
goto L600;
L300:
#ifdef DEBUG
printf("\nTRANSL: 4_ char:");
#endif
nchar += 4;
s_wsfe(&io___21);
do_fio(&c__1, (char *)&ltrout, (ftnlen)sizeof(integer));
e_wsfe();
spflag = 1;
if (nblank == 0) {
goto L310;
}
i__1 = nblank;
for (i__ = 1; i__ <= i__1; ++i__) {
++nchar;
s_wsfe(&io___22);
do_fio(&c__1, ispace, (ftnlen)1);
e_wsfe();
L310:
;
}
goto L600;
L400:
#ifdef DEBUG
printf("\nTRANSL: 5 char:");
#endif
nchar += 5;
s_wsfe(&io___23);
do_fio(&c__1, (char *)&ltrout, (ftnlen)sizeof(integer));
e_wsfe();
spflag = 1;
goto L600;
L500:
nchar = 0;
s_wsfe(&io___24);
do_fio(&c__1, (char *)&ltrout, (ftnlen)sizeof(integer));
e_wsfe();
spflag = 1;
goto L600;
L550:
#ifdef DEBUG
printf("\nTRANSL: 5_ char:");
#endif
nchar += 5;
s_wsfe(&io___25);
do_fio(&c__1, (char *)&ltrout, (ftnlen)sizeof(integer));
e_wsfe();
spflag = 0;
if (nblank == 0) {
goto L560;
}
spflag = 1;
i__1 = nblank;
for (i__ = 1; i__ <= i__1; ++i__) {
++nchar;
s_wsfe(&io___26);
do_fio(&c__1, ispace, (ftnlen)1);
e_wsfe();
L560:
;
}
L600:
if (nchar < 52) {
goto L700;
}
s_wsfe(&io___27);
e_wsfe();
nchar = 0;
L700:
txlast = (real) ixl;
lstelm = elmhat;
lstltr = *ltr;
return 0;
} /* transl_ */

228
tmpcsave/trelis.c
Wyświetl plik

@ -1,228 +0,0 @@
/* trelis.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
#define DEBUG 1
/* Common Block Declarations */
/* Table of constant values */
static integer c__1 = 1;
/* Subroutine */ int trelis_(integer *isave, integer *pathsv, integer *lambda, integer *imax)
{
/* Initialized data */
static int init = 0;
static real pthtrl[200][25];
static integer lmdsav[200][25];
static integer n = 0;
static integer ndelay = 200;
static integer ipnod[25] = { 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1 };
static integer ncall = 0;
static integer nmax = 0;
static integer mmax = 0;
/* Local variables */
static integer i, i1, k, k0, kd, ip, ieq, ltr, ndel;
static real wait;
static integer isavg, ipmax;
static real xsavg, xmmax, xnmax;
static integer ltrsv[200], ndlavg;
static real xdlavg;
static integer ndelst;
extern /* Subroutine */ int transl_(integer *);
if (init == 0) {
init = 1;
for (i = 0; i < 200; i++) {
for (k=0;k<25;k++){
pthtrl[i][k] = 5.f;
lmdsav[i][k] = 5;
}
}
}
/* THIS SUBROUTINE STORES THE SAVED NODES AT EACH */
/* STAGE AND FORMS THE TREE OF SAVED PATHS LINKING */
/* THE NODES. DECODING IS ACCOMPLISHED BY FINDING */
/* THE CONVERGENT PATH IF IT OCCURS WITHIN A MAXIMUM */
/* DELAY SET BY THE PARAMETER NDELAY. IF CONVERGENCE */
/* TO A SINGLE PATH DOES NOT OCCURS, THEN DECODING IS */
/* DONE BY READING THE LETTER ON THE PATH WITH HIGHEST */
/* PROBABILITY */
/* Parameter adjustments */
// --lambda;
// --pathsv;
/* Function Body */
/* KEEP AVERAGE OF ISAVE, NDEL FOR DATA ANALYSIS: */
++ncall;
if (blkend.iend != 1) {
goto L10;
}
isavg = xsavg;
ndlavg = xdlavg;
blkend.iend = 0;
#ifdef DEBUG
printf("\nAVG NO OF PATHS SAVED: %d %d", isavg, ndlavg);
printf("\nPERCENT OF TIME PATHS = 25: %f %f", xmmax, xnmax);
printf("\nPERCENT OF TIME DELAY = 200:");
#endif
L10:
xsavg = (xsavg * (ncall - 1) + *isave) / ncall;
xdlavg = (xdlavg * (ncall - 1) + ndel) / ncall;
if (ndel != ndelay) {
goto L20;
}
++nmax;
xnmax = (real) nmax;
xnmax /= ncall;
L20:
if (*isave != 25) {
goto L30;
}
++mmax;
xmmax = (real) mmax;
xmmax /= ncall;
L30:
/* STORE PATHSV AND CORRESPONDING LAMBDA IN THE TRELLIS USING A CIRCULAR BUFFER OF LENGTH NDELAY : */
++n;
if (n == ndelay + 1) {
n = 1;
}
i1 = *isave;
for (i = 0; i <*isave; ++i) {
pthtrl[n][i] = (real) pathsv[i];
lmdsav[n][i] = lambda[i];
#ifdef DEBUG
printf("\nTRELIS: n:%d pthtrl[%d][%d]:%f lmdsav[%d][%d]:%d",n, n,i,pthtrl[n][i], n,i,lmdsav[n][i]);
#endif
}
/* PERFORM DYNAMIC PROGRAM ROUTINE TO FIND CONVERGENT PATH: */
k = 0;
for (i = 0; i < *isave; ++i) {
ipnod[i] = i;
}
L190:
++k;
if (k == ndelay) {
goto L700;
}
for (ip = 0; ip < *isave; ++ip) {
i = n - k + 1;
if (i <= 0) {
i = ndelay + i;
}
ipnod[ip] = pthtrl[ipnod[ip]][i];
if (ip == *imax) {
ipmax = ipnod[ip];
}
}
/* IF ALL NODES ARE EQUAL,THEN PATHS CONVERGE: */
for (ieq = 1; ieq <= *isave; ++ieq) {
if (ipnod[0] != ipnod[ieq]) {
goto L190;
}
}
/* PATHS CONVERGE; SET NDEL: */
ndel = k + 1;
/* IF POINT OF CONVERGENCE IS SAME AS IT WAS ON */
/* LAST CALL, THEN NO NEED TO RE-DECODE SAME NODE: */
if (ndel == ndelst + 1) {
goto L800;
}
/* IF POINT OF CONVERGENCE OCCURS AT SAME DELAY AS LAST CALL, THEN TRANSLATE: */
if (ndel != ndelst) {
goto L350;
}
i = n - ndel + 1;
if (i <= 0) {
i = ndelay + 1;
}
ltr = lmdsav[ipnod[0]][i];
#ifdef DEBUG
printf("\nTRELIS: IF POINT OF CONVERGENCE OCCURS AT SAME DELAY AS LAST CALL, THEN TRANSLATE:");
#endif
transl_(&ltr);
goto L800;
/* OTHERWISE,POINT OF CONVERGENCE HAS OCCURED */
/* EARLIER ON THIS CALL, SO NEED TO TRANSLATE */
/* EVERYTHING ON THE CONVERGENT PATH FROM */
/* PREVIOUS POINT OF CONVERGENCE TO THIS POINT: */
L350:
k0 = 0;
ip = ipnod[0];
i1 = ndelst;
for (k = ndel; k <= i1; ++k) {
++kd;
i = n - k + 1;
if (i <= 0) {
i = ndelay + i;
}
ltrsv[kd - 1] = lmdsav[ip][i];
ip = pthtrl[ip][i];
}
/* REVERSE ORDER OF DECODED LETTERS, SINCE THEY */
/* WERE OBTAINED FROM THE TRELLIS IN REVERSE; */
/* TRANSLATE EACH: */
i1 = kd;
for (i = 1; i <= i1; ++i) {
ltr = ltrsv[kd - i];
#ifdef DEBUG
printf("\nTRELIS: REVERSE ORDER OF DECODED LETTERS - TRANSLATE EACH:");
#endif
transl_(&ltr);
/* L500: */
}
goto L800;
L700:
/* PATHS HAVE NOT CONVERGED AT MAXIMUM ALLOWABLE */
/* DELAY, SO TRANSLATE WHAT IS ON HIGHEST */
/* PROBABILITY PATH: */
ndel = ndelay;
i = n - ndelay + 1;
if (i <= 0) {
i = ndelay + i;
}
ltr = lmdsav[ipmax][i];
#ifdef DEBUG
printf("\nTRELIS: PATHS HAVE NOT CONVERGED AT MAXIMUM ALLOWABLE DELAY SO TRANSLATE:");
#endif
transl_(&ltr);
/* PRUNE AWAY NODES WHICH ARE NOT ON THIS PATH: */
for (k = 0; k < *isave; ++k) {
if (ipnod[k] == ipmax) {
goto L750;
}
lambda[k] = 0;
L750:
;
}
L800:
ndelst = ndel;
return 0;
} /* trelis_ */

84
tmpcsave/trprob.c
Wyświetl plik

@ -1,84 +0,0 @@
/* trprob.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
#include <stdio.h>
//#define DEBUG 1
/* Common Block Declarations */
int trprob_(integer ip, integer *lambda, real *dur, integer *ilrate, real p[][25])
{
static integer i, k, n;
static real pin[30];
static integer kelm;
static real psum, ptrx;
static integer ielem, irate;
extern int ptrans_(integer *, integer *, integer *, integer *, real *, real *, real *, integer *);
extern doublereal xtrans_(integer *, real , integer *);
/* THIS SUBROUTINE COMPUTES THE TRANSITION PROBABILITY */
/* FROM SAVED PATH IP TO EACH STATE N AND STORES THE */
/* RESULT IN P(IP, N). */
/* VARIABLES: */
/* IP - INPUT SAVED PATH IDENTITY */
/* LAMBDA - INPUT SAVED LTR STATE IDENTITY */
/* DUR - INPUT SAVED ELEMENT DURATION */
/* ILPATE - INPUT SAVED DATA RATE IDENTITY */
/* P - OUTPUT TRANSITION PROBABILITY MATRIX */
/* THE FOLLOWING FUNCTION SUBROUTINES ARE USED: */
/* XTRANS - RETURNS THE KEYSTATE TRANSITION PROBABILITY */
/* CONDITIONED ON ELEMENT TYPE AND DATA RATE */
/* PTRANS - RETURNS THE PATH-CONDITIONAL STATE TRANSITION PROB */
/* LOOK UP ELEMENT TYPE FOR LTR STATE LAMBDA: */
if (*lambda == 0) {
for (n = 0; n < 30; ++n) {
p[n][ip] = 0.f;
}
#ifdef DEBUG
printf("\nTRPROB: lambda:%d", *lambda);
#endif
return 0;
}
ielem = blklam.ilami[blklam.ielmst[*lambda - 1] - 1];
/* COMPUTE KEYSTATE TRANSITION PROBABILITY: */
ptrx = xtrans_(&ielem, *dur, ilrate);
/* FOR EACH STATE, COMPUTE STATE TRANSITION PROBABILITY: */
psum = 0.f;
for (k = 0; k < 6; ++k) { // ELEMENT STATE
for (i = 0; i < 5; ++i) { // DATA RATE STATE
n = i*6 + k;
ptrans_(&k, &i, lambda, ilrate, &ptrx, &psum, pin, &n);
}
}
for (n = 0; n < 30; ++n) {
p[n][ip] = pin[n] / psum;
#ifdef DEBUG
printf("\nTRPROB: p[%2d][%2d]:%f pin[%2d]:%f psum:%f lambda:%d", n,ip, p[n][ip], n, pin[n], psum, *lambda);
#endif
}
L200:
return 0;
} /* trprob_ */

84
tmpcsave/xtrans.c
Wyświetl plik

@ -1,84 +0,0 @@
/* xtrans.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "common.h"
doublereal xtrans_(integer *ielem, real d0, integer *irate)
{
/* Initialized data */
static integer kimap[6] = { 1,3,1,3,7,14 };
static real aparm[3] = { 3.f,1.5f,1.f };
/* System generated locals */
real ret_val;
/* Builtin functions */
double exp(doublereal);
/* Local variables */
static real b0, b1, p0, p1, alpha;
static integer mscale;
static real rscale;
/* THIS FUNCTION IMPLEMENTS THE CALCULATION OF KEYSTATE */
/* TRANSITION PROBABILITY, CONDITIONED ON ELEMENT TYPE, */
/* CURRENT DURATION, AND DATA RATE. */
/* VARIABLES: */
/* IELEM- INPUT CURRENT ELEMENT TYPE */
/* D0- INPUT CURRENT ELEMENT DURATION */
/* IRATE - INPUT CURRENT DATA RATE */
/* TABLES IN COMMON CONTAIN DENSITY PARMS FOR EACH ELEMENT TYPE, DATA RATE. */
/* SCALE DURATION AND OBTAIN DENSITY PARAMETER: */
// mscale = kimap[(0 + (0 + (*ielem - 1 << 2))) / 4];
mscale = kimap[*ielem];
rscale = 1200.f / *irate;
b0 = d0 / (mscale * rscale);
b1 = (d0 + 5.f) / (mscale * rscale);
if (*ielem == 5) {
goto L20;
}
if (*ielem == 4) {
goto L10;
}
alpha = mscale * aparm[0];
goto L100;
L10:
alpha = aparm[1] * 7.f;
goto L100;
L20:
alpha = aparm[2] * 14.f;
L100:
if (b1 <= 1.f) {
goto L200;
}
if (b0 < 1.f && b1 > 1.f) {
goto L300;
}
ret_val = exp(-alpha * (b1 - b0));
goto L400;
L200:
p1 = 1.f - exp(alpha * (b1 - 1.f)) * .5f;
p0 = 1.f - exp(alpha * (b0 - 1.f)) * .5f;
ret_val = p1 / p0;
goto L400;
L300:
p1 = exp(-alpha * (b1 - 1.f)) * -.5f;
p0 = 1.f - exp(alpha * (b0 - 1.f)) * .5f;
ret_val = p1 / p0;
L400:
return ret_val;
} /* xtrans_ */

23
tmpsave/autocr.f
Wyświetl plik

@ -1,23 +0,0 @@
SUBROUTINE AUTOCR (S5, RS)
DIMENSION S5(512), RS(512), S(1000), RS1(500)
DATA S/1000*0./, XN/0./
XN=XN+1
DO 100 I=1, 500
S(I) = S5(I)
RS1(I)=0.
100 CONTINUE
DO 200 I=1, 500
DO 300 K=1, 500
RS1(I)= RS1(I)+S(K+I-1)*S(K)
300 CONTINUE
200 CONTINUE
DO 400 I=1, 500
RS(I)=(RS(I)*(XN-1.)+RS1(I))/XN
400 CONTINUE
RETURN
END

49
tmpsave/bpfdet.f
Wyświetl plik

@ -1,49 +0,0 @@
SUBROUTINE BPFDET(ZIN, Z)
C
C THIS SUBROUTINE IMPLEMENTS THE BANDPASS FILTER AND
C ENVELOPE DETECTOR FUNCTIONS. THE BPF IS A SIMPLE CASCADE
C OF TWO 2-POLE DIGITAL RESONATORS AT A CENTER FREQ OF
C 1000 HZ. THE LPF OF THE ENVELOPE DETECTOR IS A
C THREE - POLE CHEBYSCHEV 100 HZ LPF.
C
C
DIMENSION A ( 4 )
DATA A/ .000030051, 2.9507982, 2.90396345, -0.953135172/
DATA CK1/1.37158/ , CK2/ 0.9409/, CG/ 0.0150/
DATA C1 / 1.2726/ ,C2/ 0.8100/ , C / 0.1900/
C OFF IS TWO 2 - POLE RESONATORS :
W3 = W2
W2 = W1
W1=C1*W2-C2*W3+C*ZIN
X3=X2
X2=X1
X1=CK1*X2-CK2*X3+CG*W1
ZBPF=X1
C ENVELOPE DETECTOR (SQUARE-LAW):
C sQUARE-
XDET=SQRT(ZBPF**2)
C LOW-PASS FILTER-
Y3=Y2
Y2=Y1
Y1=Y0
Y0=XDET*A(1)
Z3=Z2
Z2=Z1
Z1=Z
Z=Y0+3.*(Y1+Y2)+Y3
Z=Z+A(2)*Z1-A(3)*Z2-A(4)*23
RETURN
END

90
tmpsave/initl.f
Wyświetl plik

@ -1,90 +0,0 @@
SUBROUTINE INITL
DIMENSION IELMST(400),ILAMI(16),ILAMX(6)
DIMENSION ELEMTR(16,6),RTRANS(5,2),ISX(6)
DIMENSION MEMFCN(400,6),LTRMAP(400)
DIMENSION MEMDEL(6,6),MEMPR(6,6),IBLANK(400)
DIMENSION IARRAY(8)
DIMENSION ITEXT(200)
CHARACTER*70 IALPH
COMMON /BLKLAM/ IELMST,ILAMI,ILAMX
COMMON /BLKRAT/ MEMDEL
COMMON /BLKELM/ ELEMTR
COMMON /BLKSPD/ RTRANS,MEMPR
COMMON /BLKMEM/ MEMFCN
COMMON /BLKS/ ISX
COMMON /BLKTRN/LTRMAP,IBLANK, IALPH
common /BLKTXT/ ITEXT
DATA ISX/1,1,0,0,0,0/
DATA MEMFCN /9,11,13,15,9,11,13,15,9,0,11,0,13,0,15,0,384*0,
&10,12,14,16,10,12,14,16,0,10,0,12,0,14,0,16,384*0,
&1,0,0,0,5,0,0,0,1,5,1,5,1,5,1,5,384*0,
&0,2,0,0,0,6,0,0,2,6,2,6,2,6,2,6,384*0,
&0,0,3,0,0,0,7,0,3,7,3,7,3,7,3,7,384*0,
&0,0,0,4,0,0,0,8,4,8,4,8,4,8,4,8,384*0/
DATA IELMST/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,384*0/
DATA ILAMI/3,4,5,6,3,4,5,6,1,2,1,2,1,2,1,2/
DATA ILAMX/1,1,0,0,0,0/
DATA LTRMAP/3,4,5,6,3,4,5,6,1,2,1,2,1,2,1,2,384*0/
DATA IALPH /'ABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890;:%&'/
C /'A','B','C','D','E','F','G','H','I','J','K','L','M'/
C 2 'N','O','P','Q','R','S','T','U','V','W','X','Y','Z',
C 2 '1','2','3','4','5','6','7'/
C 2 '8','9','0',';',':','%','&',0,0,'^',',','AS','SN',
C 2 0,0,0,0,'NR','NO','GA','OK','AR','SK',0,0,0,0,
C 2 'IMI',0,0,0,0,'BT',0,0,0,'EEE'/
DATA IBLANK/400*0/
DATA ELEMTR/.55,.6,.5,.5,.56,.5,.5,.5,8*0.,
&.45,.5,.5,.5,.45,.5,.5,.5,8*0.,
&8*0.,.581,.54,.923,.923,.923,.93,.92,.95,
&8*0.,.335,.376,.062,.062,.062,.062,.04,.04,
&8*0.,.067,.069,.012,.012,.012,.012,.009,.009,
&8*0.,.015,.015,.003,.003,.003,.003,.001,.001/
DATA RTRANS /.1,.2,.4,.2,.1,.15,.2,.3,.2,.15/
DATA MEMDEL /0,0,2,2,5,10,0,0,2,2,5,10,
&2,2,0,0,0,0,2,2,0,0,0,0,2,2,0,0,0,0,
&2,2,0,0,0,0/
DATA MEMPR /0,0,1,2,1,2,0,0,1,2,1,2,1,1,0,0,0,0,
&1,1,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0/
OPEN (UNIT=20, FILE='MORSEM')
DO 10 I=1,300
READ (20,30) (IARRAY(K),K=1,8)
30 FORMAT(8I3)
DO 11 K=1,6
11 MEMFCN(I,K) = IARRAY(K+2)
LTRMAP(I) = IARRAY(1)
IELMST(I) = IARRAY(2)
IF ((IELMST(I).EQ.7).OR.(IELMST(I).EQ.3)) IBLANK(I) = 1
IF ((IELMST(I).EQ.8).OR.(IELMST(I).EQ.4)) IBLANK(I) = 2
10 CONTINUE
ENDFILE 20
OPEN(UNIT=20,FILE='OUTPUT')
DO 50 I=1,300
WRITE(20,40) (MEMFCN(I,K),K=1,6)
40 FORMAT(10X,6(I3,2X))
50 CONTINUE
ENDFILE 20
OPEN(UNIT=20,FILE='TEXT')
DO 60 I=1,108
READ(20,70) ITEXT(I)
70 FORMAT(I2)
60 CONTINUE
ENDFILE 20
RETURN
END

67
tmpsave/inputl.f
Wyświetl plik

@ -1,67 +0,0 @@
SUBROUTINE INPUTL
DIMENSION ESEP(6),EDEV(6)
COMMON /BLK1/TAU/BLK6/ DMEAN,XDUR,ESEP,EDEV
COMMON /BLK2/ WC, WCHIRP, ASIGMA, BSIGMA, PHISGM, RSIGM, TCHIRP,
&GAMMA
DATA TAU/.000125/
DATA ESEP/1,3,1,3,7,14/
DATA EDEV/6*0./
DATA XDUR/0./
PRINT 100
100 FORMAT(1X,'INPUT KEYING PARMS: RATE,MEAN ELEM DURATIONS')
C READ 200,RATE,(ESEP(K), K=1,6)
PRINT 150
150 FORMAT(1X,'INPUT ELEM DURATION STD OEVIATIONS')
C READ 200,(EDEV(K),K=1,6)
200 FORMAT(F7.1)
PRINT 300
300 FORMAT(1X,'INPUT SIG PARMS- AVAR,BVAR,FCHIRP,TCHIRP,PHIVAR')
C READ 210,AVAR,BVAR,FCHIRP,TCHIRP,PHIVAR
210 FORMAT(F7.1,F7.1,F7.1,F7.1,F7.1)
PRINT 400
400 FORMAT(1X,'INPUT SIG PARMS: GAMMA,FREQ,NOISE')
C READ 220, GAMMA,FC,RNOISE
220 FORMAT(F7.1,F7.1,F7.1)
C ADDED BY MAURI
AVAR = 4.
BVAR = 4.
RATE = 20.
FC = 4000.
FCHIRP = 20.
RNOISE = 4.
PHIVAR = 4.
TCHIRP = 2.
GAMMA = 1.
ASIGMA=SQRT(AVAR)
BSIGMA=SQRT(BVAR)
PHISGM=SQRT(PHIVAR)
RSIGM=SQRT(RNOISE)
DMEAN=1200./RATE
RC=6.28319*FC
WCHIRP = 6.2831*FCHIRP
C IF (ESEP(1).NE.0.) GO TO 500
ESEP(1)=1.
ESEP(2)=3.
ESEP(3)=1.
ESEP(4)=3.
ESEP(5)=7.
ESEP(6)=14.
500 RETURN
END

95
tmpsave/kalfil.f
Wyświetl plik

@ -1,95 +0,0 @@
SUBROUTINE KALFIL (Z, IP, RN, ILX, IXS, KELEM, JNODE, ISRATE, DUR,
&ILRATE, PIN, LKHDJ)
C 2 ILRATE, PIN, LKHDJ)
C
C THIS SUBROUTINE COMPUTES THE ARRAY OF KALMAN FILTER
C RECURSIONS USED TO DETERMINE THE LIKELIHOODS.
C
C VARIABLES:
C Z - INPUT MEASUREMENT
C IP - INPUT PATH IDENTITY
C RN - INPUT NOISE POWER ESTIMATE
C ILX - INPUT SAVED KEYSTATE ON PATH IP
C IXS - INPUT KEYSTAT OF NEW NODE
C KELEM - INPUT ELEM STATE OF NEW NODE
C ISRATE - INPUT SPEED STATE OF NEW NODE
C OUR - INPU CURRENT DURATION OF ELEMENT ON IP
C ILRATE - INPUT SPEED STATE ON PATH IP
C PIN - TRANS PROB FROM PATH IP TO NODE N
C LKHDJ - OUTPUT CALCULATED LIKELIHOOD VALUE
C
C SUBROUTINES USED
C MODEL - OBTAINS THE SIGNAL-STATE-DEPENDENT LINEAR
C MODEL FOR THE KALMAN FILTER RECURSIONS
C
REAL LKHDJ
DIMENSION YKKIP(25), PKKIP(25)
DIMENSION YKKSV(750), PKKSV(750)
COMMON/BLKSV1/YKKIP, PKKIP, YKKSV, PKKSV
DATA YKKIP/25*.5/,PKKIP/25*.10/
DATA PINMIN/.00010/
C IF TRANSITION PROBABILITY IS VERY SMALL, DON'T
C BOTHER WITH LIKELIHOOD CALCULATION:
C
IF (PIN.GT.PINMIN) GO TO 100
LKHDJ = 0.
GO TO 400
C
C OBTAIN STATE-DEPENDENT MODEL PARAMETERS:
C
100 CALL MODEL(DUR, KELEM, ILRATE, ISRATE, IXS, PHI, QA, HZ)
C
C GET PREVIOUS ESTIMATES FOR PATH IP
C
YKK = YKKIP(IP)
PKK = PKKIP(IP)
C
C IMPLEMENT KALMAN FILTER FOR THIS TRANSITION
C
YPRED=PHI*YKK
PPRED=PHI*PKK*PHI+QA
PZ=HZ*PPRED+RN
PZINV=1./PZ
G=PPRED*HZ*PZINV
PEST=(1.-G*HZ)*PPRED
ZR=Z-HZ*YPRED
YKKSV(JNODE)=YPRED+G*ZR
PKKSV(JNODE)=PEST
IF(YKKSV(JNODE) .LE. .01) YKKSV(JNODE)=.01
A=0.5*PZINV*ZR**2
IF(A .LE. 1000.) GO TO 200
LKHDJ=0.
GO TO 400
200 LKHDJ=(1./SQRT(PZ))*EXP(-A)
GO TO 400
PRINT 3, Z,HZ,QA,PHI,PZ,ZR,G,PEST,YKK,YKKSV(JNODE),LKHDJ
3 FORMAT(1X,11(F6.3,1X),/)
400 RETURN
END

113
tmpsave/key.f
Wyświetl plik

@ -1,113 +0,0 @@
SUBROUTINE KEY(DUR,X)
DIMENSION ESEP (6),EDEV(6),MORSE(10,40)
DIMENSION IOUTP(500)
DIMENSION ITEXT(200)
DIMENSION ISYMBL(6)
C CHARACTER*6 ISYMBL
COMMON/BLKEND/IEND
COMMON/BLK1/TAU/PLK6/DMEAN,XDUR,ESEP,EDEV
COMMON/BLKTXT/ITEXT
DATA IK/0001000000000/
DATA LTR/20/,NELM/0/, N/0/, NLTR/1/
DATA MORSE/1,3,2,0,0,0,0,0,0,0,
&2,3,1,3,1,3,1,0,0,0,2,3,1,3,1,3,1,0,0,0,
&2,3,1,3,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,
&1,3,1,3,2,3,1,0,0,0,2,3,2,3,1,0,0,0,0,0,
&1,3,1,3,1,3,1,0,0,0,1,3,1,0,0,0,0,0,0,0,
&1,3,2,3,2,3,2,0,0,0,2,3,1,3,2,0,0,0,0,0,
&1,3,2,3,1,3,1,0,0,0,2,3,2,0,0,0,0,0,0,0,
&2,3,1,0,0,0,0,0,0,0,2,3,2,3,2,0,0,0,0,0,
&1,3,2,3,2,3,1,0,0,0,2,3,2,3,1,3,2,0,0,0,
&1,3,2,3,1,0,0,0,0,0,1,3,1,3,1,0,0,0,0,0,
&2,0,0,0,0,0,0,0,0,0,1,3,1,3,2,0,0,0,0,0,
&1,3,1,3,1,3,2,0,0,0,1,3,2,3,2,0,0,0,0,0,
&2,3,1,3,1,3,2,0,0,0,2,3,1,3,2,3,2,0,0,0,
&1,3,1,3,2,3,2,3,2,0,1,3,1,3,1,3,2,3,2,0,
&1,3,1,3,2,3,2,3,2,0,1,3,1,3,1,3,2,3,2,0,
&1,3,1,3,1,3,1,3,2,0,1,3,1,3,1,3,1,3,1,0,
&2,3,1,3,1,3,1,3,1,0,2,3,2,3,1,3,1,3,1,0,
&2,3,2,3,2,3,1,3,1,0,2,3,2,3,2,3,2,3,1,0,
&2,3,2,3,2,3,2,3,2,0,40*0/
C 2 2,3,1,3,1,3,1,0,0,0,2,3,1,3,1,3,1,0,0,0,
C 2 2,3,1,3,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,
C 2 1,3,1,3,2,3,1,0,0,0,2,3,2,3,1,0,0,0,0,0,
C 2 1,3,1,3,1,3,1,0,0,0,1,3,1,0,0,0,0,0,0,0,
C 2 1,3,2,3,2,3,2,0,0,0,2,3,1,3,2,0,0,0,0,0,
C 2 1,3,2,3,1,3,1,0,0,0,2,3,2,0,0,0,0,0,0,0,
C 2 2,3,1,0,0,0,0,0,0,0,2,3,2,3,2,0,0,0,0,0,
C 2 1,3,2,3,2,3,1,0,0,0,2,3,2,3,1,3,2,0,0,0,
C 2 1,3,2,3,1,0,0,0,0,0,1,3,1,3,1,0,0,0,0,0,
C 2 2,0,0,0,0,0,0,0,0,0,1,3,1,3,2,0,0,0,0,0,
C 2 1,3,1,3,1,3,2,0,0,0,1,3,2,3,2,0,0,0,0,0,
C 2 2,3,1,3,1,3,2,0,0,0,2,3,1,3,2,3,2,0,0,0,
C 2 1,3,1,3,2,3,2,3,2,0,1,3,1,3,1,3,2,3,2,0,
C 2 1,3,1,3,2,3,2,3,2,0,1,3,1,3,1,3,2,3,2,0,
C 2 1,3,1,3,1,3,1,3,2,0,1,3,1,3,1,3,1,3,1,0,
C 2 2,3,1,3,1,3,1,3,1,0,2,3,2,3,1,3,1,3,1,0,
C 2 2,3,2,3,2,3,1,3,1,0,2,3,2,3,2,3,2,3,1,0,
C 2 2,3,2,3,2,3,2,3,2,0,40*0/
DATA ISYMBL/1H.,1H-,1H ,1H;,1H/,1H~/
C CHARACTER*500 IOUTP
BETA=1000.*TAU*DUR
IF(BETA .LT. XDUR) GO TO 200
NELM=NELM+1
IELM=MORSE(NELM,LTR)
IF(IELM.NE. 0) GO TO 100
NELM=0
Y=RAND(IK)
IELM=4
IF(Y .GT. .9) IELM=5
IF((Y.LE. .9).AND.(Y.GT. .3)) IELM=6
Y=RAND(IK)
Y=35*(Y-0.001)+1.
IY=Y
LTR=IY+1
GO TO 100
NLTR=NLTR+1
LTR=ITEXT(NLTR)
IF(LTR .EQ. 0) IELM=4
IF(LTR .EQ. 37)IELM=5
IF(LTR .EQ. 38)IELM=6
NLTR=NLTR+1
LTR=ITEXT(NLTR)
100 N=N+1
IOUTP(N)=ISYMBL(IELM)
IF(N .LT. 300) GO TO 110
N=0
NLTR=0
TEND=1
PRINT 90
90 FORMAT(1X,'END OF RUN; INPUT DATA WAS:')
DO 10 K=1,10
K1=(K-1)*50+1
K2=K*50
PRINT 105,(IOUTP(L), L=K1,K2)
105 FORMAT(1X,50A1)
10 CONTINUE
READ 105,KAIT
110 XM=ESEP(IELM)*DMEAN
XSIGM=EDEV(IELM)*DMEAN
Y=RAND(IK)
Y=2.*(Y-.5)
XDUR=XM+Y+XSIGM
IF(XDUR .LT. 20.) XDUR=20.
X=1.
IF(IELM .GE. 3) X=0
200 RETURN
END

Some files were not shown because too many files have changed in this diff Show More