    ;=========================================================================
    ; SFFT.TXT
    ; Keith Larson
    ; TMS320 DSP Applications
    ; (C) Copyright 1996,1997,1998
    ; Texas Instruments Incorporated
    ;
    ; This is unsupported freeware with no implied warranties or
    ; liabilities.  See the C3x DSK disclaimer document for details
    ;
    ; This and all other DSK applications can be downloaded from
    ; the Texas Instruments FTP site.
    ;
    ; Further reading and other information
    ; -------------------------------------
    ; SFFT.TXT Details how the SFFT works
    ; Designer Notebook page 22 'Fast Logrithms on a Floating Point Device'
    ; APPHELP1.TXT and APPHELP2.TXT included with the DSK software
    ;
    ; Texas Instruments FTP site.
    ;
    ;   ftp://ftp.ti.com/mirrors/tms320bbs              Main FTP site
    ;   ftp://ftp.ti.com/mirrors/tms320bbs/c3xfiles     TMS320C3x code examples
    ;   ftp://ftp.ti.com/mirrors/tms320bbs/c4xfiles     TMS320C4x code examples
    ;   ftp://ftp.ti.com/mirrors/tms320bbs/c3xdskfiles  DSK file subdirectory
    ;=========================================================================

    OVERVIEW
    --------
    SFFT.ASM uses a technique known as a 'Sliding FFT' (or SFFT), to efficiently
    calculate the spectrum of a signal on a sample by sample basis.  This unique
    ability makes the SFFT particularily well suited for applications where
    signal analysis, filtering, modulation, demodulation, or other forms of
    signal manipulation in the frequency domain must be performed in real time.
    The SFFT algorithm is also quite similar to the DFT.

    The SFFT is calculated on a sample by sample basis, essentialy being the
    equivalent of overlapped FFT's having an overlap of 1 sample.  In this
    case, the past frequency data is efficiently reused to calculate the
    frequency spectra of the next sample window by adding the frequency domain
    spectra of a new sample, while simultaneously subtracting the frequency
    domain spectra of the oldest sample.

    The SFFT can be expressed in very simple terms that do not require first-
    hand knowledge of the DFT or FFT.  In addition, the SFFT can be used to
    derive the DFT equation, which may be useful to DSP beginners, or by DSP
    experts looking for a different approach to solve a problem.

    SFFT THEORY: A BETTER WAY TO USE THE IMPULSE RESPONSE
    -----------------------------------------------------
    The SFFT can be based on the following simple concepts

     1) The property of superposition allows two or more signals to be linearly
        added to create a new signal.  It can also be concluded that a sampled
        time domain signal is the summation of a series of individual input
        samples or impulses of varying magnitude (fig 1).  Similarily, signals,
        or impulses, can be subtracted.

        If an input signal sample buffer (fig 1) of data is kept in memory, a
        sliding rectangular window of data samples (fig 2 and fig 4) can be
        constructed by adding the newest sample and subtracting the oldest
        sample (fig 3) from the previous original windowed signal (fig 2).

        The diagram below shows how the addition and subtraction of samples
        can 'slide' a window of data samples from those shown in figure 2 to
        those shown in figure 4.

           <- Older                        Newer ->

     Input          *                         *
     signal       * | *                     * |
     sample     * | | | *                 * | |
     buffer   * | | | | | *             * | | |
       T <--+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-->  Fig 1
            *               * | | | | *      T=0
                              * | | |
                                * | *
                                  *
     Original       *
     Windowed     * | *                     *
     Signal       | | | *                 * |
                  | | | | *             * | |
       T <--+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-->  Fig 2
                            * | | | | *
                              * | | |
                                * | *
                                  *
     new-old                                  *
     sample                                   |
     window      Subtract old                 |
                 T=N                          |
       T <--+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+--> Fig 3
                  |                          T=0
                  |                          Add new
                  *

     Next           *                         *
     Windowed       | *                     * |
     Signal         | | *                 * | |
                    | | | *             * | | |
       T <--+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+--> Fig 4
                            * | | | | *
                  | |         * | | |       | |
                  | |           * | *       | |
                    |             *           |
                    |  Window is time         |
                --->|  shifted one sample --->|

     2) The frequency domain response of an impulse, or single sample point
        where all other data points are zero, results in a flat frequency
        response with a magnitude in each frequency bin equal to the
        impulse input magnitude.

        Conversely, the impulse is the additive result of many sinusoidal
        frequency components.  The time when the impulse occurs within the
        sample window is determined by the phase angles of the individual
        component frequencies.  An impulses time of arrival is essentialy
        determined by a linear phase shift between each frequency bin.

     3) In the frequency domain, the addition of frequency samples also
        follow the rules of superposition.

        Essentialy the spectra of figure 3, the new-old sample window, is
        added to the spectra of Figure 2, the orginal windowed signal, to
        create the new spectra of Figure 4.  The difference is that complex
        data is used in the frequency domain to represent the phase information
        of the individual component frequencies.

     4) The summation of a series of simple impulse transforms, which have
        correspondingly simple frequency domain transforms, will result in
        the composite frequency domain transform of the signal.

     5) A sliding rectangular window is created by subtracting off the Nth
        oldest sample, which in the frequency domain will have gone through
        a multiple of 2*pi radian rotations.

        Note: In some applications complex time domain inputs may also useful.
              For this application, only the Real data from an ADC is used.

    MORE DETAIL
    -----------
    If an impulse sample occurs at T=0 the frequency response calculation is
    further simplified since the response contains only REAL and no IMAG
    components.  The transform of an impulse at T=0 is simply to store the
    magnitude of the impulse into each REAL bin, and zero the IMAG bin.

    If T!=0, the time shift creates a phase shift or complex vector rotation
    within each frequency bin.  The phase rotation angle being proportional to
    the time shift and the frequency of interest.

      If the time shift is one sample period, as used in the SFFT, special
    conditions can be applied.  At low frequencies, the amount of phase shift
    from sample to sample is low, or in the case of 0 hz, zero radians of
    phase.  At higher frequencies the phase rotation is greatest.  At the
    Nyquist frequency the vector rotation is PI/2 radians per sample, which
    corresponds to 2 samples per sine wave cycle.  Vector rotation for bins
    between DC and the Nyquist rate are proportional to the bin frequency.

      A Fourier transform also produces both negative and positive frequencies,
    which are mirror images of each other.  A practical aspect is that only
    positive frequencies need to be computed.  This is suitable for spectrum
    analysis and filtering.  The ranges for n and the resulting complex
    rotation vectors (twiddle factors) for each bin are therefore.

        Positive frequencies      0   <= n < N/2
        Negative frequencies     -N/2 <= n < 0

                                     -j*2*pi*n/N
        complex(R_phase,I_phase) = exp

        REAL_tw[n] = cos(n*2*pi/N);
        IMAG_tw[n] = sin(n*2*pi/N);

      The basic SFFT operation is a vector rotate of each previous bin value,
    add the newest sample, and subtract the oldest sample.  This is a very
    simple operation, but at the expense that all bins must be computed
    before the next input sample is ready.

       NewBinVal = (New      -        Old) + (OldBinval * vect_rotate)

                                                        -j*2*pi*n/N
       Bin[n]    = (Sample[0]-Sample[N-1]) + (Bin[n] * exp           )

    VISUALIZING THE SFFT
    --------------------
    The easiest way to visualize the SFFT is to consider that each new sample
    occurs at T=0, making each new sample all REAL in the frequency domain.
    Then by considering that the past summation has been time shifted by one
    sample, a vector rotation proportional to the frequency is applied.  A
    schematic representation for an SFFT bin is shown below.

      FREQUENCY BIN DIAGRAM (equivalent to an IIR filter)
      ---------------------------------------------------

       IN --+--------------*-------------------------//------->
            |  +--------+  | +-+                More bins
            +--|N delay |----|X|---*-----------------//------->
               +--------+  | +-+   |
                           |  ^    |
                     2^N   |  |    |  +---+
    	  K2 = K1 --------+    +->| - |<---+
                           |          +---+    |
                           |            |      |
                           |          +---+    |
                           +--------->| + |    |
                                      +---+    |
                                        |      |
                      -j*2*pi*n/N     +---+    |  ^
               K1 * exp          ---->| X |    |  | R,I
                                      +---+    |
                      complex vector    |      |
                      "rotation rate"   +------+
                      Twiddle etc...    |
                                        V
                                     Fbin OUT

      Where: Vector_rotation_rate[n-th Freq] = 2*PI * n / (N*Fs);
             K1 & K2 force convergence (see below)

    Fbin CONVERGENCE & STABILITY
    ----------------------------
    One aspect of the SFFT that must be considered is that there is a feedback
    loop which will affect the stability of the bin values.  This is similar to
    an IIR filter where in the Z domain, a pole sites on the unit circle.  To
    maintain stability, and keep the bin values from growing out of control,
    the magnitude of the complex vector rotation twiddles need to be set to
    slightly less than one, placing the pole inside the unit circle.   This
    causes the impulse energy magnitude in each bin to exponentialy decay
    towards zero.

    The effect of adding a stability factor is that by Nth bin rotation an
    impulse will have decayed to K1^N of its original magnitude.  To exactly
    subtract the Nth oldest sample, the Nth oldest sample is scaled by a
    second coefficient K2=K1^N.

    A side effect of the exponential decay is that the SFFT is now windowed by
    an exponentialy decaying window.  To minimize this effect it is important
    to keep K1 close to 1.000.  0.999 for example.

    SFFT WINDOWING
    --------------
    Unlike the FFT and DFT, SFFT windowing cannot be performed in the time
    domain since the input window is moving in time, and therefore so must the
    window function.  Instead the SFFT windowing operation is performed in the
    frequency domain using a technique known as convolution.

    The desirable effect of windowing is a multiplicative process in the time
    domain whereby the sharp discontinuities at the endpoints that accompany a
    rectangular data window are smoothed out.  Without a smoothing window these
    abrupt changes smear the frequency spectrum over many bins.

    In the frequency domain, the coefficients of most windowing functions are
    simple and do not require large storage arrays.  For the raised cosine
    window function, the coefficients are particularily simple (-.5,+1.0,-.5)
    and are easily imedded into the code as addition and subtraction.

    However, frequency domain or 'convolutional window filtering' is applied to
    the REAL and IMAG data seperately before the REAL/IMAG data is combined
    into a magnitude.  Generaly this is not a not a problem since the operation
    is fast and only occurs during output.  Another benefit is that other
    window functions are rapidly and easily implimented by picking different
    convolution coefficients.

                            RAISED COSINE WINDOW
                            --------------------
                Time domain                     Frequency domain
                -----------                     ----------------

         :               *               :
         :           *       *           :                   +1.0
         :        *             *        :                     ^
         :- - - -*- - - - - - - -*- - - -:                     |
         :      *                 *      :                     |
         :   *                       *   :                     |
         *-------------------------------*-->   <--+---0-0-0-+-+-+-0-0-0----+-->
          0                        N             -N/2        |   |         N/2
                                                             V   V
            W[n]=1-cos(2*PI*n/N)/2                         -.5 -.5

    USING SFFT.ASM FOR SPECTRUM ANALYSIS
    ------------------------------------
    If the 'SPECT_EN' variable is set to 1 (TRUE), this will configure the DSK
    analog output to be the computed spectrum of the analog input beginning at
    BIN_START and ending at BIN_END.  The output is then viewed using an
    oscilloscope which is triggered on a positive synch pulse.  The DAC output
    voltage is proportional to the log magnitude of each frequency bin.

    To help pass impulses with minimal magnitude errors, each DAC output sample
    can be repeated up to DAC_RPT times.  Also, the AIC TA register value can
    be programed to have a very high passband.  This will increase the DAC
    output distortion if used for audio, but is acceptable for visual purposes.

    Also, the BIN_START and BIN_END values do not need to begin at zero or end
    at SFFTSIZE/2.  This can be used to show that the frequency bins repeat
    in the frequency domain as predicted by the discrete fourier transform.
    The only restrictions are how much memory and CPU processing power is
    available.

    USING SFFT.ASM FOR HILBERT TRANSFORMS & ARBITRARY PHASE ANGLES FILTERS
    ----------------------------------------------------------------------
    If 'SPECT_EN' is set to 0, the output will be configured to be the summation
    of the reconstructed REAL and IMAG components.

      An arbitrary output phase angle is implimented by performing a complex
    multiplication of the REAL and IMAG components by a complex vector
    determined by the parameter 'ANGLE'.  If ANGLE=90 degrees, the Hilbert
    transform is reconstructed from the pass band SFFT bins covering BIN_START
    to BIN_END.  If ANGLE=0.0, no phase shift occurs.

      The 0 degree and matched 90 degree phase shift Hilbert Transform are
    particularily useful in telecommunications applications where the quadrature
    outputs can be used to shift the spectrum of a signal, or for various radio
    and modem modulation schemes.

    RAISED COSINE WINDOWED FILTERS
    ------------------------------
    By applying the raised cosine window to the summation of bin values, the
    REAL or IMAG filter response ripple is improved.

    The method implemented uses a series of coefficients that are applied to
    each frequency bin and then added much like an FIR filter except in the
    frequency domain.

    The coefficient values result from a) the convolution of the response of a
    raised cosine function with the signal response and b) the multiplication
    of a rectangular bandpass filter, also applied in the frequency domain.
    A group delay, or time shift, will also be seen which is equal to N/2
    plus the time it takes a signal to make it through the ADC/DAC conversion
    process.

      In the examples shown below, it can be seen that for a given passband
    bandwidth, the number of bins required is actualy 'WIDTH+2' and that the
    sign of the coefficients alternates +,-,+,-...  The endpoints, which are
    also scaled by 50%, are the result of the window coefficients and essentialy
    define the edge characteristics of the filter.

              1.0
               ^
               |
               |
               |              Length = 1 bin, raised cosine window function
           <-+-+-+->
             |   |
             V   V
           -.5  -.5

              1.0
               ^
               |  0.5
               |   ^
               |   |          Length = 2 bins
           <-+-+-+-+->
             |   |
             V   |
           -.5   |
                 V
               -1.0

             1.0  1.0
               ^   ^
               |   |
               |   |
               |   |          Length = 3 bins
           <-+-+-+-+-+->
             |   |   |
             V   |   V
           -.5   |  -.5
                 V
               -1.0

             1.0  1.0
               ^   ^
               |   |  0.5
               |   |   ^
               |   |   |      Length = 4 bins
           <-+-+-+-+-+-+->
             |   |   |
             V   |   |
           -.5   |   |
                 V   V
               -1.0  -1.0

    NON WINDOWED SFFT
    -----------------
    A special case that may be of importance is when the SFFT is used to
    compute the all pass 0' and 90' Hilbert transforms of a non-windowed
    but synchronized signal.  Note that frequency bin spreading will occur
    if the signal is not harmonicaly related to the sample window.

    For REAL summations, the input is exactly reconstructed by scaling the
    0-th or DC bin by 50%.  This scaling compensates for a 2:1 rise in signal
    level since all bin data energy, except for the 0th bin, is split equaly
    between the positive and negative frequencies.

    At the 0th bin, there is no IMAG information since no phase shift is
    applied to that bin.  A DC component for an IMAG reconstruction therefore
    does not exist.

                            N/2 SFFT R/I bins
                   +----+----+----+----+----+----+
                   |    |    |    |    |    |N   |
       IN ---->----|  0 |  1 |  2 | .. | .. |- -1|
                   |    |    |    |    |    |2   |
                   +----+----+----+----+----+----+
                    | |  | |  | |  | |  | |  | |
                   +-+| +-+| +-+| +-+| +-+| +-+|
                   |+|--|+|--|+|--|+|--|+|--|+|----->>---- REALSUM
                   +-+| +-+| +-+| +-+| +-+| +-+|
                      |    |    |    |    |    |
                     +-+  +-+  +-+  +-+  +-+  +-+
                     |+|--|+|--|+|--|+|--|+|--|+|--->>---- IMAGSUM
                     +-+  +-+  +-+  +-+  +-+  +-+

                   :   :
                -->:   :<-- 0-th bin
                   :   :

    PERFORMANCE
    -----------
    Since the SFFT needs only to compute the 'bins of interest' within the span
    of one time sample, narrow band analysis or filtering is very efficient,
    even when the effective FFT size is very large.  However, large numbers of
    bins and or high sampling rates may become impractical for a single
    processor.  In which case a traditional block style FFT or filter may be
    more practical.

    For example, in a filter application, only a few frequency bins may be
    required, with the unused bins effectively being zero since they are not
    needed for reconstuction.  The maximum sampling rate or the number of bins
    that can be calculated, is as shown.

      Ts(min) = (SFFT_cycles_per_bin * bins + loop_overhead) * nS/cycle

      Ts(min) = (7 * N/2 + 52) 40 nS

      NOTE: The loop overhead value is the time consumed by interrupt routines,
          data formating, input and output.  SFFT.ASM is not highly optimized
          since it is for educational purposes.

          Immediate optimization gains can be achieved by inlining the three
          major functions 'Input', 'SFFT', and 'Output', removing 3 calls and
          3 returns, or 24 cycles, from the loop overhead.

    LOOP UNROLLING FOR HIGH SPEED FILTERING
    ---------------------------------------
    The inner loop of the SFFT consumes 5 computational cycles but executes in
    6 cycles.  The conflict occurs from a data bus bandwidth limitation and
    results from the STF||STF operation immediately preceding a double load of
    data for the MPYF3 instruction.

    This 'null' cycle can be effectively filled by moving the filter summations
    within the loop.  This is essentially 'free' since the summation can be
    done entirely within registers, requiring no data path access.

    The +1,-1 convolutional filter coefficients for raised cosine windowing
    can be hard coded within the loop by performing subtractions that invert
    the sum each time it goes through the loop.  This avoids fetching
    coefficients from the data bus.

    Overall, the forward and reverse SFFT are computed at 6-7 cycles per
    bin depending on if both REAL and IMAG outputs are required.  The general
    case 'educational' example SFFT.ASM is slightly slower while SFFT2.ASM
    which is written for filtering.

    FITTING THE CODE AND DATA INTO MEMORY
    -------------------------------------
    If the effective desired SFFT/FFT size is 512 points, then only 256
    positive frequencies need to be computed.  With R/I twiddle and R/I SFFT
    data associated with each bin, 1024 words of memory will be required.  In
    addition, 512 words of input buffer data will be needed.

    To maximize speed, the inner loop of the SFFT uses dual access on-chip
    memory to access data at the rate of two data moves per CPU cycle.  To
    avoid program fetch conflicts, the SFFT code is loaded into the second
    on-chip SRAM block, which also holds the data buffer.

    If off chip memory is available, excellent performance can be achieved by
    placing as much SFFT bin data on-chip as possible.  The input window sample
    buffer and code can be external since the main code loop will easily fit
    inside the cache and the sample buffer is only accessed twice per SFFT
    cycle.

     Note: The SFFT only needs to calculate the difference of the input
           of the most recent and the oldest data sample one time.  This
           value is reused for all bin calculations and is kept in a register

    If circular or bit-reversed data storage were used, this would force the
    data and twiddle buffers to be on 2^N word boundaries.  In addition, the
    circular addressing registers would be consumed.  Since the overhead of
    checking and reloading of the buffer pointers is minimal, and allows non
    2^N sizes, explicit pointer testing is used in SFFT.ASM.

    USING THIS CODE WITH 'C'
    ------------------------
    If you want to use the functions in this code with a high level language
    such as 'C', you will need to perform context save and restore operations
    at the beginning and end of each function.

    TLC32040 ADC & DAC CONSIDERATIONS
    ---------------------------------
    The application file SFFT.ASM is written to use a TLC32040 Analog Interface
    Chip (AIC) connected as used in a TMS320C31 DSP Starter Kit or DSK
    (TMDS3200031).  Further documentation for the DSK is either available in
    the DSK kit, or by download from the Texas Instruments FTP site at

    ftp://ftp.ti.com/mirrors/tms320bbs              Main TMS320 FTP mirror site
    ftp://ftp.ti.com/mirrors/tms320bbs/c3xdskfiles  C3x DSK files subdirectory

    SFFT SUMMARY
    ------------

     - A time signal is comprised of a series of samples.

     - Each sample is an impulse.

     - The time signal is a time summation of a series of impulses.

     - The frequency spectra of a single impulse at T=0 is trivial to calculate.
       Essentialy being only a REAL component in each frequency bin whose
       magnitude is that of the impulse

     - The frequency spectra of a signal is the summation of the individual
       impulse responses.

     - A shift in time is essentialy a shift in phase or 'phase rotate' in
       the frequency domain

     - It is easiest to consider each new impulse as occuring at T=0, and
       performing the time shift on the past summation of samples as a whole.

     - At each bin the amount of 'phase rotation' or 'twiddle factor'  that
       is applied to each bin is proportional to the frequency of the bin.  At
       DC (n=0) the phase shift is zero and at Fnyq (n=N/2), the phase shift is
       pi radians.

     - After phase rotating each bin, simply add the new sample/impulse
       value. (Dont forget to start with each bin magnitude as zero)

     - At this point the fourier transform is a forever expanding series in
       both the time and frequency domains.

     - The Nth oldest sample is rotated n multiples of 2*pi radians, making the
       Nth oldest sample completely REAL with no IMAG component.

       At N samples of age, phase rotation = N*(n*2*pi/N)=n*2*pi

     - A 'sliding' rectangular window is created by subtracting the T=Nth
       oldest sample while adding the newest T=0 sample.  Essentialy at T=N
       each frequency bin has rotated N times and is back to 0 radians of
       phase and can be properly subtracted.

    CONCLUSION
    ----------
    SFFT.ASM is written for the DSP beginner, but is packed with features that
    also make it useful to the experienced DSP veteran.  SFFT.ASM implements a
    continuous time Fourier transform which can be used to construct filters,
    analyze spectra, or used as a general-purpose DSP teaching platform.