Python scipy.signal.decimate() Examples

def harvest_get_downsampled_signal(x, fs, target_fs):
    decimation_ratio = np.round(fs / target_fs)
    offset = np.ceil(140. / decimation_ratio) * decimation_ratio
    start_pad = x[0] * np.ones(int(offset), dtype=np.float32)
    end_pad = x[-1] * np.ones(int(offset), dtype=np.float32)
    x = np.concatenate((start_pad, x, end_pad), axis=0)

    if fs < target_fs:
        raise ValueError("CASE NOT HANDLED IN harvest_get_downsampled_signal")
    else:
        try:
            y0 = sg.decimate(x, int(decimation_ratio), 3, zero_phase=True)
        except:
            y0 = sg.decimate(x, int(decimation_ratio), 3)
        actual_fs = fs / decimation_ratio
        y = y0[int(offset / decimation_ratio):-int(offset / decimation_ratio)]
    y = y - np.mean(y)
    return y, actual_fs 

def harvest_get_downsampled_signal(x, fs, target_fs):
    decimation_ratio = np.round(fs / target_fs)
    offset = np.ceil(140. / decimation_ratio) * decimation_ratio
    start_pad = x[0] * np.ones(int(offset), dtype=np.float32)
    end_pad = x[-1] * np.ones(int(offset), dtype=np.float32)
    x = np.concatenate((start_pad, x, end_pad), axis=0)

    if fs < target_fs:
        raise ValueError("CASE NOT HANDLED IN harvest_get_downsampled_signal")
    else:
        try:
            y0 = sg.decimate(x, int(decimation_ratio), 3, zero_phase=True)
        except:
            y0 = sg.decimate(x, int(decimation_ratio), 3)
        actual_fs = fs / decimation_ratio
        y = y0[int(offset / decimation_ratio):-int(offset / decimation_ratio)]
    y = y - np.mean(y)
    return y, actual_fs 

def harvest_get_downsampled_signal(x, fs, target_fs):
    decimation_ratio = np.round(fs / target_fs)
    offset = np.ceil(140. / decimation_ratio) * decimation_ratio
    start_pad = x[0] * np.ones(int(offset), dtype=np.float32)
    end_pad = x[-1] * np.ones(int(offset), dtype=np.float32)
    x = np.concatenate((start_pad, x, end_pad), axis=0)

    if fs < target_fs:
        raise ValueError("CASE NOT HANDLED IN harvest_get_downsampled_signal")
    else:
        try:
            y0 = sg.decimate(x, int(decimation_ratio), 3, zero_phase=True)
        except:
            y0 = sg.decimate(x, int(decimation_ratio), 3)
        actual_fs = fs / decimation_ratio
        y = y0[int(offset / decimation_ratio):-int(offset / decimation_ratio)]
    y = y - np.mean(y)
    return y, actual_fs 

def __call__(self, pkg):
        pkg = format_package(pkg)
        wav = pkg['chunk']
        wav = wav.data.numpy()
        factor = random.choice(self.factors)
        x_lr = decimate(wav, factor).copy()
        x_lr = torch.FloatTensor(x_lr)
        x_ = F.interpolate(x_lr.view(1, 1, -1),
                           scale_factor=factor,
                           align_corners=True,
                           mode='linear').view(-1)
        if self.report:
            if 'report' not in pkg:
                pkg['report'] = {}
            pkg['report']['resample_factor'] = factor
        pkg['chunk'] = x_
        return pkg 

def resample(recording, resample_rate):
    '''
    Resamples the recording extractor traces. If the resampling rate is multiple of the sampling rate, the faster
    scipy decimate function is used.

    Parameters
    ----------
    recording: RecordingExtractor
        The recording extractor to be resampled
    resample_rate: int or float
        The resampling frequency

    Returns
    -------
    resampled_recording: ResampleRecording
        The resample recording extractor

    '''
    return ResampleRecording(
        recording=recording,
        resample_rate=resample_rate
    ) 

def get_traces(self, channel_ids=None, start_frame=None, end_frame=None):
        start_frame_not_sampled = int(start_frame / self.get_sampling_frequency() *
                                      self._recording.get_sampling_frequency())
        start_frame_sampled = start_frame
        end_frame_not_sampled = int(end_frame / self.get_sampling_frequency() *
                                    self._recording.get_sampling_frequency())
        end_frame_sampled = end_frame
        traces = self._recording.get_traces(start_frame=start_frame_not_sampled,
                                            end_frame=end_frame_not_sampled,
                                            channel_ids=channel_ids)
        if np.mod(self._recording.get_sampling_frequency(), self._resample_rate) == 0:
            traces_resampled = signal.decimate(traces,
                                               q=int(self._recording.get_sampling_frequency() / self._resample_rate),
                                               axis=1)
        else:
            traces_resampled = signal.resample(traces, int(end_frame_sampled - start_frame_sampled), axis=1)
        return traces_resampled.astype(self._dtype) 

def decimate(x, q):
    '''decimate x by a factor of q with some settings that I like'''
    return signal.decimate(x, q, 20*q, ftype='fir', axis=0, zero_phase=False) 

def wavwrite(srcfile, fs, training):
    try:
        mat = io.loadmat(srcfile)
    except ValueError:
        print('Could not load %s' % srcfile)
        return

    dat = mat['dataStruct'][0, 0][0]
    if ds_factor != 1:
        dat = signal.decimate(dat, ds_factor, axis=0, zero_phase=True)

    mn = dat.min()
    mx = dat.max()
    mx = float(max(abs(mx), abs(mn)))
    if training and mx == 0:
        print('skipping %s' % srcfile)
        return
    if mx != 0:
        dat *= 0x7FFF / mx
    dat = np.int16(dat)

    winsize = win_dur * 60 * fs
    stride = 60 * fs
    for elec in range(16):
        aud = dat[:, elec]
        for win in range(nwin):
            dstfile = srcfile.replace('mat', str(win) + '.' + str(elec) + '.wav')
            beg = win * stride
            end = beg + winsize
            clip = aud[beg:end]
            audiolab.wavwrite(clip, dstfile, fs=fs, enc='pcm16') 

def upsample_wav(wav, args, model):
  
  # load signal
  x_hr, fs = librosa.load(wav, sr=args.sr)

  x_lr_t = decimate(x_hr, args.r)
  
  # pad to mutliple of patch size to ensure model runs over entire sample
  x_hr = np.pad(x_hr, (0, args.patch_size - (x_hr.shape[0] % args.patch_size)), 'constant', constant_values=(0,0))
  
  # downscale signal
  x_lr = decimate(x_hr, args.r)

  # upscale the low-res version
  P = model.predict(x_lr.reshape((1,len(x_lr),1)))
  x_pr = P.flatten()

  # crop so that it works with scaling ratio
  x_hr = x_hr[:len(x_pr)]
  x_lr_t = x_lr_t[:len(x_pr)] 

  # save the file
  outname = wav + '.' + args.out_label
  librosa.output.write_wav(outname + '.lr.wav', x_lr_t, fs / args.r) 
  librosa.output.write_wav(outname + '.hr.wav', x_hr, fs)  
  librosa.output.write_wav(outname + '.pr.wav', x_pr, fs)  

  # save the spectrum
  S = get_spectrum(x_pr, n_fft=2048)
  save_spectrum(S, outfile=outname + '.pr.png')
  S = get_spectrum(x_hr, n_fft=2048)
  save_spectrum(S, outfile=outname + '.hr.png')
  S = get_spectrum(x_lr, n_fft=2048/args.r)
  save_spectrum(S, outfile=outname + '.lr.png')

# ---------------------------------------------------------------------------- 

def channel_preprocessing(sig, dec, fc, Fs):
    '''deinterleave IQ samples, tune to channel frequency and decimate'''
    IQ = deinterleave_IQ(sig)
    IQ_tuned = frequency_shift(IQ, fc, Fs)
    IQd = decimate(IQ_tuned, dec)
    return IQd 

def find_channel_offset(s1, s2, nd, nl):
    '''use cross-correlation to find channel offset in samples'''
    B1 = signal.decimate(s1, nd)
    B2 = np.pad(signal.decimate(s2, nd), (nl, nl), 'constant')
    xc = np.abs(signal.correlate(B1, B2, mode='valid'))
    return (np.argmax(xc) - nl)*nd 

def test_shape(self):
        # Regression test for ticket #1480.
        z = np.zeros((30, 30))
        d0 = signal.decimate(z, 2, axis=0, zero_phase=False)
        assert_equal(d0.shape, (15, 30))
        d1 = signal.decimate(z, 2, axis=1, zero_phase=False)
        assert_equal(d1.shape, (30, 15)) 

def test_basic_FIR(self):
        x = np.arange(12)
        y = signal.decimate(x, 2, n=1, ftype='fir', zero_phase=False).round()
        assert_array_equal(y, x[::2]) 

def test_basic_IIR(self):
        x = np.arange(12)
        y = signal.decimate(x, 2, n=1, ftype='iir', zero_phase=False).round()
        assert_array_equal(y, x[::2]) 

def test_bad_args(self):
        x = np.arange(12)
        assert_raises(TypeError, signal.decimate, x, q=0.5, n=1)
        assert_raises(TypeError, signal.decimate, x, q=2, n=0.5) 

def test_shape(self):
        # Regression test for ticket #1480.
        z = np.zeros((10, 10))
        d0 = signal.decimate(z, 2, axis=0)
        assert_equal(d0.shape, (5, 10))
        d1 = signal.decimate(z, 2, axis=1)
        assert_equal(d1.shape, (10, 5)) 

def test_basic(self):
        x = np.arange(6)
        assert_array_equal(signal.decimate(x, 2, n=1).round(), x[::2]) 

def _test_phaseshift(self, method, zero_phase):
        rate = 120
        rates_to = [15, 20, 30, 40]  # q = 8, 6, 4, 3

        t_tot = int(100)  # Need to let antialiasing filters settle
        t = np.arange(rate*t_tot+1) / float(rate)

        # Sinusoids at 0.8*nyquist, windowed to avoid edge artifacts
        freqs = np.array(rates_to) * 0.8 / 2
        d = (np.exp(1j * 2 * np.pi * freqs[:, np.newaxis] * t)
             * signal.windows.tukey(t.size, 0.1))

        for rate_to in rates_to:
            q = rate // rate_to
            t_to = np.arange(rate_to*t_tot+1) / float(rate_to)
            d_tos = (np.exp(1j * 2 * np.pi * freqs[:, np.newaxis] * t_to)
                     * signal.windows.tukey(t_to.size, 0.1))

            # Set up downsampling filters, match v0.17 defaults
            if method == 'fir':
                n = 30
                system = signal.dlti(signal.firwin(n + 1, 1. / q,
                                                   window='hamming'), 1.)
            elif method == 'iir':
                n = 8
                wc = 0.8*np.pi/q
                system = signal.dlti(*signal.cheby1(n, 0.05, wc/np.pi))

            # Calculate expected phase response, as unit complex vector
            if zero_phase is False:
                _, h_resps = signal.freqz(system.num, system.den,
                                          freqs/rate*2*np.pi)
                h_resps /= np.abs(h_resps)
            else:
                h_resps = np.ones_like(freqs)

            y_resamps = signal.decimate(d.real, q, n, ftype=system,
                                        zero_phase=zero_phase)

            # Get phase from complex inner product, like CSD
            h_resamps = np.sum(d_tos.conj() * y_resamps, axis=-1)
            h_resamps /= np.abs(h_resamps)
            subnyq = freqs < 0.5*rate_to

            # Complex vectors should be aligned, only compare below nyquist
            assert_allclose(np.angle(h_resps.conj()*h_resamps)[subnyq], 0,
                            atol=1e-3, rtol=1e-3) 

def fast_xambg(refChannel, srvChannel, rangeBins, freqBins, shortFilt=True):
    ''' Fast Cross-Ambiguity Fuction (frequency domain method)
    
    Parameters:
        refChannel: reference channel data
        srvChannel: surveillance channel data
        rangeBins:  number of range bins to compute
        freqBins:   number of doppler bins to compute (should be power of 2)
        shortFilt:  (bool) chooses the type of decimation filter to use.
                    If True, uses an all-ones filter of length 1*(decimation factor)
                    If False, uses a flat-top window of length 10*(decimation factor)+1
    Returns:
        xambg: the cross-ambiguity surface. Dimensions are (nf, nlag+1, 1)
        third dimension added for easy stacking in dask

    '''
    if refChannel.shape != srvChannel.shape:
        raise ValueError('Input vectors must have the same length')

    # calculate decimation factor
    ndecim = int(refChannel.shape[0]/freqBins)

    # pre-allocate space for the result
    xambg = np.zeros((freqBins, rangeBins+1, 1), dtype=np.complex64)

    # complex conjugate of the second input vector
    srvChannelConj = np.conj(srvChannel)    

    if shortFilt:
        # precompute short FIR filter for decimation (all ones filter with length
        # equal to the decimation factor)
        dtaps = np.ones((ndecim + 1,))
    else:
        # precompute long FIR filter for decimation. (flat top filter of length
        # 10*decimation factor).  
        dtaps = signal.firwin(10*ndecim + 1, 1. / ndecim, window='flattop')

    dfilt = signal.dlti(dtaps, 1)

    # loop over range bins 
    for k, lag in enumerate(np.arange(-1*rangeBins, 1)):
        channelProduct = np.roll(srvChannelConj, lag)*refChannel
        #decimate the product of the reference channel and the delayed surveillance channel
        xambg[:,k,0] = signal.decimate(channelProduct, ndecim, ftype=dfilt)[0:freqBins]

    # take the FFT along the first axis (Doppler)
    # xambg = np.fft.fftshift(np.fft.fft(xambg, axis=0), axes=0)
    xambg = np.fft.fftshift(fft(xambg, axis=0), axes=0)
    return xambg 

def sinusoid_analysis(X, input_sample_rate, resample_block=128, copy=True):
    """
    Contruct a sinusoidal model for the input signal.

    Parameters
    ----------
    X : ndarray
        Input signal to model

    input_sample_rate : int
        The sample rate of the input signal

    resample_block : int, optional (default=128)
       Controls the step size of the sinusoidal model

    Returns
    -------
    frequencies_hz : ndarray
       Frequencies for the sinusoids, in Hz.

    magnitudes : ndarray
       Magnitudes of sinusoids returned in ``frequencies``

    References
    ----------
    D. P. W. Ellis (2004), "Sinewave Speech Analysis/Synthesis in Matlab",
    Web resource, available: http://www.ee.columbia.edu/ln/labrosa/matlab/sws/
    """
    X = np.array(X, copy=copy)
    resample_to = 8000
    if input_sample_rate != resample_to:
        if input_sample_rate % resample_to != 0:
            raise ValueError("Input sample rate must be a multiple of 8k!")
        # Should be able to use resample... ?
        # resampled_count = round(len(X) * resample_to / input_sample_rate)
        # X = sg.resample(X, resampled_count, window=sg.hanning(len(X)))
        X = sg.decimate(X, input_sample_rate // resample_to, zero_phase=True)
    step_size = 2 * round(resample_block / input_sample_rate * resample_to / 2.)
    a, g, e = lpc_analysis(X, order=8, window_step=step_size,
                           window_size=2 * step_size)
    f, m = lpc_to_frequency(a, g)
    f_hz = f * resample_to / (2 * np.pi)
    return f_hz, m 

def on_render(self):
        # get data from GUI
        downsample_fac = int(eval_expr(self.ui.downsample_fac.text()))
        win_size = int(eval_expr(self.ui.win_size.text()))
        fft_size = int(eval_expr(self.ui.fft_size.text()))
        overlap_fac = float(eval_expr(self.ui.overlap_fac.text()))
        clip_min, clip_max = float(eval_expr(self.ui.clip_min.text())), float(eval_expr(self.ui.clip_max.text()))
        x_tick_num, x_res = int(eval_expr(self.ui.x_num_ticks.text())), float(eval_expr(self.ui.x_resolution.text()))
        x_tick_rotation = int(eval_expr(self.ui.x_tick_rotation.text()))
        y_ticks_num, y_res = int(eval_expr(self.ui.y_num_ticks.text())), float(eval_expr(self.ui.y_resolution.text()))


        if downsample_fac > 1:
            downsampled_data = sig.decimate(self.data, downsample_fac, ftype='fir')
            downsampled_fs = self.fs / downsample_fac
        else:
            downsampled_data = self.data
            downsampled_fs = self.fs

        ft = Stft(downsampled_data, downsampled_fs, win_size=win_size, fft_size=fft_size, overlap_fac=overlap_fac)
        result = ft.stft(clip=(clip_min, clip_max))

        x_ticks, x_tick_labels = create_ticks_optimum(ft.freq_axis(), num_ticks=x_tick_num, resolution=x_res)
        y_ticks, y_tick_labels = create_ticks_optimum(ft.time_axis(), num_ticks=y_ticks_num, resolution=y_res)

        fig = self.ui.mpl.fig
        fig.clear()
        ax = fig.add_subplot(111)

        img = ax.imshow(result, origin='lower', cmap='jet', interpolation='none', aspect='auto')

        ax.set_xticks(x_ticks)
        ax.set_xticklabels(x_tick_labels, rotation=x_tick_rotation)
        ax.set_yticks(y_ticks)
        ax.set_yticklabels(y_tick_labels)

        if self.ui.x_grid.isChecked():
            ax.xaxis.grid(True, linestyle='-', linewidth=1)

        if self.ui.y_grid.isChecked():
            ax.yaxis.grid(True, linestyle='-', linewidth=1)

        ax.set_xlabel('Frequency [Hz]')
        ax.set_ylabel('Time [s]')

        fig.colorbar(img)
        fig.tight_layout()

        self.ui.mpl.draw()

        self.ui.sampling_freq.setText('%d' % downsampled_fs)
        self.ui.data_length.setText('%.2f' % ft.t_max)
        self.ui.freq_res.setText('%s' % (downsampled_fs * 0.5 / np.float32(ft.fft_size))) 

def pre_proc(dataset, database_name, is_debug=False):
	"""
	
	:param Mit2Segment list dataset: ECG segments.
	:return None:
	"""
	
	clear_id_set, noise_id_set = [], []
	for segment in dataset:
		if is_debug:
			print("now process %s	id=%s." % (segment.database_name, str(segment.global_id)))
		# denoising and write to txt file
		segment.denoised_ecg_data = denoise_ecg(segment.raw_ecg_data)
		segment.write_ecg_segment(rdf=1)
		# ecg data list to .mat
		list2mat(segment, is_debug=True)
		# compute RRI, RAMP and EDR.
		eng.computeFeatures(segment.base_file_path)
		if os.path.exists(segment.base_file_path + "/Rwave.mat"):
			RwaveMat = sio.loadmat(segment.base_file_path + "/Rwave.mat")
			Rwave = np.transpose(RwaveMat['Rwave'])
			Rwave = np.reshape(Rwave, len(Rwave))
			# RR intervals
			RR_intervals = np.diff(Rwave)
			# store RR intervals
			np.save(segment.base_file_path + "/RRI.npy", RR_intervals)
			# RRI validity check
			rri_flag = rricheck(segment.denoised_ecg_data, RR_intervals)
			if rri_flag:
				clear_id_set.append(segment.global_id)
			else:
				noise_id_set.append(segment.global_id)
				continue
			# compute R peaks amplitude(RAMP)
			RAMP = compute_r_peak_amplitude(segment.denoised_ecg_data, Rwave)
			# smoothing filtering
			RRI = smooth(RR_intervals, 3)
			RAMP = smooth(RAMP, 3)
			# spline interpolation
			RRI = RRI / ECG_RAW_FREQUENCY * 1000.0
			RRI = interp_cubic_spline(RRI, fs=4)
			RAMP = interp_cubic_spline_qrs(Rwave, RAMP, fs=4)
			# store RRI and RAMP
			np.save(segment.base_file_path + "/RRI.npy", RRI)
			np.save(segment.base_file_path + "/RAMP.npy", RAMP)
			# EDR
			EDRMat = sio.loadmat(segment.base_file_path + "/EDR.mat")
			EDR = np.transpose(EDRMat['EDR'])
			EDR = np.reshape(EDR, len(EDR))
			# downsampling
			EDR = decimate(EDR, 25)
			np.save(segment.base_file_path + "/EDR.npy", EDR)
			# print(".............")
		else:
			noise_id_set.append(segment.global_id)
	print(len(noise_id_set))
	print(len(clear_id_set))
	np.save(database_name[0] + "_" + database_name[1] + "_clear_id.npy", np.array(clear_id_set))
	np.save(database_name[0] + "_" + database_name[1] + "_noise_id.npy", np.array(noise_id_set)) 

def sinusoid_analysis(X, input_sample_rate, resample_block=128, copy=True):
    """
    Contruct a sinusoidal model for the input signal.

    Parameters
    ----------
    X : ndarray
        Input signal to model

    input_sample_rate : int
        The sample rate of the input signal

    resample_block : int, optional (default=128)
       Controls the step size of the sinusoidal model

    Returns
    -------
    frequencies_hz : ndarray
       Frequencies for the sinusoids, in Hz.

    magnitudes : ndarray
       Magnitudes of sinusoids returned in ``frequencies``

    References
    ----------
    D. P. W. Ellis (2004), "Sinewave Speech Analysis/Synthesis in Matlab",
    Web resource, available: http://www.ee.columbia.edu/ln/labrosa/matlab/sws/
    """
    X = np.array(X, copy=copy)
    resample_to = 8000
    if input_sample_rate != resample_to:
        if input_sample_rate % resample_to != 0:
            raise ValueError("Input sample rate must be a multiple of 8k!")
        # Should be able to use resample... ?
        # resampled_count = round(len(X) * resample_to / input_sample_rate)
        # X = sg.resample(X, resampled_count, window=sg.hanning(len(X)))
        X = sg.decimate(X, input_sample_rate // resample_to)
    step_size = 2 * round(resample_block / input_sample_rate * resample_to / 2.)
    a, g, e = lpc_analysis(X, order=8, window_step=step_size,
                           window_size=2 * step_size)
    f, m = lpc_to_frequency(a, g)
    f_hz = f * resample_to / (2 * np.pi)
    return f_hz, m 

def sinusoid_analysis(X, input_sample_rate, resample_block=128, copy=True):
    """
    Contruct a sinusoidal model for the input signal.

    Parameters
    ----------
    X : ndarray
        Input signal to model

    input_sample_rate : int
        The sample rate of the input signal

    resample_block : int, optional (default=128)
       Controls the step size of the sinusoidal model

    Returns
    -------
    frequencies_hz : ndarray
       Frequencies for the sinusoids, in Hz.

    magnitudes : ndarray
       Magnitudes of sinusoids returned in ``frequencies``

    References
    ----------
    D. P. W. Ellis (2004), "Sinewave Speech Analysis/Synthesis in Matlab",
    Web resource, available: http://www.ee.columbia.edu/ln/labrosa/matlab/sws/
    """
    X = np.array(X, copy=copy)
    resample_to = 8000
    if input_sample_rate != resample_to:
        if input_sample_rate % resample_to != 0:
            raise ValueError("Input sample rate must be a multiple of 8k!")
        # Should be able to use resample... ?
        # resampled_count = round(len(X) * resample_to / input_sample_rate)
        # X = sg.resample(X, resampled_count, window=sg.hanning(len(X)))
        X = sg.decimate(X, input_sample_rate // resample_to)
    step_size = 2 * round(resample_block / input_sample_rate * resample_to / 2.)
    a, g, e = lpc_analysis(X, order=8, window_step=step_size,
                           window_size=2 * step_size)
    f, m = lpc_to_frequency(a, g)
    f_hz = f * resample_to / (2 * np.pi)
    return f_hz, m 

def sinusoid_analysis(X, input_sample_rate, resample_block=128, copy=True):
    """
    Contruct a sinusoidal model for the input signal.

    Parameters
    ----------
    X : ndarray
        Input signal to model

    input_sample_rate : int
        The sample rate of the input signal

    resample_block : int, optional (default=128)
       Controls the step size of the sinusoidal model

    Returns
    -------
    frequencies_hz : ndarray
       Frequencies for the sinusoids, in Hz.

    magnitudes : ndarray
       Magnitudes of sinusoids returned in ``frequencies``

    References
    ----------
    D. P. W. Ellis (2004), "Sinewave Speech Analysis/Synthesis in Matlab",
    Web resource, available: http://www.ee.columbia.edu/ln/labrosa/matlab/sws/
    """
    X = np.array(X, copy=copy)
    resample_to = 8000
    if input_sample_rate != resample_to:
        if input_sample_rate % resample_to != 0:
            raise ValueError("Input sample rate must be a multiple of 8k!")
        # Should be able to use resample... ?
        # resampled_count = round(len(X) * resample_to / input_sample_rate)
        # X = sg.resample(X, resampled_count, window=sg.hanning(len(X)))
        X = sg.decimate(X, input_sample_rate // resample_to, zero_phase=True)
    step_size = 2 * round(resample_block / input_sample_rate * resample_to / 2.)
    a, g, e = lpc_analysis(X, order=8, window_step=step_size,
                           window_size=2 * step_size)
    f, m = lpc_to_frequency(a, g)
    f_hz = f * resample_to / (2 * np.pi)
    return f_hz, m 

def sinusoid_analysis(X, input_sample_rate, resample_block=128, copy=True):
    """
    Contruct a sinusoidal model for the input signal.

    Parameters
    ----------
    X : ndarray
        Input signal to model

    input_sample_rate : int
        The sample rate of the input signal

    resample_block : int, optional (default=128)
       Controls the step size of the sinusoidal model

    Returns
    -------
    frequencies_hz : ndarray
       Frequencies for the sinusoids, in Hz.

    magnitudes : ndarray
       Magnitudes of sinusoids returned in ``frequencies``

    References
    ----------
    D. P. W. Ellis (2004), "Sinewave Speech Analysis/Synthesis in Matlab",
    Web resource, available: http://www.ee.columbia.edu/ln/labrosa/matlab/sws/
    """
    X = np.array(X, copy=copy)
    resample_to = 8000
    if input_sample_rate != resample_to:
        if input_sample_rate % resample_to != 0:
            raise ValueError("Input sample rate must be a multiple of 8k!")
        # Should be able to use resample... ?
        # resampled_count = round(len(X) * resample_to / input_sample_rate)
        # X = sg.resample(X, resampled_count, window=sg.hanning(len(X)))
        X = sg.decimate(X, input_sample_rate // resample_to, zero_phase=True)
    step_size = 2 * round(resample_block / input_sample_rate * resample_to / 2.)
    a, g, e = lpc_analysis(X, order=8, window_step=step_size,
                           window_size=2 * step_size)
    f, m = lpc_to_frequency(a, g)
    f_hz = f * resample_to / (2 * np.pi)
    return f_hz, m 

def __call__(self, pkg):
        pkg = format_package(pkg)
        wav = pkg['chunk']
        tensor_wav = wav
        re_factor = self.sr // 8000
        wav = wav.data.numpy().reshape(-1).astype(np.float32)
        inwav = wav
        #wav = resample(wav, len(wav) // re_factor)
        #wav = librosa.core.resample(wav, self.sr, 8000,
        #                            res_type='kaiser_fast')
        wav = decimate(wav, self.sr // 8000)
        wav = np.array(wav * (2 ** 15), dtype=np.int16)
        total_frames = int(np.ceil(len(wav) / self.FRAME_SIZE))
        P_ = total_frames * self.FRAME_SIZE - len(wav)
        orilen = len(wav)
        if P_ > 0:
            wav = np.concatenate((wav, 
                                  np.zeros((P_,), dtype=np.int16)),
                                 axis=0)
        owav = []
        T = len(wav)
        data = [wav[t:t + self.FRAME_SIZE] for t in range(0, T,
                                                          self.FRAME_SIZE)]
        for frame in data:
            enc = self.c2.encode(frame)
            dec = self.c2.decode(enc)
            owav.extend(dec.tolist())
        owav = np.array(owav, dtype=np.int16)
        owav = owav[:orilen]
        #owav = np.array(owav, dtype=np.float32) / (2 ** 15)
        # resample up to original srate
        owav = resample(owav, len(owav) * re_factor)
        #owav = librosa.core.resample(owav, 8000, self.sr, res_type='kaiser_fast')
        owav = owav.astype(np.float32) / (2 ** 15)
        #owav = owav / (2 ** 15)
        owav = norm_energy(owav, inwav)
        if self.report:
            if 'report' not in pkg:
                pkg['report'] = {}
            pkg['report']['kbps'] = self.kbps
        #tensor_wav.data = torch.from_numpy(owav)
        pkg['chunk'] = torch.from_numpy(owav)
        return pkg