r"""
        Compute the spectral density function.  The spectral density is
        the discrete time Fourier transform of the autocovariance
        function.  In particular,

        .. math::

            f(w) = \sum_k \gamma(k) \exp(-ikw)

        where gamma is the autocovariance function and the sum is over
        the set of all integers.

        Parameters
        ----------
        two_pi : Boolean, optional
            Compute the spectral density function over :math:`[0, \pi]` if
            two_pi is False and :math:`[0, 2 \pi]` otherwise.  Default value is
            True
        res : scalar or array_like(int), optional(default=1200)
            If res is a scalar then the spectral density is computed at
            `res` frequencies evenly spaced around the unit circle, but
            if res is an array then the function computes the response
            at the frequencies given by the array

        Returns
        -------
        w : array_like(float)
            The normalized frequencies at which h was computed, in
            radians/sample
        spect : array_like(float)
            The frequency response

def spectral_density(self, two_pi=True, res=1200):
        r"""
        Compute the spectral density function.  The spectral density is
        the discrete time Fourier transform of the autocovariance
        function.  In particular,

        .. math::

            f(w) = \sum_k \gamma(k) \exp(-ikw)

        where gamma is the autocovariance function and the sum is over
        the set of all integers.

        Parameters
        ----------
        two_pi : Boolean, optional
            Compute the spectral density function over :math:`[0, \pi]` if
            two_pi is False and :math:`[0, 2 \pi]` otherwise.  Default value is
            True
        res : scalar or array_like(int), optional(default=1200)
            If res is a scalar then the spectral density is computed at
            `res` frequencies evenly spaced around the unit circle, but
            if res is an array then the function computes the response
            at the frequencies given by the array

        Returns
        -------
        w : array_like(float)
            The normalized frequencies at which h was computed, in
            radians/sample
        spect : array_like(float)
            The frequency response

        """
        from scipy.signal import freqz
        w, h = freqz(self.ma_poly, self.ar_poly, worN=res, whole=two_pi)
        spect = h * conj(h) * self.sigma**2

        return w, spect