/*****************************************************************************
*                                                                            *
*   NPSD.SPL     Copyright (C) 2008 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Evaluates N samples of the Power Spectral Density           *
*                                                                            *
*   Revisions:   11 Mar 2008  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_NPSD

    NPSD

    Purpose: Calculates N samples of the Power Spectral Density

    Syntax:  NPSD(s, N)

                s - a series, the input.

                N - an optional integer, the number of samples used to
                    compute the PSD. Defaults to the series length.

    Returns:
             A frequency domain series.

    Example:

             W1: gsin(100, 1.0, 0.2)*5;setvunits("V")
             W2: npsd(w1)
             W3: npsd(w1, 200)
             W4: npsd(w1, 20);overp(w2, lred);overp(w3, lgreen);setsym(1, 1)

             W1 contains 100 samples of a 0.2 Hz sinewave with an
             amplitude of 5 volts.

             W2 evaluates the PSD.  W3 computes the PSD by zero padding
             the input with an additional 100 zeros. W4 computes the
             PSD by time aliasing the input to produce an input of 20
             samples. As shown by the overplots, the samples of W4
             perfectly intersect the samples of W2 and W3.

    Remarks:
             NPSD uses the FFT method to evaluate N uniformly spaced
             samples of the PSD. If N > length(s), the input is zero-padded.
             If N < length(s), the input is time-aliased.

             For a given N, NPSD returns N/2 + 1 output samples.

             See PSD for a discussion of the output units of the PSD.

  See Also:
            FFT
            POWSPEC
            PSD
            SPECTRUM
#endif


/* calculate an N point PSD */
ITERATE npsd(s, n)
{
        local p;

        if (argc < 2)
        {
                if (argc < 1) error("npsd - input series required");
                
                n = length(s);
        }

        if (n < length(s))
        {
                p = psd(rowsum(ravel(s, n)) / ceil(length(s) / n));
                p *= length(s) / n;
        }
        else
        {
                p = psd(s, n);
        }
        
        return(p);
}