Approved For Release 209Q/fi8fta :£IA-ApP9fe(M788R001 8001 60001 -6
(u)
and
the
Stanford mean solar magnetic
field
are updated
on
a weekly
or
monthly basis
as
the data
are available.
ELF data
from the
"old" system
are recorded
twice
daily;
ELF data
from the
"new"
system are being
recorded
at
the
following
intervals:
half-hourly averages
(48
times
a
day)
and
two sets
of half-daily averages
(twice
a
day)
at 00:00
and
12:00 UT,
and at 04:30 and
16:30 UT.
Whole-day averages
are also being recorded.
(U)
To
summarize
the
status
of
the satellite-downlink data-
acquisition
system:
•
The downlink-geophysical
system is
in place and
in operation
at
the Los Altos
site.
•
Partial data acquisition from the downlink
is
in place.
•
Software
for
the final data-acquisition system will be
completed and
implemented
in August.
E.
(U)
Geophysical Data/RV Performance Correlation Analysis
(U)
Data
for RV sessions are
to be analyzed
for statistical cor-
relation with respect
to
the battery
of geophysical data
sets
listed
in
Table
1,
and
those acquired by Time Research Institute via measurement
(Figure
3).
Recording of RV and geophysical data
is now
in progress.
The overall system for data acquisition and analysis
of RV performance/
geophysical data
is depicted
in Figure
4.
When enough data have been
collected
toward
the
end
of
the contract period,
analysis will be performed
The
tasks described
in earlier
sections are
in preparation
for
this
task,
and
therefore have
received
the bulk of
the effort.
Certain subtasks
in
the analysis
task, which require
a
longer
leadtime,
have,
however,
already
been completed
in preparation
for
the analysis.
(U)
The primary statistical program that will
be used
to
scan
the
data
for possible relationships
is called EPOCH ANALYSIS.
This program
reads
two
files
simultaneously.
The
first file
is an event
file,
the
second
a data
file.
The program first reads an event,
then scans
the
data
temporally backward and
forward
in
time around
the
event.
This
information
is
stored,
a
second event
is
read,
and
so
forth.
When all
Approved For Release 2dUbM6iASS4fpfeE-Bb788Rooi sooi 60001 -6
Approved For Releas^|(^^0g/^5§lj^jlg^-00788R001 8001 60001 -6
LOS ALTOS STATION
SRI STATION
UNCLASSIFIED
FIGURE
(U)
REAL-TIME GEOPHYSICAL DATA-ACQUISITION SYSTEM
(U)
the events and
surrounding data have been read,
a printout
is created
that
lists appropriate cross-correlation statistics between event
and
data
elements.
(U)
Preliminary
scans
of data
generated during an approximate
100-site
series with one remote viewer have been carried
out.
The
session quality was
graded
on
a
scale
of
to
,
and correlations
between solar magnetic
field
and
solar
sunspot number were
investigated.
Some correlation between RV performance
and
solar
sunspot number was
found,
which,
if substantiated
by
further data, would
indicate
the
possibility
that performance might improve immediately after
a peak in
the
sunspot number,
and would deteriorate just before
the
sunspot number
peaks
in
its
27-to-29-day cycle.
This result
is based
on data points
that
are
too small
in number
to be
taken seriously at
the point,
however,
and
is mentioned only
to give
an example
of
the
types
of correlations
that will
be
sought out and examined during
the course
of
study.
Approved For ReleasIfeNfiAAS SiUMiBOb-00788R001 8001 60001 -6
Approved For Release 20y/p|^:^lg-gj)p|^(g788R001 8001 60001 -6
Approved For Release 2(jbbN^QlA:SS4fitB-B3b788R001 8001 60001 -6
Approved For Release
6-00788R001 8001 60001 -6
£.2000/08/08
'
IV
SUMMARY
(U)
(U)
Approximately 80 percent
of
the project's data-acquisition
tasks have been completed,
and all of
the equipment and hardware have
been purchased and delivered.
(U)
System calibration
should
be
completed
in
the near
future,
and ELF and downlink-data acquisition will have begun
in their
final
formats
.
(U)
At
the above point,
the
focus
of effort will
turn
to analysis
of past and present geophysical and ELF data,
soon
to be
followed by
the
initiation of correlation studies
of
these data against RV performance.
Approved For Releas^lfgq^pg/^SSfPf&l^-007881^001 8001 60001 ~6
Approved For Release 20mstf&gim
788R001 8001 60001 -6
THIS PAGE INTENTIONALLY LEFT BLANK
Approved For Release 20UM(!!1AS81PV^0788R001 8001 60001 ~6
Approved For Release 2000/08/08
: CIA-RDP96-00788R001 8001 60001 -6
UNCLASSIFIED
Appendix A
FAST FOURIER TRANSFORM ROUTINE FOR ELF DATA
(U)
Approved For Release ||j|^^0^:^^^|g^-00788R001 8001 60001 -6
Approved For Release 2000/08/08
: CIA-RDP96-00788R001 8001 60001 -6
UNCLASSIFIED
i
!
THIS PAGE INTENTIONALLY LEFT BLANK
Approved For Release 20
mmstm
0788*00 ' 8001 60001
Approved For Releas^jj^tf)g/^jg^j?g^6-00788R001 8001 60001 -6
Appendix A
FAST FOURIER ROUTINE FOR ELF DATA
(U)
*
*
*
PROGRAM DESCRIPTION
*
*
*
************************************************************************
THIS
PROGRAM USES
THE
FAST FOURIER TRANSFORM
( FFT > ALGORITHM TO
CALCULATE
THE
SPECTRAL
MAGNITUDES
OF AN ARRAY
OF CLOSE-PACKED
REAL
INPUT
DATA POINTS.
THE PROGRAM
IS ORGANIZED AS
FOLLOWS:
1.
AT
THE
START,
DATA POINTS ARE
STORED
IN
THE
REAL ARRAY
X<1>,
WHERE
RUNS FROM Z TO
NX-1.
NX MUST
BE AN
INTEGRAL
POWER OF
2,
OF
THE
FORM
NX
-
2^MX
(IE,
MX
=
LOG2(NX>>.
SUCCESSIVE
ELEMENTS
OF
X
REPRESENT SUCCESSIVE
SAMPLES
OF AN
INPUT SIGNAL,
SAMPLED AT
REGULAR
INTERVALS
OF
TIME
DT.
2.
THE
REAL ARRAY
X<
)
IS TREATED
FOR THE
FFT AS
A COMPLEX ARRAY
OF
NX/2
PAIRS
OF
REAL
AND
IMAGINARY ELEMENTS.
THAT
IS,
THE
REAL
ELEMENTS
OF
THE ARRAY ARE
THE
EVEN-NUMBERED
INPUT SAMPLES,
AND
THE
IMAGINARY
ELEMENTS
OF
THE ARRAY ARE
THE ODD-NUMBERED
INPUT
SAMPLES.
AFTER
THE
FFT
IS CALCULATED,
AN ADDITIONAL
STEP
IS
USED
TO EXTRACT THE
SPECTRUM OF
THE
REAL
INPUT DATA.
TRANSFORMING
CLOSE-PACKED DATA
IN THIS WAY,
EVEN THOUGH
IT REQUIRES AN
ADDITIONAL
STEP
FOR THE
REAL TRANSFORM EXTRACTION,
IS FASTER THAN
TRANSFORMING A COMPLEX ARRAY
OF
NX
REAL- 1 MAG I NARY
PAIRS
<2*NX
ELEMENTS),
WHERE ALL
OF
THE
INITIAL
IMAGINARY VALUES ARE
ZERO.
t
3.
THE
FFT
IS CALCULATED AS
FOLLOWS:
A.
THE
INPUT PAIRS
IN
X ARE
REARRANGED
IN
BIT-REVERSED
ORDER.
B.
THE
FFT
IS CALCULATED
USING
DECIMATION
IN
TIME,
WITH
ANGLE ARGUMENTS
IN
EACH SUB-DFT APPEARING
IN
NATURAL
( IE
,
INCREASING
)
ORDER.
C.
AFTER THE
LAST
PASS THROUGH
THE
FFT ALGORITHM,
THE
ARRAY
X CONTAINS
SPECTRAL
VALUES
IN
NORMAL
ORDER,
WITH
EACH
EVEN
POINT A
REAL
VALUE AND
EACH
SUBSEQUENT ODD
POINT THE CORRESPONDING
IMAGINARY
VALUE.
4.
FOLLOWING
THE
FFT THE
REAL TRANSFORM
IS
EXTRACTED.
NOTE THAT
IF
THE
SAMPLING
INTERVAL
(THE
TIME
BETWEEN
SUCCESSIVE
SAMPLES)
IS
DT,
THEN
THE
TOTAL
SIGNAL
INTERVAL
PROCESSED
IS
NX*DT
(THE
FFT
ASSUMES THAT THE MISSING
RIGHT
END POINT
IS
THE
SAME
AS
THE
FIRST
POINT).
THE
FFT GENERATES
A SPECTRUM CONTAINING VALUES AT
INCREMENTS
OF
DF
IN
FREQUENCY,
WHERE
DF
“
1/(NX*DT).
IF
WE WERE
TO TRANSFORM AN ARRAY
OF
NX COMPLEX
DATA POINTS
(WHOSE
IMAGINARY
VALUES WERE
ZERO,
SINCE
WE
ARE
CONCERNED WITH A
REAL
SIGNAL)
WE
WOULD GENERATE
NX COMPLEX
SPECTRAL
POINTS
RANGING
IN
FREQUENCY
FROM Z TO
( NX- 1 >*DF
Hz,
HOWEVER,
THE
POINTS FROM
( NX/2+1 )*DF
Hz
TO
< NX- 1 )*DF
Hz
ARE
MERELY
THE
COMPLEX CONJUGATES
OF
THE
LOWER
POINTS,
AND CONTAIN
NO ADDITIONAL
INFORMATION.
THIS
IS
BECAUSE
THE
NYQUIST
RATE,
OR
THE
HIGHEST UNALIASED
FREQUENCY THAT CAN
BE
SAMPLED,
IS
1 /2*NX*DF
.
WHEN WE TRANSFORM
NX
REAL CLOSE-PACKED
POINTS AND
THEN
EXTRACT THE
REAL
TRANSFORM,
WE
GET ONLY
THE
LOWER
NX/2 COMPLEX
SPECTRAL
VALUES,
WHICH ARE ALL
THAT ARE
NEEDED.
Approved For ReleasUMGl8A£SlEIGD96-00788R001 8001 60001 -6
Approved For Release 20C^|^^^l^j)^|-^788R001 8001 60001 -6
5.
NEXT,
THE
SPECTRUM
IS CONVOLVED WITH A SHORT WINDOW FUNCTION.
THE
REAL AND
IMAGINARY
(EVEN AND ODD)
POINTS ARE CONVOLVED
SEPARATELY,
SINCE
THE WINDOW FUNCTION
IS A SEQUENCE
OF
REAL
NUMBERS.
WINDOWING
IS NECESSARY
IF
THE
SIGNAL
FILTERS WE ARE
EFFECTIVELY SYNTHESIZING ARE
TO
HAVE
A USEFUL
SHAPE.
WITHOUT
WINDOWING WE WOULD
FIND
THAT
EACH
SPECTRAL
FILTER WOULD
HAVE
A
NARROW PASSBAND BUT SIGNIFICANT SIDELOBE
RESPONSES:.
THAT
IS,
THE
MAGNITUDE
OF
A SPECTRAL
LINE XU) WOULD
DEPEND
NOT ONLY
ON
SIGNAL
COMPONENTS
NEAR
1*DF
IN
FREQUENCY,
BUT ALSO SIGNIFICANTLY
ON
COMPONENTS AT OTHER FREQUENCIES AS WELL.
WINDOWING
BROADENS
THE
SHAPE
OF
THE
PASSBAND
OF
EACH SPECTRAL
FILTER BUT DECREASES
THE
SIDELOBE
RESPONSES.
THE AMOUNT
OF
BROADENING AND iSUPPRESSION
OF
SIDELOBES DEPENDS
ON
THE WINDOW ORDER,
OR
THE
LENGTH
OF
THE
WINDOW FUNCTION WHICH
IS CONVOLVED WITH
THE
RAW SPECTRUM.
6.
FINALLY,
THE MAGNITUDE
OF
THE
SIGNAL AT EACH SPECTRAL
FREQUENCY
IS CALCULATED AS
THE
ROOT SUM
OF
THE
SQUARES
OF
THE
REAL AND
IMAGINARY
SPECTRAL COMPONENTS
(IE,
THE VECTOR MAGNITUDE
OF
EACH COMPLEX
SPECTRAL
POINT) AND
SCALED TO
MAKE! THE
PROCESSING
GAIN
INDEPENDENT
OF
THE
SIZE
OF
THE
BLOCK LENGTH
NiX.
WE COULD
ALSO CALCULATE
THE
PHASE
OF
THE
SPECTRAL
POINTS
BUT THIS
INFORMATION
ISN’T TOO VALUABLE
FOR OUR
USE.
NOTE
THAT THE
HIGHEST MEANINGFUL
SPECTRAL
FREQUENCY
DEPENDS
ON
THE ANTI-
ALIASING
LOW-PASS FILTER THAT
IS
USED WHEN THE
SIGjNAL
IS
SAMPLED.
THERE
IS
NOT MUCH MEANING TO SPECTRAL POINTS
ABOVE^ THE
FILTER
CUTOFF
FREQUENCY,
WHERE
SIGNALS ARE ATTENUATED
AND: FREQUENCY
ALIASING
BECOMES A PROBLEM.
THUS,
IN THIS
PROGRAM! WE
DO NOT
CALCULATE
SPECTRAL COMPONENTS ABOVE ABOUT
Hz
(THE
NOMINAL
CUTOFF
FREQUENCY
OF
THE
LOW-PASS
FILTER).
THIS PROGRAM
IS WRITTEN
IN
BASIC
FOR THE
IBM
PC COMPUTER.
HOWEVER,
IT
IS EASILY ADAPTED
TO OTHER MACHINES.
FOUR THINGS
TO WATCH OUT
FOR WHEN
TRANSFERRING
THE
PROGRAM TO ANOTHER VERSION
OF
BASIC ARE:
!
1.
ALL ARRAYS
START WITH
INDEX
0.
THAT
IS,
THE ARRAY
X(NX)
RUNS
FROM
X<0)
TO X(NX-l).
ARRAY REFERENCES WILL
NEED to
BE CHANGED
IF
THE
PROGRAM
IS
TO
BE
USED ON MACHINES WHERE
THE
FIRST ELEMENT
OF ARRAYS
HAS
INDEX
1.
2.
ALL VARIABLES WHOSE
NAMES
START WITH
THE
LETTERS
"I"
THROUGH
''N”
ARE
IMPLICITLY
INTEGERS
(DEFINT
I-N
STATEMENT):.
THESE
VARIABLES
MAY
NEED TO
BE
RENAMED.
IT
IS
IMPORTANT! THAT LOOP
COUNTERS AND ARRAY
INDEX VARIABLES
BE
INTEGERS
BECAUSE
INTEGER
ARITHMETIC
(WHERE AVAILABLE)
IS FASTER THAN FLOATING-POINT
ARITHMETIC.
3.
THIS PROGRAM
USES
THE
INTEGER DIVISION OPERATOR
THIS
MAY
BE
CHANGED TO
H /"
IN OTHER VERSIONS
OF
BASIC.
4.
THIS PROGRAM
USES
LOGICAL OPERATORS
("NOT",
-AND",
-OR”)
ON
INTEGER VARIABLES
(NOT LOGICAL
VARIABLES)
IN THE
BIT-REVERSAL
ROUTINE.
IF
THESE OPERATORS ARE
NOT AVAILABLE
THE! ALTERNATE
VERSION
OF
STATEMENTS
3000-3190 USING ONLY
INTEGER ARITHMETIC
MUST
BE
SUBSTITUTED.
Approved For Release 20g))|^(gti/fgg|(I(^0788R001 8001 60001 -6
Approved For Release 2000/08/08
: CIA-RDP96-00788R001 8001 60001 -6
UNCLASSIFIED
************************************************************************
*
CHOICE
OF
SAMPLING
TIME
AND BLOCK
SIZE
*
*
*
************************************************************************
THE
SAMPLING
TIME
DT AND
THE
BLOCK
SIZE
NX TOGETHER DETERMINE
THE
NUMBER
OF
SPECTRAL
POINTS CALCULATED AND THEIR SPACING
IN FREQUENCY
.
THE
CHOICE
OF
SAMPLING TIME ALSO AFFECTS THE
FREQUENCIES AT WHICH POWER-
LINE
HARMONIC
INTERFERENCE
THAT
IS
PASSED
BY THE ANTI-ALIASING LOW-PASS
FILTER WILL APPEAR
IN
THE
SPECTRUM-
THE
FOLLOWING
DISCUSSION ASSUMES
THAT THE
Hz
LOW-PASS
FILTER WILL
BE
USED TO ATTENUATE
HIGH
FREQUENCY
SIGNAL
COMPONENTS
BEFORE
SAMPLING.
IF
THE
INPUT
SIGNAL
THAT
IS
SAMPLED CONTAINED
ONLY
FREQUENCIES BELOW
Hz
IT WOULD
BE
SUFFICIENT TO SAMPLE
AT THE
NYQUIST
RATE
OF
60 SAMPLES
PER
SECOND-
HOWEVER,
SINCE
THE
FILTER
DOES
NOT
HAVE
INFINITE
ATTENUATION ABOVE
Hz
IT
IS
NECESSARY TO SAMPLE
AT A SOMEWHAT HIGHER
RATE,
AND
THEN DISCARD
THOSE
SPECTRAL
POINTS
REPRESENTING
SIGNALS ABOVE
Hz.
1-
ONE
CONSIDERATION
IN CHOOSING
THE
SAMPLING
RATE
IS
THE ACTUAL
RATE AVAILABLE WITH A GIVEN A/D CONVERTER-
WITH
THE
CMC
BUSSTER
D 1 6 A/D CONVERTER,
SAMPLING
TIMES CAN
BE AS
SMALL AS 0-0005
s.
HOWEVER,
IF
MORE
THAN
ONE
SIGNAL
IS
TO
BE
DIGITIZED AT ONE
TIME,
THE
MINIMUM SAMPLING
TIME
IS 0.01
s
PER CHANNEL,
OR A SAMPLING
RATE
OF
100 SAMPLES/SECOND-
THIS
SEEMS
LIKE
A REASONABLE
CHOICE.
2.
A SECOND CONSIDERATION
IS
THE
EFFECT OF
THE
SAMPLING
RATE
ON
THE
FREQUENCIES
OF ALIASED
POWER
LINE
HARMONICS.
ONLY
HARMONICS
AT
60 AND
120
Hz
ARE
LIKELY
TO
BE A PROBLEM.
WITH
DT
- 0,01
S,
THE
NYQUIST
RATE
IS
Hz,
SO
THE
SPECTRUM WILL CONTAIN POINTS AT
FREQUENCIES
FROM 0 TO
Hz-
ANY
Hz
SIGNAL
THAT
IS DIGITIZED
WILL APPEAR
IN
THE
SPECTRUM AT 40
Hz,
WHICH
POINT WILL
BE
THROWN
OUT,
SO
Hz
INTERFERENCE WON’T
BE A PROBLEM.
HOWEVER,
INTERFERENCE AT
120
Hz WILL APPEAR AT
Hz
IN
THE
SPECTRUM,
AND
THIS
FACT MUST
BE
KEPT
IN MIND WHEN ANALYZING
THE
DATA.
(HIGHER
HARMONICS APPEAR AS:
180 AT 20
Hz,
240 AT
Hz,
300 AT 0
Hz,
360 AT
Hz,
AND
SO ON.
)
GIVEN
A SAMPLING
TIME
DT
- 0.01
s,
WE
CAN CALCULATE
THE
SPECTRAL
SPACING
FOR DIFFERENT CHOICES
OF
BLOCK
SIZE
NX,
AND
WE
FIND
THE
FOLLOWING:
NX
DF
NUMBER
OF
POINTS
<
1.5625
Hz
128
0.7813
Hz
256
0.3906
Hz
512
0.1953
Hz
153
1024
0.0977
Hz
307
204 8
0.0488
Hz
614
Approved For
Approved For Release 200Q(^^^l^pP9|-00788R001 8001 60001 -6
»*
WINDOWING
*
; *
*
******* *********************** a******************************************
(AFTER
THE
SPECTRUM OF
THE
INPUT SIGNAL
HAS
BEEN
FOUND,
IT
IS CONVOLVED
(WITH A SHORT WINDOW SEQUENCE
TO
IMPROVE
THE
SHAPE
OF
THE
SYNTHESIZED
(FILTERS.
THOUGH THE
SPECTRAL VALUES ARE COMPLEX,
THE WINDOW SEQUENCE
IS
(A SEQUENCE
OF
REAL
NUMBERS,
SO THE CONVOLUTION
IS
DONE
SEPARATELY
FOR
(THE
REAL AND
IMAGINARY SEQUENCES
IN
THE
SPECTRUM.
THE CONVOLUTION
IS
OF
(THE
FORM
(
X(j)
<- W3 *CX(j-3)+X(j+3>]
+W2* [X(j-2)+X(j+2)3
»
+Wl*tX(j-l )+X< j+1 >J
+W0*X(j>
(WHERE
X(j)
IS
THE
REAL OR
IMAGINARY SUBSEQUENCE
OF
THE
SPECTRUM,
AND
THE
(CONVOLUTION
IS
ILLUSTRATED FOR A 4TH ORDER WINDOW SEQUENCE.
NOTE
THAT
(IF
THE CONVOLUTION
IS TO
BE
DONE
IN-PLACE,
A SMALL ARRAY MUST
BE
USED TO
(HOLD
THE
PREVIOUS VALUES
OF
X(j-3),
X(j-2>,
X(j-l),
and
X(j>
FOR
USE
IN
(SUBSEQUENT CONVOLUTIONS.
(THE
VALUES TO
USE
FOR DATA POINTS
BEYOND
THE
ENDS
OF
THE 'SPECTRUM ARE
(THE COMPLEX CONJUGATES
OF
THE
POINTS REFLECTED ABOUT THE
ENDS,
SO THE
(AUXILL IARY ARRAY OF
PREVIOUS VALUES MUST
BE
INITIALIZED
4(S
FOLLOWS:
I
Re( -3 )
-
Re( 3
)
Im<-3>
- -Im(3>
t
Re( -2 )
-
R®< 2
)
Im( -2
)
- -Im<2>
(
Re( - 1
)
- Red)
Im<-1)
- -Im(l>
POINTS AFTER THE
FAR END
OF
THE
SPECTRUM ARE
REFLECTED
IH
THE
SAME WAY,
BUT SINCE WE WON'T CALCULATE
SPECTRAL
POINTS ABOVE
Hz THERE
IS
NO
NEED
TO WINDOW POINTS THERE.
; THE
FOLLOWING WINDOW SEQUENCES ARE
THE MOST
IMPORTANT
FOR; OUR
USE.
sTHESE ARE ALL MINIMUM -SIDELOBE WINDOWS WHICH SUPPRESS
SIDiELOBE
RESPONSES
l AS MUCH AS
POSSIBLE
FOR A GIVEN WINDOW ORDER,
AND HAVE
BE
!EN NORMALIZED
l F OR UNITY
GAIN.
NOTE
THAT
THE
ORDER-1 WINDOW CORRESPOND^ TO
NO WINDOW
; AT ALL.
THE ORDER-2 WINDOW
IS
SIMILAR TO THE
HAMMING
FUNCTION.
ORDER
W0
. 00000
1 .00000
. 0000000
1 .0000000
W1
-0.42875
-0.5895613
-0.6727190
W2
0.0922278
0.1878524
W3
-0.0146337
HIGHEST
-13
dB
-43
dB
-72
dB
-98
dB
SIDELOBE
3-dB WIDTH
0.89
DF
1.30
DF
1.61
DF
1.86
DF
6-dB WIDTH
1.20
DF
1.81
DF
2.25
DF
2.62
DF
(FOR OUR
USE,
I
THINK
THE
3RD-ORDER WINDOW
IS
THE MOST APPROPRIATE,
AND
(THE
PROGRAM
USES THIS WINDOW FUNCTION.
SINCE W3
- 0
IN THIS CASE,
THE
(ACTUAL CONVOLUTION CONTAINS ONLY
3 TERMS.
SINCE W0
-
,
THE
LAST TERM
(IS OUST X(j),
WITHOUT A MULTIPLICATION.
THUS
THE WINDOWING
OF
EACH
(POINT
INCLUDES
5 ADDITIONS AND
2 MULTIPLICATIONS
FOR
EACH
OF
THE
REAL
(AND
IMAGINARY COMPONENTS.
IF
ONLY VERY
SHORT BLOCKS
(SAY
NX
-
64)
ARE
(TO
BE
PROCESSED,
THE
2ND-ORDER WINDOW MIGHT
BE MORE APPROPRIATE
SINCE
(THE
FILTERS ARE
A
BIT NARROWER
IN THAT CASE.
Approved For Release 20(Jjp|(01|B:^J(^l|ff)^p0788R001 8001 60001 -6
Approved For Release,2900/Q8/ sCIRIED ots UNCIAS 3 WE ee We He We He He We WW ei ie a he * 3* PROGRAM bal ihe * gH HH ee RO III ICI IC Ye ROR IOI IIa I 4 sDEFINITIONS: 19288 DEFINT I-N "VARIABLES STARTING “I"-"N" ARE INTEGERS 1918 NX = 1824 "DATA BLOCK SIZE 1828 NXM1 = NX-1 1838 DIM X(NX) "INPUT DATA ARRAY X( 1), 1 = 8, «22, NX-1 1848 OT = 9.81 "INPUT SAMPLE TIME IN SECONDS 1858 DF = 1./¢DT*NX) "SPECTRAL LINE SPACING IN Hz 1868 NMAX = 38./0F "NUMBER OF HIGHEST SPECTRAL LINE <= 38 Hz 1199 NF = NX\2 "FFT SIZE (NR OF COMPLEX POINTS) 1118 MF = LOG(NF)/LOG(2.) "NF = 2°MF 1129 NFD2 = NFN2 1138 NFM1 = NF-1 1149 DIM CACNF) "COSINE AND SINE ARRAYS, 1159 DIM SACNF) * ANGLES FROM & TO 18@*(NX-1)/NX 1162 PI2 = 6.2831853 1229 Wl = -%.5859613 ‘WINDOW COEFFICIENTS 1218 W2 = %.8922278 ge ee IH HH IR FO III III IOI ICI ICICI I I ROI I IOI ICI ICI III IIe ie it 3 s;CALCULATE SINE AND COSINE ARRAYS: 3 3THIS NEED BE DONE ONLY ONCE. NOTE THAT THE FFT ROUTINE USED SINES AND sCOSINES FROM 8 TO <18@ DEGREES IN STEPS OF 368/NF DEGREES, AND THE REAL sTRANSFORM EXTRACTION ROUTINE USES SINES AND COSINES UP TO <9@ DEGREES IN sSTEPS OF 368/NX DEGREES. THUS THE ODD NUMBERED POINTS IN THESE ARRAYS sFOR ANGLES ABOVE 9% DEGREES ARE NOT USED. IF MEMORY SPACE IS A PROBLEM, +THE ARRAYS CAN BE CONDENSED AND ONLY POINTS FROM @ TO <18@ DEGREES IN 3STEPS OF 36@/NF CALCULATED FOR USE BY THE FFT. IN THIS CASE, THE VALUES sUSED BY THE REAL TRANSFORM EXTRACTION WILL HAVE TO BE CALCULATED AS sNEEDED. IT IS ALSO POSSIBLE TO USE A TABLE CONTAINING VALUES FOR ANGLES sONLY FROM @ TO 45 DEGREES, BUT THIS MAKES THE LOOKUP PROCEDURE SOMEWHAT sMORE COMPLICATED, INVOLVING A POSSIBLE CHANGE OF SIGN OR EXCHANGE OF SIN sAND COS VALUES DEPENDING ON THE SEMIQUADRANT OF THE ANGLE. 15292 DA = PI2/NX "ANGLE INCREMENT 1512 FOR J = @ TO NFM1 "CALCULATE SINES AND COSINES 1528 A = DA*) 1539 CACJ) = COS(A) 1548 SA(J) = SINCA) 1558 NEXT J A FE ee ee ee ee ee ie ieee ie eI RI I iI ie ei Rie ir ie ie ie ie ie ae iI ae ie ie ae ake ae 3 3SAMPLE INPUT DATA INTO ARRAY X: THIS CODE MUST BE WRITTEN TO SAMPLE THE INPUT SIGNAL EVERY &.@1 SECOND sAND STUFF SUCCESSIVE VALUES INTO X. ONLY NX VALUES NEED BE TAKEN. sENTER HERE TO START PROCESSING A NEW BLOCK OF DATA. 32288 CODE VET TO BE WRITTEN 25 Approved For Releasd PPq@UB ARS S HE? 6-00788R001800160001-6
Approved For Release 2O(0j^^^l^g|^p|-gJ)788ROO1 8001 60001 -6
x ***»»****»«**«*»********** a********************************** ***********
;
? B IT
REVERSAL:
I
; THE
REAL
COMPACT SAMPLES
IN
X ARE
NOW TREATED AS A COMPLEX ARRAY
OF
NF
;<REAL, IMAGINARY)
PAIRS.
PAIRS AT BIT-REVERSED ADDRESSES ARE
NOW SWAPPED
; PRIOR
TO
THE
FFT ROUTINE.
?THE
FOLLOWING
ROUTINE
USES
LOGICAL
OPERATORS.
3000
K
-
I
’INITIAL
BIT-REVERSED ADDRESS
3010
FOR 0*1
TO NFM1
'SCAN ALL ADDRESSES EXCEPT
THE
FIRST
3020
L
-
NFD2
’INITIAL
BIT MASK
3030
IF
NOT ( K AND
L) GOTO 3070
’THIS
BIT SET?
3040
K
-
K AND
NOT
L
’
YES
- CLEAR
IT,
3050
L
-
L\2
’
SHIFT MASK
RIGHT,
3060
GOTO
3030
'
AND CONTINUE
3070
K
-
K OR
L
'BIT WAS CLEAR,
SET
IT
3080
IF
K
<-
GOTO 3190
'SKIP
IF ADDRESS
SAME
OR ALREADY
DONE
3090
•
0 + 0
'ARRAY ADDRESS
OF
LOWER
REAL COMPONENT
3100
K2
«
K+K
{ARRAY ADDRESS
OF
UPPER REAL COMPONENT
3110
T
-
X< 0 2
)
'SWAP
REALS
3120
X ( 02
)
-
X< K2
)
3130
X ( K2
>
-
T
3140
T
=
X< 0 2+ 1
>
'SWAP
IMAGINARIES
3150
X ( 0 2+ 1
>
-
X< K2+ 1
)
3180
X( K2+ 1
)
-
T
3190
NEXT
iTHE FOLLOWING VERSION USES ONLY
INTEGER ARITHMETIC
BUT
IS A BIT LESS
iEFFICIENT.
; 3000
FOR
-
TO
NFM1
SCAN ALL ADDRESSES
EXCEPT THE
FIRST
; 3010
K1
-
‘INITIALIZE
DIVISOR
; 3020
K
- 0
•INITIALIZE
BIT-REVERSED ADDRESS
; 3030
FOR
I
-
TO MF
; BIT-RE VERSE
ALL
MF
BITS
; 3040
K2
=
Kl\2
i 3050
K
-
( K-K2 >*2+Kl
} 3060
K 1
-
K2
; 3070
NEXT
I
1 3080
IF
K O
0 GOTO
3190
SKIP
IF ADDRESS
SAME
OR ALREADY
DONE
i 3090
-
0+0
ARRAY ADDRESS
OF
LOWER
REAL COMPONENT
? 3100
K2
-
K+K
{ARRAY ADDRESS
OF
UPPER REAL COMPONENT
; 3 1 10
T
-
X( 02
>
’SWAP
REALS
; 3 1 20
X< 02 >
-
X(K2>
; 3 1 30
X ( K2
>
-
T
; 3 1 40
T
-
X ( 0 2+ 1
)
•SWAP
IMAGINARIES
; 3150
X( 02+ 1
)
-
X ( K2 + 1
)
; 3 1 80
X ( K2+ 1
)
«
T
; 3 1 90
NEXT
Approved For Release 20UMdA83W6B0788R001 8001 60001 '6
Approved For Release
6-00788R001 8001 60001 -6
FFT ROUTINE:
THIS ROUTINE
USED
DECIMATION
IN
TIME.
THE
COMPLEX
DATA
INPUT
IS
IN BIT-
REVERSED
ORDER.
1.
THE
OUTER
LOOP MAKES
MF
PASSES THROUGH
THE
DATA.
THE
INITIAL
NUMBER
OF
SUB-DFT'S
IS
1,
AND
DOUBLES
ON
EACH
PASS.
THE
INITIAL
INCREMENT
IN ANGLE
BETWEEN
SUB-DFT'S
IS
180 DEGREES
,
AND
HALVES
ON
EACH
PASS.
2.
THE
MIDDLE
LOOP
PROCESSES
THE
SUB-DFT'S.
AT
EACH
SUB-DFT
ANGLES
INCREASE
IN NATURAL
ORDER AS 0
<I=»1>?
0,
(I»2>j
0,
45,
90,
135
(
I “3
)
;
0,
22.5,
45,
...
CI«4>;
AND
SO
ON.
3.
THE
INNER
LOOP CALCULATES
THE
BUTTERFLIES
FOR
A GIVEN ANGLE
IN
EACH
SUB-DFT.
EACH
DEC I MAT ION- 1 N-T I ME
BUTTERFLY
HAS
THE
FORM
Re< j
)
<-
Re(j >+[ReCk)*C+Im(k)*S]
ImCj)
<-
Im(j )+EIm(k)*C-Re(k)*S]
Re< k
)
<-
Re( j >-CRe( k )*C+Im(k>*SI
Im( k >
<-
Im( j )- C Im( k )*C-Re( k >*S]
WHERE
C AND
S ARE
THE
COSINE
AND
SINE
OF
THE ANGLE.
3200
3210
LM
-
LM2
- LM+LM
'INITIAL
NUMBER
OF
SUB-DFT'S
*
3220
IDA
-
NF
' INITIAL ANGLE
INCREMENT
-
180 DEG
3230
FOR
I
«
TO
MF
'OUTER
LOOP
3240
IA
- 0
' INITIAL
SUB-DFT ANGLE
- 0 DEGREES
3250
FOR
L
**
TO
LM
'MIDDLE
LOOP
3260
3270
C
-
CA( IA)
S
«
SAC IA)
'LOOK
UP
COS AND
SIN
FOR THIS ANGLE
3280
FOR
0 «L- 1
TO
NFM1
STEP
LM2
i INNER
LOOP
3290
OR
-
0+0
’ARRAY ADDRESS
OF
LOWER
REAL COMPONENT
3300
-
UR+1
'ARRAY ADDRESS
OF
LOWER
IMAG COMPONENT
3310
KR
» 0R+LM2
'ARRAY ADDRESS
OF
UPPER
REAL
COMPONENT
3320
KI
-
KR+1
’ARRAY ADDRESS
OF
UPPER
IMAG COMPONENT
3330
TR
-
XC KR )*C+X( KI
) *S
'BUTTERFLY
3340
3350
3360
3370
3380
T I
»
XCKI
> *C-X ( KR
)
X(KR)
-
X ( OR > -TR
XCKI
)
«
XCOI )-TI
XCOR)
»
XC OR >+TR
XCOI
)
«
XCOI )+TI
*S
3390
NEXT
'END
OF
INNER
LOOP
3400
IA
-
IA+IDA
'BUMP ANGLE
3410
NEXT
L
'END
OF
MIDDLE
LOOP
3340
3430
LM
*
LM2
LM2
« LM2+LM2
•DOUBLE
NUMBER
OF
SUB-DFT’S
3440
IDA
«
I DA\2
'HALVE ANGLE
INCREMENT
3450
NEXT
I
•BOTTOM
OF
OUTER
LOOP
Approved For Releas4JM<CX^
Approved For Release 20(^(^^6|:^l^|^p|^788R001 8001 60001 -6
EXTRACT CLOSE-PACKED
REAL
DATA TRANSFORM:
NOV VE
MUST EXTRACT THE
SPECTRUM
OF
THE
INPUT DATA AS
IF
IT HAD
BEEN
CALCULATED
FROM A SEQUENCE
OF
NX COMPLEX
INPUT VALUES ALL
OF WHOSE
IMAGINARY COMPENENTS WERE
ZERO .
THE ALGORITHM FOR THIS
t S AS
FOLLOWS:
1.
FOR POINTS j-1,2,
NF/2-1
CALCULATE
Re( j
>
I m< j
)
Re< k
)
Im( k
)
<-
Re( j >+Re( k )- CRe( j >-Re( k
) ] *S+ II m( j >+Im(k)]*C
<-
Im(j >-Im(k)-[R©(j
) -Re( k
) 3 *C - C I m{j j
) + Im( k
) 3 *S
<-
R«( j
) + Re(k>+CRe( j
) -Re< k
) ] *S- C I m<! j >+Im<k)3*C
<-
-Im(j
) + Im< k >- [Re( j >-Re< k
) 3 *C- C Im(
j
j >+Im<k)3*S
WHERE
k
-
NF
-j,
A
- j*360/NX
DEGREES,
C
- COS(A),
AND
S
-
S I N( A
)
.
WE ALSO
HAVE
THE
SPECIAL
POINTS:
Re(0>
<-
CRe( 0
> + Im( 0
) 3
I m( 0 )
<-
Re( NF/2
)
<-
R©< NF /2
) *2
Im(NF/2)
<-
-Im(NF/2)*2
AND
(THOUGH THESE VALUES ARE
NOT CALCULATED):
Re( NF
)
<-
CRe( 0 )-Im( 0 > ] *2
Im(NF)
IS
UNKNOWN
3500
X<0>
-
X( 0 )*2
.
'SPECIAL
CASE
FOR
X<0>
3510
X(
>
-
0,
3520
FOR
J
-
TO NFD2-1
•LOOP
FOR GENERAL
CASE
3530
JR
-
J+J
•ARRAY ADDRESS
OF
LOWER
REAL COMPONENT
3540
JI
-
OR+1
'ARRAY ADDRESS
OF
LOWER
I MAG COMPONENT
3550
KR
-
NX-JR
'ARRAY ADDRESS
OF
UPPER
REAL
COMPONENT
3560
KI
-
KR+1
•ARRAY ADDRESS
OF
UPPER
IMAG COMPONENT
3570
C
-
CA( J
)
‘LOOK
UP
COS AND
SIN
FOR
THIS: ANGLE
3580
S
-
SA( J
)
3590
RS
-
X( JR ) + X< KR
)
'REAL
SUM
3600
RD
-
X( JR )-X< KR
)
'REAL
DIFF
j
3610
QS
-
X( JI )+X( KI
)
•
I MAG
SUM
3620
QD
-
X( JI )-X< KI
I MAG
DIFF
3630
TR
- QS*C-RD*S
'TERM
FOR
REALS
3640
T I
-
RD*C + QS*S
•TERM FOR
IMAGS
3650
X < JR )
-
RS + TR
3660
X(JI>
- QS-TI
3670
X( KR
)
-
RS-TR
3680
X ( K
I
>
-
-QS-TI
3690
NEXT
J
BOTTOM OF
LOOP
i
3700
X ( NF
>
-
X ( NF )*2.
•SPECIAL
CASE
FOR
X( NF/2
)
j
3710
X( NF + 1
>
-
X( NF + 1 )*2
.
Approved For Release 20UN6tASSIfNfiD0788R001 8001 60001 -6
Approved For Releas
6-00788R001 8001 60001 -6
g**********************************A*************************************
jwiNDOW SPECTRAL
DATA TO
IMPROVE
FILTER
SHAPE:
• SEE
THE
DISCUSSION
OF WINDOWING ABOVE
.
THIS
ROUTINE
IS
FOR A 3RD-ORDER
; WI NDOW FUNCTION.
3800
R2
»
X ( 4 >
INITIALIZE
PREVIOUS VALUE ARRAY,
Re<2)
3810
R1
-
X( 2
)
'
Red)
3820
R0
=
X<0)
'
Re( 0
)
3830
Q2
»
-X ( 5
)
’
Im< 2
>
3840
Q 1
*
-X< 3
)
Im(l)
3850
Q0
-
-X(
)
'
I m( 0
)
-
0.
3860
FOR 00
- 0 TO NMAX*2
STEP
'WINDOW POINTS O
Hz
3870
*
00+1
'ADDRESS
OF
IMAG COMPONENT
3880
X(O0)
-
(R2+XIO0+4
> >*W2+<Rl+X(O0+2> >*W1+R0
'REAL COMPONENT
3890
X(01)
*
<Q2+X(01+4
)
) *W2+ ( Ql+X(01+2
) )*W1 +Q0
'IMAG COMPONENT
3900
R2
-
R
'SHIFT PREVIOUS ARRAY
3910
R 1
=•
R0
3920
R0
-
X ( 00+2
>
'
AND
INSERT
NEW VALUE
3930
Q2
=
Q1
3940
Q1
=
Q0
3950
Q0
=
X ( 0 1 +2
>
3960
NEXT
'END
OF
LOOP
. ************************************************************************
; CALCULATE MAGNITUDES:
JtHE
MAGNITUDE
OF
EACH
SPECTRAL
VECTOR
IS CALCULATED AS
i
MAGNITUDE!,))
»
SQRT [Re( j >*Re( j
) + Im(j)*Itn(j )]*SF
? WHERE
SF
«
1/NX
IS
A SCALING
FACTOR*
USED
FOR
REASONS:
j
1.
THE
SCALING FACTOR
SETS
THE
PROCESSING
GAIN*
WITH
SF
-
1/NX,
!
THE
GAIN FOR
A SINGLE-COMPONENT
SIGNAL
IS
1.
THAT
IS,
IF
THE
t
INPUT
SIGNAL
IS A SINE-WAVE AT FREQUENCY
n*DF WITH PEAK AMPLITUDE
;
1.0*
THE
SPECTRUM AT FREQUENCY
n*DF WILL ALSO
BE
1.0.
;
2*
THE
SCALING FACTOR MAKES
THE
PROCESSING GAIN
INDEPENDENT OF
;
THE
BLOCK
SIZE
NX*
SO THE
SAME
SIGNAL
CAN
BE
PROCESSED WITH
i
DIFFERENT BLOCK
SIZES
IN ORDER
TO GENERATE
DIFFERENT SPECTRA WITH
;
VARYING
FREQUENCY
RESOLUTION,
BUT
THE
MAGNITUDE
OF A NARROW-BAND
i
SIGNAL
COMPONENT WILL
BE
THE
SAME
IN
EACH
SPECTRUM.
; NOTE THAT MAGNITUDES ARE
ONLY CALCULATED
FOR SPECTRAL
POINTS
<=
Hz*
; MAGNITUDES ARE
STORED AT THE
BOTTOM
OF
THE ARRAY
X*
IN
THE
FIRST NMAX
; POI NTS
.
; TH I S
ROUTINE
CAN
BE
INCLUDED
IN
THE WINDOWING
ROUTINE ABOVE WITH A SMALL
; INCREASE
IN PROGRAM EFFICIENCY*
BUT
IT
IS
SHOWN
SEPARATELY
HERE
FOR
; CLARITY
.
4000
SF
-
l./NX
'SCALE
FACTOR
4010
FOR
- 0 TO
NMAX
'PROCESS
POINTS
<»
Hz
4020
- 0+0
'ADDRESS
OF
REAL COMPONENT
4030
K2
«
02+1
'ADDRESS
OF
IMAG COMPONENT
4040
X(0)
- SQRT<X<02)*X<02)+X<K2>*X<K2>>*SF
'CALCULATE
MAGNITUDE
4050
NEXT
'BOTTOM OF
LOOP
; END
OF
PROGRAM.
Approved For Releas4JM<C£^$$H^^