Python scipy 模块，fftpack() 实例源码

def bandPass(df , fmin = None, fmax = None , unit = "Hz" ) :
   """
      Return filtered signal
   """
   from scipy.fftpack import rfft, irfft, rfftfreq  #Warning convention of scipy.fftpack != numpy.fft   !!!
   if type(df) == pd.Series :
      df = pd.DataFrame(df)
      ise = True
   else :
      ise = False
   filtered = pd.DataFrame( index = df.index  )
   W = rfftfreq( df.index.size, d = dx( df ) )
   for col in df.columns :
      tmp = rfft(df[col])
      if fmin is not None :
         tmp[ ( W < fmin) ] = 0
      if fmax is not None :
         tmp[ ( W > fmax) ] = 0
      filtered[col] = irfft(tmp)
   if ise :
      return filtered.iloc[:,0]
   else :
      return filtered

def arfft(x):
    if use_numpy_arfft:
        y = numpy.fft.rfft(x)
        y = abs(y)
        return y
    else:
        assert (len(x) % 2) == 0

        y = scipy.fftpack.rfft(x)

        if type(y[0]) == numpy.float32:
            cty = numpy.complex64
        elif type(y[0]) == numpy.float64:
            cty = numpy.complex128
        else:
            assert False

        # y = [ Re0, Re1, Im1, Re2, Im2, ..., ReN ]
        #y1 = numpy.sqrt(numpy.add(numpy.square(y[1:-1:2]),
        #                          numpy.square(y[2:-1:2])))
        y0 = abs(y[0])
        yn = abs(y[-1])
        y = abs(y[1:-1].view(cty))
        y = numpy.concatenate(([y0], y, [yn]))
        return y

def rfft(x):
    if use_numpy_rfft:
        y = numpy.fft.rfft(x)
        return y
    else:
        assert (len(x) % 2) == 0

        y = scipy.fftpack.rfft(x)

        if type(y[0]) == numpy.float32:
            cty = numpy.complex64
        elif type(y[0]) == numpy.float64:
            cty = numpy.complex128
        else:
            assert False

        # y = [ Re0, Re1, Im1, Re2, Im2, ..., ReN ]
        y1 = y[1:-1].view(cty)
        y = numpy.concatenate((y[0:1], y1, y[-1:]))
        return y

# apply automatic gain control.
# causes each winlen window of samples
# to have average absolute value of 1.0.
# winlen is in units of samples.

def DST(x):
    """
    Converts Scipy's DST output to Matlab's DST (scaling).
    """
    X = scipy.fftpack.dst(x,type=1,axis=0)
    return X/2.0

def IDST(X):
    """
    Inverse DST. Python -> Matlab
    """
    n = X.shape[0]
    x = np.real(scipy.fftpack.idst(X,type=1,axis=0))
    return x/(n+1.0)

def make_wave(self):
        """Transforms to the time domain.

        returns: Wave
        """
        ys = scipy.fftpack.dct(self.amps, type=3) / 2
        return Wave(ys, self.framerate)

def make_dct(self):
        amps = scipy.fftpack.dct(self.ys, type=2)
        return Dct(amps, self.framerate)

def make_wave(self):
        """Transforms to the time domain.

        returns: Wave
        """
        N = len(self.hs)
        ys = scipy.fftpack.idct(self.hs, type=2) / 2 / N
        #NOTE: whatever the start time was, we lose it when
        # we transform back
        #ts = self.start + np.arange(len(ys)) / self.framerate
        return Wave(ys, framerate=self.framerate)

def make_dct(self):
        """Computes the DCT of this wave.
        """
        N = len(self.ys)
        hs = scipy.fftpack.dct(self.ys, type=2)
        fs = (0.5 + np.arange(N)) / 2
        return Dct(hs, fs, self.framerate)