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


#if @HELP_NSPECTRUM

    NSPECTRUM

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

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

                     s - The input series or array.

                    N  - An optional integer, the SPECTRUM 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 SPECTRUM. 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 real series or array, the N point spectrum of the input.

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

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

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

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

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

    Remarks:
             NSPECTRUM computes N equally spaced samples of the SPECTRUM 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 SPECTRUM (where the number of samples is
             equla to length(s)), however, the computation is generally
             faster since a shorter SPECTRUM is computed.

             See SPECTRUM for a discussion of the spectrum normalization.

             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
             NFFT
             Spectrum
#endif



/* compute an N point spectrum by zero padding or time aliasing */
SERIES ITERATE nspectrum(s, n, wintype, fixamp)
{
        if (argc < 4)
        {
                if (argc < 3)
                {
                        if (argc < 2)
                        {
                                if (argc < 1) error("nspectrum - 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)) * n / length(s);
        }

        /* peform spectrum calculation */
        return(spectrum(s, n));
}