python类irfft()的实例源码

def make_residual_func(samples, indices, **params):
    'closure for residual func'
    fft_size = 2
    while fft_size < indices[-1]:
        fft_size *= 2
    freqs = rfftfreq(fft_size, 5)
    ind_from = int(round(1/(params['t_max']*freqs[1])))
    ind_to = ind_from+params['n_harm']
    def make_series(x):
        'Calculates time series from parameterized spectrum'
        nonlocal freqs, ind_from, ind_to, params
        spectrum = zeros_like(freqs, 'complex')
        sign_x0 = 0 if x[0] == 0.5 else abs(x[0]-0.5)/(x[0]-0.5)
        spectrum[0] = rect(sign_x0*exp(params['scale'][0]*abs(x[0]-0.5)), 0)
        for i in range(ind_from, ind_to):
            spectrum[i] = rect(
                exp(params['scale'][1]*x[1]+params['slope']*log(freqs[i])),
                params['scale'][2]*x[2+i-ind_from]
            )
        return irfft(spectrum)
    def residual_func(x):
        'calculates sum of squared residuals'
        nonlocal samples, indices
        series = make_series(x)
        sum_err = 0
        for position, ind in enumerate(indices):
            sum_err += (series[ind]-samples[position])**2
        return sum_err
    return make_series, residual_func

def single_step_propagation(self):
        """
        Perform single step propagation. The final Wigner function is not normalized.
        :return: self.wignerfunction
        """
        expV = self.get_expV(self.t)

        # p x -> theta x
        self.wignerfunction = fft.rfft(self.wignerfunction, axis=0)
        self.wignerfunction *= expV

        # theta x  ->  p x
        self.wignerfunction = fft.irfft(self.wignerfunction, axis=0)

        # p x  ->  p lambda
        self.wignerfunction = fft.rfft(self.wignerfunction, axis=1)
        self.wignerfunction *= self.get_expK(self.t)

        # p lambda  ->  p x
        self.wignerfunction = fft.irfft(self.wignerfunction, axis=1)

        # p x -> theta x
        self.wignerfunction = fft.rfft(self.wignerfunction, axis=0)
        self.wignerfunction *= expV

        # theta x  ->  p x
        self.wignerfunction = fft.irfft(self.wignerfunction, axis=0)

        # increment current time
        self.t += self.dt

        return self.wignerfunction

def single_step_propagation(self):
        """
        Perform single step propagation. The final Wigner function is not normalized.
        :return: self.wignerfunction
        """
        # p x -> theta x
        self.wignerfunction = fft.rfft(self.wignerfunction, axis=0)
        self.wignerfunction *= self.expV

        # theta x  ->  p x
        self.wignerfunction = fft.irfft(self.wignerfunction, axis=0)

        # p x  ->  p lambda
        self.wignerfunction = fft.rfft(self.wignerfunction, axis=1)
        self.wignerfunction *= self.expK

        # p lambda  ->  p x
        self.wignerfunction = fft.irfft(self.wignerfunction, axis=1)

        # p x -> theta x
        self.wignerfunction = fft.rfft(self.wignerfunction, axis=0)
        self.wignerfunction *= self.expV

        # theta x  ->  p x
        self.wignerfunction = fft.irfft(self.wignerfunction, axis=0)

        return self.wignerfunction

def istft(stft_signal, size=1024, shift=256,
          window=signal.blackman, fading=True, window_length=None):
    """
    Calculated the inverse short time Fourier transform to exactly reconstruct
    the time signal.

    :param stft_signal: Single channel complex STFT signal
        with dimensions frames times size/2+1.
    :param size: Scalar FFT-size.
    :param shift: Scalar FFT-shift. Typically shift is a fraction of size.
    :param window: Window function handle.
    :param fading: Removes the additional padding, if done during STFT.
    :param window_length: Sometimes one desires to use a shorter window than
        the fft size. In that case, the window is padded with zeros.
        The default is to use the fft-size as a window size.
    :return: Single channel complex STFT signal
    :return: Single channel time signal.
    """
    assert stft_signal.shape[1] == size // 2 + 1

    if window_length is None:
        window = window(size)
    else:
        window = window(window_length)
        window = np.pad(window, (0, size - window_length), mode='constant')

    window = _biorthogonal_window_loopy(window, shift)

    # Why? Line created by Hai, Lukas does not know, why it exists.
    window *= size

    time_signal = scipy.zeros(stft_signal.shape[0] * shift + size - shift)

    for j, i in enumerate(range(0, len(time_signal) - size + shift, shift)):
        time_signal[i:i + size] += window * np.real(irfft(stft_signal[j]))

    # Compensate fade-in and fade-out
    if fading:
        time_signal = time_signal[
                      size - shift:len(time_signal) - (size - shift)]

    return time_signal

def istft(stft_signal, signal_len, size=1024, shift=256,
          window=signal.blackman, fading=True, window_length=None):
    """
    Calculated the inverse short time Fourier transform to exactly reconstruct
    the time signal.

    :param stft_signal: Single channel complex STFT signal
        with dimensions frames times size/2+1.
    :param size: Scalar FFT-size.
    :param shift: Scalar FFT-shift. Typically shift is a fraction of size.
    :param window: Window function handle.
    :param fading: Removes the additional padding, if done during STFT.
    :param window_length: Sometimes one desires to use a shorter window than
        the fft size. In that case, the window is padded with zeros.
        The default is to use the fft-size as a window size.
    :return: Single channel complex STFT signal
    :return: Single channel time signal.
    """
    assert stft_signal.shape[1] == size // 2 + 1

    if window_length is None:
        window = window(size)
    else:
        window = window(window_length)
        window = np.pad(window, (0, size - window_length), mode='constant')

    window = _biorthogonal_window_loopy(window, shift)

    # Why? Line created by Hai, Lukas does not know, why it exists.
    window *= size

    time_signal = scipy.zeros(stft_signal.shape[0] * shift + size - shift)

    for j, i in enumerate(range(0, len(time_signal) - size + shift, shift)):
        time_signal[i:i + size] += window * np.real(irfft(stft_signal[j]))
    time_signal = time_signal[0:signal_len]
    # Compensate fade-in and fade-out
    #if fading:
    #    time_signal = time_signal[
    #                  size - shift:len(time_signal) - (size - shift)]

    return time_signal

def single_step_propagation(self):
        """
        Perform single step propagation. The final Wigner functions are not normalized.
        """
        ################ p x -> theta x ################
        self.wigner_ge = fftpack.fft(self.wigner_ge, axis=0, overwrite_x=True)
        self.wigner_g = fftpack.fft(self.wigner_g, axis=0, overwrite_x=True)
        self.wigner_e = fftpack.fft(self.wigner_e, axis=0, overwrite_x=True)

        # Construct T matricies
        TgL, TgeL, TeL = self.get_T_left(self.t)
        TgR, TgeR, TeR = self.get_T_right(self.t)

        # Save previous version of the Wigner function
        Wg, Wge, We = self.wigner_g, self.wigner_ge, self.wigner_e

        # First update the complex valued off diagonal wigner function
        self.wigner_ge = (TgL*Wg + TgeL*Wge.conj())*TgeR + (TgL*Wge + TgeL*We)*TeR

        # Slice arrays to employ the symmetry (savings in speed)
        TgL, TgeL, TeL = self.theta_slice(TgL, TgeL, TeL)
        TgR, TgeR, TeR = self.theta_slice(TgR, TgeR, TeR)
        Wg, Wge, We = self.theta_slice(Wg, Wge, We)

        # Calculate the remaning real valued Wigner functions
        self.wigner_g = (TgL*Wg + TgeL*Wge.conj())*TgR + (TgL*Wge + TgeL*We)*TgeR
        self.wigner_e = (TgeL*Wg + TeL*Wge.conj())*TgeR + (TgeL*Wge + TeL*We)*TeR

        ################ Apply the phase factor ################
        self.wigner_ge *= self.expV
        self.wigner_g *= self.expV[:(1 + self.P_gridDIM//2), :]
        self.wigner_e *= self.expV[:(1 + self.P_gridDIM//2), :]

        ################ theta x -> p x ################
        self.wigner_ge = fftpack.ifft(self.wigner_ge, axis=0, overwrite_x=True)
        self.wigner_g = fft.irfft(self.wigner_g, axis=0)
        self.wigner_e = fft.irfft(self.wigner_e, axis=0)

        ################ p x  ->  p lambda ################
        self.wigner_ge = fftpack.fft(self.wigner_ge, axis=1, overwrite_x=True)
        self.wigner_g = fft.rfft(self.wigner_g, axis=1)
        self.wigner_e = fft.rfft(self.wigner_e, axis=1)

        ################ Apply the phase factor ################
        self.wigner_ge *= self.expK
        self.wigner_g *= self.expK[:, :(1 + self.X_gridDIM//2)]
        self.wigner_e *= self.expK[:, :(1 + self.X_gridDIM//2)]

        ################ p lambda  ->  p x ################
        self.wigner_ge = fftpack.ifft(self.wigner_ge, axis=1, overwrite_x=True)
        self.wigner_g = fft.irfft(self.wigner_g, axis=1)
        self.wigner_e = fft.irfft(self.wigner_e, axis=1)

        #self.normalize_wigner_matrix()

def istft(stft_signal, size=1024, shift=256,
          window=signal.blackman, fading=True, window_length=None):
    """
    Calculated the inverse short time Fourier transform to exactly reconstruct
    the time signal.

    :param stft_signal: Single channel complex STFT signal
        with dimensions frames times size/2+1.
    :param size: Scalar FFT-size.
    :param shift: Scalar FFT-shift. Typically shift is a fraction of size.
    :param window: Window function handle.
    :param fading: Removes the additional padding, if done during STFT.
    :param window_length: Sometimes one desires to use a shorter window than
        the fft size. In that case, the window is padded with zeros.
        The default is to use the fft-size as a window size.
    :return: Single channel complex STFT signal
    :return: Single channel time signal.
    """
    assert stft_signal.shape[1] == size // 2 + 1

    if window_length is None:
        window = window(size)
    else:
        window = window(window_length)
        window = np.pad(window, (0, size - window_length), mode='constant')

    window = _biorthogonal_window_loopy(window, shift)

    # Why? Line created by Hai, Lukas does not know, why it exists.
    window *= size

    time_signal = scipy.zeros(stft_signal.shape[0] * shift + size - shift)

    for j, i in enumerate(range(0, len(time_signal) - size + shift, shift)):
        time_signal[i:i + size] += window * np.real(irfft(stft_signal[j]))

    # Compensate fade-in and fade-out
    if fading:
        time_signal = time_signal[
                      size - shift:len(time_signal) - (size - shift)]

    return time_signal