' 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