' Doppler correction for movement during sampling time of FFT
' Control-line circular flight path
DIM fmax(2000), Fav(2000), anglemax(2000)
pi = 4 * ATN(1): dtr = pi / 180: rtd = 1 / dtr: ftm = .254: CLS
PRINT " This QuickBasic code determines RPM and airspeed from flights recorded"
PRINT " on audio tape and subsequently analysed using Richard Hornes"
PRINT " Spectrogram computer program, as described elsewhere on my website"
PRINT " at www.supercoolprops.eftel.com"
PRINT
PRINT " The following settings are required:"
PRINT " 1. Line length 17.67m (F2A)"
PRINT " 2. Mike placed 40 feet back from circle edge"
PRINT " 3. Spectrogram settings:"
PRINT " a. FFT Size: 2048"
PRINT " b. Sample rate: 11058 Hz"
PRINT
PRINT "Warning: this program is a bit slow to run on older computers"
PRINT
R = 17.67 ' F2A line length .... radius of flight circle in feet
L = 40 * ftm ' distance of mike from edge of circle
' NB: L is insensitive variable
Vs = 340 ' velocity of sound in m/s
fm = 1 ' static frequency, no motion
Ns = 2048 ' FFT samples number of points
Ss = 11058 ' sample rate
St = Ns / Ss ' sample time in seconds
V = 80 ' guessed speed of aircraft m/s
' Doppler calculations by Joe Supercool
' no wind correction
' microphone must be at same height as model and back 20 metres outside
' the circle
' Need Richard Hornes Spectrogram program to get Doppler-shifted frequencies
' Enter zero to exit
INPUT " Harmonic number ? "; harmonic ' Count up off of Spectrogam plot
13 Fcoming = 843
Fgoing = 557
INPUT " enter 0 to exit...input Fcoming = ", Fcoming
IF Fcoming = 0 THEN GOTO 14
INPUT " Fgoing = ", Fgoing
PRINT "Fcoming = ", Fcoming
PRINT "Fgoing = ", Fgoing
nn = 0: speedmean = 0
FOR jjj = 1 TO 20 ' loop to converge guessed speed to actual
D = V * St ' distance travelled in sample time
angl = D / R ' arc angle travelled
IF D > 15 THEN GOTO 2 ELSE 3
2 PRINT "Distance travelled is too great"
PRINT "Reduce sample size, increase sample rate"
STOP
3 i = 0
FOR dummy = 20 * dtr TO 250 * dtr STEP 2 * dtr ' sweep thru arcs
f = 0: fprev = 0: i = i + 1
FOR theta = 0 TO angl STEP angl / 40
theta1 = theta + dummy
k = 1: IF theta1 > pi THEN k = -1
phi = pi - theta1
S = SQR(R ^ 2 + (R + L) ^ 2 - 2 * R * (R + L) * COS(phi))
dum = -((R + L) ^ 2 - R ^ 2 - S ^ 2) / (2 * R * S)
epsilon = ATN(SQR((1 / dum) ^ 2 - 1))
eta = pi / 2 - epsilon
fo = fm * Vs / (Vs - k * V * COS(eta)) ' Doppler shifted frequency
f = f + fo ' accumulate to get mean-over-angles
IF fo > fprev THEN fmax = fo
IF fo > fprev THEN anglemax = theta1
fprev = f
NEXT theta
Fav(i) = f / 41: fmax(i) = fmax: anglemax(i) = anglemax
NEXT dummy
n = i
' Find upper fmax
max = 1
FOR i = 1 TO n
IF fmax(i) > max THEN ii = i
IF fmax(i) > max THEN max = fmax(i)
NEXT i
' Find lower fmax
min = 3
FOR i = 1 TO n
IF fmax(i) < min THEN iii = i
IF fmax(i) < min THEN min = fmax(i)
NEXT i
' Find upper fmaxav .. average frequency
maxav = 1
FOR i = 1 TO n
IF Fav(i) > maxav THEN jj = i
IF Fav(i) > maxav THEN maxav = Fav(i)
NEXT i
' Find lower fminav
minav = 3
FOR i = 1 TO n
IF Fav(i) < minav THEN iii = i
IF Fav(i) < minav THEN minav = Fav(i)
NEXT i
uppercorrn = max / maxav
lowercorrn = min / minav
Fcomingx = Fcoming / harmonic * uppercorrn ' get fundamental
Fgoingx = Fgoing / harmonic * lowercorrn
f = 2 / (1 / Fcomingx + 1 / Fgoingx)
RPM = f * 60
V = (1 - f / Fcomingx) * 340
speed = V * 60 * 60 / 1000
nn = nn + 1
speedmean = speedmean + speed
NEXT jjj
RPM = INT(RPM / 100) * 100
PRINT
PRINT USING " RPM = ##### Speed kph = ###"; RPM; speedmean / nn
PRINT : GOTO 13
14 PRINT "end program"
END