/*****************************************************************************
*                                                                            *
*   NFFT.SPL      Copyright (C) 2009 DSP Development Corporation             *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Performs an N point FFT via zero padding or time aliasing  *
*                                                                            *
*   Revisions:    21 Apr 2009  RRR  Creation - from NFFT.MAC                 *
*                                                                            *
*****************************************************************************/


#if @HELP_NFFT

    NFFT

    Purpose: Calculates an N point FFT by zero padding or time aliasing.

    Syntax:  NFFT(s, N, wintype, fixamp)

                     s - The input series or array.

                    N  - An optional integer, the FFT length. If N > length(s),
                         the input series is zero padded to length N.
                         If N < length(s), the input series is time
                         aliased to produce N samples of the FFT. Defaults
                         to length(s).

               wintype - An optional integer, the windowing function.
                            0: Hamming
                            1: Hanning
                            2: Rectangular (none, default)
                            3: Kasier
                            4: Flattop
                            5: Blackman

                         Defaults to no window (rectangular) if not specified.


               fixamp - An optional integer, the windowing amplitude
                        correction method.
                            0: do not correct amplitude (default)
                            1: correct amplitude
                            2: correct RMS
                            3: correct mean squared


    Returns: A series or array

    Example:
             W1: 1..12
             W2: nfft(w1, 20)
             W3: nfft(w1, 4)
             W4: decimate(fft(w1), 3)

             W2 contains 20 samples of the FFT of W1 with 8 zeros
             appended for an input length of 20.

             W3 contains 4 samples of the FFT of W1. This is
             numerically equivalent to decimating the full FFT by 3
             as shown in W4, but because a 4 point FFT is calculated,
             the computation is performed more quickly.

    Example:
             W1: 1..12
             W2: nfft(w1, 20, 3)
             W3: nfft(w1, 4, 3)
             W4: decimate(fft(kaiser(w1)), 3)

             Same as the first example except a Kaiser window is employed.

    Remarks:
             NFFT computes N equally spaced samples of the FFT by zero
             padding if N is greater than the series length or by time
             aliasing the input series if N is less than the series length.

             For N < length(s), the result is numerically equivalent to
             decimating the full FFT (where the number of samples is
             equal to length(s)), however, the computation is generally
             faster since a shorter FFT is computed.

             If a windowing function is specified, the window is applied
             to the entire input series.

             If fixamp == 1, the window correction factor is the mean of
             the spectral window. This assures that the spectrum of a
             sinusoid of amplitude 1.0 has a peak of 1.0.

             If fixamp == 2, the correction is applied as follows:

             w = winfun(s) * rms(s) / rms(winfun(s))

             where winfun is Blackman, Flattop, Hamming, Hanning or Kaiser.
             This assures that:

             sqrt(area(psd(w))) == rms(s)    approximately


             If fixamp == 3, the correction is applied as follows:

             w = winfun(s) / sqrt(win * win / length(win))

             where win is the windowing function.


    See Also:
             Blackman
             FFT
             Flattop
             Hamming
             Hanning
             Kaiser
             Nspectrum
             Spectrum
#endif


/* compute an N point FFT by zero padding or time aliasing */
SERIES ITERATE nfft(s, n, wintype, fixamp)
{
        if (argc < 4)
        {
                if (argc < 3)
                {
                        if (argc < 2)
                        {
                                if (argc < 1) error("nfft - input series required");
                                
                                n = length(s);
                        }
                        
                        wintype = -1;
                }
                
                fixamp = 0;
        }

        if (wintype != -1)
        {
                /* apply windowing function */
                s = winfunc(wintype, s, fixamp);
        }

        if (n < length(s))
        {
                /* time aliased FFT */
                s = rowsum(ravel(s, n));
        }
        
        return(fft(s, n));
}