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_decimate():
    from numpy import arange,sin
    from scipy.signal import decimate, resample, cheby1, lfilter, filtfilt
    t = arange(0,30)
    q = 2
    n = 8
    print(t)
    print(decimate(t,2))
    print(resample(t, len(t)/2))
    t2 = sin(t)
    print(t2)
    print(len(t2))
    d2 = decimate(t2,2, ftype='fir')
    print(d2)
    print(len(d2))
    b, a = cheby1(n, 0.05, 0.8 / q)
    print(b,a)
    y = filtfilt(b, a, t2)
    print(y)
    sl = [slice(None)] * y.ndim

    sl[-1] = slice(None, None, -q)
    print(sl)
    print(y[sl])
    print(t2[sl])
    #r2 = resample(t2, len(t2)/2)
    #print(r2)
    #print(len(r2))
    assert(False)

def demodulate(args):
    samples = args[0]
    decim_r1 = args[1]
    decim_r2 = args[2]
    # DEMODULATION CODE

    # LIMITER goes here

    # low pass & down sampling
    lp_samples = signal.decimate(lowpass_filter(samples,32),int(decim_r1))


# polar discriminator

    A = lp_samples[1:lp_samples.size]
    B = lp_samples[0:lp_samples.size-1]

    dphase = ( A * np.conj(B) )

    dphase.resize(dphase.size+1)
    dphase[dphase.size-1] = dphase[dphase.size-2]

    rebuilt = signal.medfilt(np.angle(dphase)/np.pi,15) #  np.cos(dphase)

    output = signal.decimate(rebuilt,int(decim_r2))

    return np.real(output)

def downsampleMat(A, rowsBy=1, colsBy=1):
	if len(A.shape) == 1:
		return signal.decimate(A, rowsBy, n=(rowsBy-1))
	if rowsBy != 1:
		A = signal.decimate(A, rowsBy, n=(rowsBy-1), axis=0)
	if colsBy != 1:
		A = signal.decimate(A, rowsBy, n=(colsBy-1), axis=1)
	return A

def DownScaleSignal(signal_t, signal, scale):
 	"""
 	downscales the time and signal vectors
 	"""
	signal_d=SGN.decimate(signal,scale,ftype='fir')
	signal_t_d = SGN.decimate(signal_t,scale,ftype='fir')

	#pulse_plot(signal_t_d, signal_d)
	return signal_t_d, signal_d

def sf_anal(infile, chunk_rate=80.0, n_steps=12, base_freq=440.0, min_level=0.001, cutoff=None):
    min_sq_level = min_level**2
    if cutoff is None: cutoff = chunk_rate

    sr, wav = load_non_wav(infile)
    chunk_size = int(round(float(sr) / chunk_rate))
    wav = high_passed(sr, wav)
    wav2 = wav * wav
    freqs = 2**(np.linspace(0.0, n_steps, num=n_steps, endpoint=False)/n_steps) * base_freq

    amp2 = wav2
    rel_cutoff = 2.0/(float(sr)/cutoff)  # relative to nyquist freq, not samplerate
    b, a = RC(Wn=rel_cutoff)
    for i in xrange(4):
        amp2 = filtfilt(b, a, amp2)

    mask = amp2>min_sq_level
    little_amp2 = decimate(
        amp2, chunk_size, ftype='fir'
    )
    little_mask = mask[np.arange(0,little_amp2.size,1)*chunk_size]

    little_corrs = []
    for freq in freqs:
        # For now, offset is rounded to the nearest sample
        offset = int(round(float(sr)/freq))
        cov = np.zeros_like(wav)
        cov[:-offset] = wav[offset:]*wav[:-offset]
        # repeatedly filter; this is effectively an 8th-order lowpass now
        smooth_cov = cov
        for i in xrange(4):
            smooth_cov = filtfilt(b, a, smooth_cov)
        
        # technically the correlation should be taken wrt the harmonic mean of the variances at
        # the two times, but we assume autocorrelation lag << smooth time
        little_corrs.append(
            decimate(
                mask * smooth_cov/np.maximum(amp2, min_sq_level),
                chunk_size,
                ftype='fir' #FIR is needed to be stable at haptic rates
            )
        )
    
    
    all_corrs = np.vstack(little_corrs)
    sample_times = (np.arange(0,little_amp2.size,1)*chunk_size).astype(np.float)/sr

    #trim "too quiet" stuff
    all_corrs = all_corrs[:,np.where(little_mask)[0]]
    sample_times = sample_times[np.where(little_mask)[0]]
    little_amp2 = little_amp2[np.where(little_mask)[0]]
    
    return dict(
        all_corrs=all_corrs,
        sample_times=sample_times,
        amp=np.sqrt(little_amp2), #RMS amp is more usual
    )

def subsample_basis_data(data, oldFreq, newFreq, filterType='Default', filterOrder=None):
    if int((oldFreq / newFreq) % int(oldFreq / newFreq)) is not 0:
        raise ArithmeticError("Subsampling can be done only with integer factors of Old Frequency / New Frequency")
    if filterType == 'Default':
        subData = signal.decimate(data, int(oldFreq / newFreq))
    else:
        if filterOrder == None:
            raise ArithmeticError("A filter order is required if the filter is to be specified")
        subData = signal.decimate(data, int(oldFreq / newFreq), filterOrder, filterType)
    return subData

def load_continuous_tsd(paths, t_min=None, t_max=None, downsample=None,
                        columns=None):
    """
    read data for a specific time interval from a list of files (or ContinuousFile objects)
    Args:
        paths: a list of pathnames or ContinuousFile objects
        t_min: the low end of the time interval to read
        t_max: the high end of the time interval to read
        downsample: if not None, it should be an integer and acts as a downsampling factor (useful to get e.g. LFP)
        columns: a list of column names for the resulting TsdFrame. If None, the labels from the ContinuousFile objects 
        are used

    Returns:
        a TsdFrame with the data 
    """
    import scipy.signal as ss
    if isinstance(paths, str):
        paths = (paths,)
    elif not is_sequence(paths):
        raise TypeError("paths must be a string or list of strings.")

    if isinstance(paths[0], str):
        cf = [ContinuousFile(p) for p in paths]
    else:
        cf = paths

    data, tstamps = cf[0].read_interval(t_min, t_max)
    if downsample:
        data = ss.decimate(data, downsample, zero_phase=True)
    data = data.reshape((-1, 1))
    columns_from_files = False
    if columns is None:
        columns = [cf[0].label]
        columns_from_files = True
    if isinstance(columns, tuple):
        columns = list(columns)
    for f in cf[1:]:
        d, ts1 = f.read_interval(t_min, t_max)
        assert len(ts1) == len(tstamps)
        if downsample:
            d = ss.decimate(d, downsample, zero_phase=True)
        data = np.hstack((data, d.reshape((-1, 1))))
        if columns_from_files:
            columns.append(f.label)

    if downsample:
        tstamps = tstamps[::downsample]
        data = data[:, :len(tstamps)]

    cont_tsd = nts.TsdFrame(tstamps, data, columns=columns)

    return cont_tsd

    def detect_pitch_hps(self, data, sampling_rate):
        freq_spectrum = np.fft.rfft(data)
        freq_spectrum = np.absolute(freq_spectrum)
        fft1 = freq_spectrum / float(100000)
        fft2 = signal.decimate(fft1, 2)
        fft3 = signal.decimate(fft1, 3)
        fft4 = signal.decimate(fft1, 4)
        fft5 = signal.decimate(fft1, 5)
        fft6 = signal.decimate(fft1, 6)

        max1 = self.detect_max_idx(fft1)
        max2 = self.detect_max_idx(fft2)
        sec_max1 = self.detect_second_max(fft1)
        sec_max2 = self.detect_second_max(fft2)

        print "maxes"
        print max1 * float(sampling_rate) / data.size
        print max2

        print "second maxes"
        print sec_max1
        print sec_max2
        
        plt.plot(fft1)
        plt.plot(fft2)
        plt.plot(fft3)
        

        n = fft6.size
        hps = np.zeros(n)

        for i in range(n):
            hps[i] = fft1[i] * fft2[i] * fft3[i] * fft4[i] * fft5[i] * fft6[i]

        #plt.plot(hps)
        plt.axis([0, 300, 0, 30])
        plt.show()

        maxFreqIdx1 = 0
        maxFreqIdx2 = 0
        maxFreqAmp1 = 0
        maxFreqAmp2 = 0
        
        for i in range(n):
            if hps[i] > maxFreqAmp1:
                maxFreqIdx2 = maxFreqIdx1
                maxFreqAmp2 = maxFreqAmp1
                maxFreqIdx1 = i
                maxFreqAmp1 = hps[i]

        #N = data.size
        return sampling_rate / maxFreqIdx1

def loadITANfolder(folder,save_folder = None,q=25,overwrite=False):
    with Timer.Timer(folder):
        if type(save_folder) is 'NoneType':
           save_folder = folder  
        save_file = save_folder + folder[-14:-1] + '.nex'
        #get files in the folder
        files = getFileList(folder)
        # load info
        fid = open(folder+files['info'][0], 'rb')
        info = read_header(fid)
        sys.stdout.flush()
        #Sampling Rate
        sample_rate = info['sample_rate']
        time_vec = openDATfile(folder+files['time'][0],'time',sample_rate)
        time_vec = time_vec[0:-1:q]
        amp_unit = '$\mu V$'
        labels = []
        nch = len(files['amp']) + len(files['adc']) # +len(files['aux'])
        data = np.zeros([time_vec.shape[0],nch])
        eng = matlab.engine.start_matlab()
        eng.cd(folder,nargout=0)
        count = 0
        for f in files['amp']:
            sys.stdout.flush()
            with Timer.Timer(f):     
                name = f[:-4]
                labels.append(name)
                aux_data = openDATfile(folder+f,'amp',sample_rate)
                data[:,count] = sig.decimate(aux_data,q)
                count +=1
                if not overwrite:
                    if os.path.isfile(folder+name+'.mat'):
                        continue
                tfile = open(folder + 'Files.txt', 'w')
                tfile.write(name +'\n')
                tfile.close()
                sio.savemat(folder+name+'.mat', {'data':aux_data})
                eng.Get_spikes_alt(sample_rate,nargout=0)
                eng.close('all', nargout=0)
        eng.save_NEX(sample_rate,labels,int(time_vec[0]),save_file,nargout=0)
        for f in files['adc']:
            sys.stdout.flush()
            with Timer.Timer(f):  
                labels.append(f[:-4])
                aux_data = openDATfile(folder+f,'adc',sample_rate)
                data[:,count] = sig.decimate(aux_data,q)
                count +=1
        Data = DataObj.DataObj(data,sample_rate/q,amp_unit,labels,time_vec,[])
    Data.save(folder+'downsampled.h5','data')
    eng.quit()    
    return Data

def load_downsample_and_dump(base_dir, N, movie_name, arr_name):
    MOV_ARR_DIRS= [os.path.join(base_dir, str(i)) for i in range(N)]
    movies = [np.load(os.path.join(d, arr_name + '.npy')) for d in MOV_ARR_DIRS]
    for i_p in range(5):
        p = i_p + 1
        print 'Downsampling factor %d:' % 2**p
        for i, movie in enumerate(movies):
            m = decimate(decimate(movie, 2**p, axis=0), 2**p, axis=1)

            print 'Saving %s %s %d...' % (movie_name, arr_name, i)
            if not os.path.isdir(os.path.join(base_dir,str(i),
                                              arr_name+'_down')):
                os.makedirs(os.path.join(base_dir,str(i),arr_name+'_down'))
            np.save(os.path.join(base_dir,str(i),arr_name+'_down/%d'%2**p),m)
            print 'Saved.'

def decimate_labeled ( x, q, **kwargs ):
    '''Labeled analogue of scipy.signal.decimate; runs along x's 'time' axis.

    Inputs:
    x - LabeledArray containing the data; must have a 'time' axis
    q - Decimation factor passed to decimate
    **kwargs - (Optional) Keyword arguments passed to decimate

    Output is a LabeledArray with the same axes as x except for 'time', which
    is subsampled accordingly.
    '''
    
    # Make sure we don't override axis in the kwargs
    kwargs.pop( 'axis', None )

    time_axis = x.axis_index( 'time' )
    
    # Decimate
    ret_array = sig.decimate( x.array, q,
                              axis = time_axis,
                              **kwargs )
    
    # Form new axes
    ret_axes = OrderedDict( x.axes )
    ret_time = np.linspace( x.axes['time'][0], x.axes['time'][-1], ret_array.shape[time_axis] )
    ret_axes['time'] = ret_time
    
    return LabeledArray( array = ret_array, axes = ret_axes )

def process_wav_files(wav_file_paths = example_wav_paths ,
                      output_dir = example_seg_paths,
                      savefig = False):

    for wav_file_path in wav_file_paths:

        # import data and down sample
        raw_sample_rate, raw_wav_data = wavfile.read( wav_file_path )
        downsample_factor = 4
        wav_data = decimate(raw_wav_data, downsample_factor)
        sample_rate = raw_sample_rate/downsample_factor

        # segment data (obtain segments and figure)
        segments, f, ax = simple_segmentation(wav_data, sample_rate)

        # read attributes from xml
        mediaID, classID = xml_attributes(wav_file_path.replace('.wav','.xml'))
 
        # add title and labels to plot
        ax.set_title(" MediaID: {0} ClassID: {1} ".format(mediaID, classID))
        ax.set_xlabel("Time (seconds)")
        ax.set_ylabel("Frequency (Hz)")

        basename = os.path.basename(wav_file_path.replace('.wav',''))

        if savefig:
            # save figure to output dir
            f.savefig(os.path.join(output_dir, basename+".pdf"))
        
        # save segments to csv
        np.savetxt(os.path.join(output_dir, basename+".csv"),
                   segments, delimiter = ",", fmt = '%3.4f')

    return None

    def DistortionProducts(self, signal, reverbShift=1., reverbTau=1., echos=5, distortFactor=4., bandpass=True):
        """ Generate distortions to reverberated signal.
        :param signal: Signal
        :type signal: numpy arrray
        :param reverbShift: Delay between reverberations
        :type reverbShift: float
        :param reverbTau: Exponential decay constant
        :type reverbTau: float
        :param echos:  Number of reverberations (~5 is usually sufficient)
        :type echos: int
        :param distortFactor: Inverse scaling factor for signal in Boltzman filter (large value is small distortion)
        :type distortFactor: float
        :returns: numpy array with distorted signal
        """
        shift = reverbShift
        tau = reverbTau
        N = echos
        sigReverb = self.signalReverb(signal, shift, tau, N)

        resampleFactor = 4.
        sigReverb = resample(sigReverb, resampleFactor*len(sigReverb))
        self.Fs = resampleFactor*self.Fs

        sigReverbD = self.boltz(sigReverb/distortFactor)
        sigReverbDn = (sigReverbD - np.mean(sigReverbD[:100]))/np.max(np.abs(sigReverbD - np.mean(sigReverbD[:100])))
        if bandpass: 
        	sigReverbD_F = self.BandPassFilter(sigReverbDn, set_numtaps=501)
        else: sigReverbD_F = sigReverbDn

        resampleFactor = 4
        sigReverbD_F = decimate(sigReverbD_F, resampleFactor)
        self.Fs = self.Fs/float(resampleFactor)
        return sigReverbD_F

def compute_model_ts(x, y, sigma, hrf_delay, theta, phi, cpd,
                     deg_x, deg_y, stim_arr, tr_length, 
                     frames_per_tr, norm_func=utils.zscore):
    
    # otherwise generate a prediction
    ts_stim = MakeFastGaborPrediction(deg_x,
                                      deg_y,
                                      stim_arr,
                                      x,
                                      y,
                                      sigma,
                                      theta,
                                      phi,
                                      cpd)
    
    # convolve it
    hrf = utils.double_gamma_hrf(hrf_delay, tr_length, frames_per_tr)
    
    # normalize it
    model = norm_func(ss.fftconvolve(ts_stim, hrf)[0:len(ts_stim)])

    # decimate it
    model = ss.decimate(model, int(frames_per_tr), 1)
    
    return model

def down_scale_signal_(signal, scale):
 	"""
 	downscales the signal vector. Re-scale the energy
 	"""
	signal_d=SGN.decimate(signal,scale,ftype='fir')
	
	return signal_d*scale

def hps(x, fs=44100, lf=255, harmonics=3, precision=2, window=lambda l:_np.kaiser(l, 7.14285)):
    """ Estimates the pitch (fundamental frequency) of the given sample array by a standard HPS implementation. """
    x -= _np.mean(x)
    N = x.size
    w = x*window(N)

    # Append zeros to the end of the window so that each bin has at least the desired precision.
    if fs/N > precision:
        delta = int(fs/precision) - N
        w = _np.append(w, _np.zeros(delta))
        N = w.size

    X = _np.log(_np.abs(_np.fft.rfft(w)))

    # Sequentially decimate 'X' 'harmonics' times and add it to itself.
    # 'precision < fs/N' must hold, lest the decimation loses all the precision we'd gain.
    hps = _np.copy(X)
    for h in range(2, 2 + harmonics):
        dec = _sig.decimate(X, h)
        hps[:dec.size] += dec*(0.8**h)

    # Find the bin corresponding to the lowest detectable frequency.
    lb = lf*N/fs

    # And then the bin with the highest spectral content.
    arg_peak = lb + _np.argmax(hps[lb:dec.size])

    # TODO: Return the full array? A ranked list of identified notes?
    return fs*arg_peak/N

def audio_to_pitch_given_boundaries(audio, filter_bank, window_boundaries):
    """
    like audio_to_pitch, but takes an iterable (window_boundaries) that contains
    2-tuples of the start and end of each window.
    """
    output = {} # output will be a dictionary indexed by pitch.
    for ratio, filters in filter_bank.items():
        if ratio == 1.0:
            resampled_audio = audio # decimate with ftype='fir' bitches when the ratio is 1. So do nothing in this case.
        else:
            resampled_audio = sig.decimate(audio, ratio, ftype='fir') # scipy's resample is different from matlab's resample...
            # use ftype='fir' instead of the default 'iir' to avoid warnings and bogus values.
        
        for pitch, filter_coeffs in filters.items():
            coeff_b, coeff_a = filter_coeffs
            out = sig.filtfilt(coeff_b, coeff_a, resampled_audio)
            out_squared = out ** 2 # energy
            
            # summarise over windows:
            summarised_output = np.zeros(len(window_boundaries))
            for (i, (start, end)) in enumerate(window_boundaries):
                summarised_output[i] = np.mean(out_squared[int(start/ratio):int(end/ratio)])
            
            output[pitch] = summarised_output

    return output

def freq_from_hps(signal, fs):
    """
    Estimate frequency using harmonic product spectrum

    Low frequency noise piles up and overwhelms the desired peaks

    Doesn't work well if signal doesn't have harmonics
    """
    signal = asarray(signal) + 0.0

    N = len(signal)
    signal -= mean(signal)  # Remove DC offset

    # Compute Fourier transform of windowed signal
    windowed = signal * kaiser(N, 100)

    # Get spectrum
    X = log(abs(rfft(windowed)))

    # Remove mean of spectrum (so sum is not increasingly offset
    # only in overlap region)
    X -= mean(X)

    # Downsample sum logs of spectra instead of multiplying
    hps = copy(X)
    for h in range(2, 9):  # TODO: choose a smarter upper limit
        dec = decimate(X, h, zero_phase=True)
        hps[:len(dec)] += dec

    # Find the peak and interpolate to get a more accurate peak
    i_peak = argmax(hps[:len(dec)])
    i_interp = parabolic(hps, i_peak)[0]

    # Convert to equivalent frequency
    return fs * i_interp / N  # Hz

def lofar(data, fs, n_pts_fft=1024, n_overlap=0,
    decimation_rate=3, spectrum_bins_left=None, **tpsw_args):
    norm_parameters = {'lat_window_size': 10, 'lat_gap_size': 1, 'threshold': 1.3}
    if not isinstance(data, np.ndarray):
        raise NotImplementedError

    if decimation_rate > 1:
        dec_data = decimate(data, decimation_rate, 10, 'fir', zero_phase=True)
        Fs = fs/decimation_rate
    else:
        dec_data = data
        Fs=fs
    freq, time, power = spectrogram(dec_data,
                                    window=('hann'),
                                    nperseg=n_pts_fft,
                                    noverlap=n_overlap,
                                    nfft=n_pts_fft,
                                    fs=Fs,
                                    detrend=False,
                                    axis=0,
                                    scaling='spectrum',
                                    mode='magnitude')
    power = np.absolute(power)
    power = power / tpsw(power)#, **tpsw_args)
    power = np.log10(power)
    power[power < -0.2] = 0

    if spectrum_bins_left is None:
        spectrum_bins_left = power.shape[0]*0.8
    power = power[:spectrum_bins_left, :]
    freq = freq[:spectrum_bins_left]

    return np.transpose(power), freq, time