SCREEN 12
VIEW (0, 60)-(639, 340)
WINDOW (-2, -80)-(5, -200)
PRINT "input text file": INPUT file$
OPEN file$ FOR INPUT AS #1
'determine the size of the file (Larry, you won't be restricted this way, I hope)
num = (LOF(1)) / 32
IF num > 8192 THEN num = 8192
FOR jj = 1 TO 15
IF num > 2 ^ jj AND num < 2 ^ (jj + 1) THEN num = 2 ^ jj
NEXT jj
IF num > 8192 AND num < 16384 THEN num = 8192
DIM f(num), mag(num), psd(num), tf(num), plv(6)
'the restriction on size to 8192 forced me to store the variables;
'you may not have to do so.
FOR I = 1 TO num
INPUT #1, f(I), r, im, mag(I)
NEXT I
'generate the transfer function for our instrument, whose frequency
'is 0.918 Hz.  I recommend this be 'settable' to accommodate different values.
FOR j = 2 TO num
tf(j) = f(j) ^ 2 / SQR((f(j) ^ 2 - .918 ^ 2) ^ 2 + 1.56 * .918 ^ 2 * f(j) ^ 2)
NEXT j
'next is the heart of the PSD calculation.  .43429 multiplier to convert
'from natural log to base 10 log.
'The 8.5 e9 calibration constant should be 'settable' according to instrument.
'The factor -4.1 is to give PSD per 1/6 th decade.
FOR j = 2 TO num
psd(j) = 10 * .43429 * LOG(6.2832 ^ 3 * f(j) ^ 3 * mag(j) ^ 2 / tf(j) ^ 2 / 8.5E+09 ^ 2) + 10 * .43429 * LOG(f(j) / f(2)) - 4.1
NEXT j
'everything else that follows is just 'plotting stuff'.
FOR I = 3 TO num
x = .43429 * LOG(1 / f(I))
xp = .43429 * LOG(1 / f(I - 1))
LINE (xp, psd(I - 1))-(x, psd(I))
NEXT I
FOR j = -80 TO -200 STEP -20
FOR I = .43429 * LOG(1 / f(num)) TO .43429 * LOG(1 / f(2))
ii = ii + 1
IF j = -80 THEN plv(ii - 1) = 10 ^ I
CIRCLE (I, j), .005
NEXT I
NEXT j
ofs = .43429 * LOG(1 / f(2))
FOR ll = 0 TO 3
LOCATE 25, ofs ^ 2.2 + 12 * ll
PRINT USING "####.##"; plv(ll)
NEXT ll
FOR jj = 0 TO 6
ul = -80 - 20 * jj
LOCATE 4 + 3 * jj, 70: PRINT STR$(ul)
NEXT jj
LOCATE 27, 30: PRINT "Period  (s)"
LOCATE 13, 76: PRINT "dB"
SLEEP 100