我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用scipy.signal.lfilter()。
def butter_highpass(highcut, fs, order): nyq = 0.5 * fs high = highcut / nyq b, a = butter(order, high, btype="highpass") return b, a # Sources for Batch Iterators # # These classes load training and test data and perform some basic preprocessing on it. # They are then passed to factory functions that create the net. There they are used # as data sources for the batch iterators that feed data to the net. # All classes band pass or low pass filter their data based on min / max freq using # a causal filter (lfilter) when the data is first loaded. # * TrainSource loads a several series of EEG data and events, splices them together into # one long stream, then normalizes the EEG data to zero mean and unit standard deviation. # * TestSource is like TrainSource except that it uses the mean and standard deviation # computed for the associated training source to normalize the EEG data. # * SubmitSource is like TestSource except that it does not load and event data.
def phormants(x, Fs): N = len(x) w = numpy.hamming(N) # Apply window and high pass filter. x1 = x * w x1 = lfilter([1], [1., 0.63], x1) # Get LPC. ncoeff = 2 + Fs / 1000 A, e, k = lpc(x1, ncoeff) #A, e, k = lpc(x1, 8) # Get roots. rts = numpy.roots(A) rts = [r for r in rts if numpy.imag(r) >= 0] # Get angles. angz = numpy.arctan2(numpy.imag(rts), numpy.real(rts)) # Get frequencies. frqs = sorted(angz * (Fs / (2 * math.pi))) return frqs
def sp_filter_butterworth_bandpass(input_data, passband, stopband, max_passband_loss, min_stopband_attenuation): # The SciPy signal processing module uses normalised frequencies, so we need to normalise the input values norm_passband = input_data.timebase.normalise_freq(passband) norm_stopband = input_data.timebase.normalise_freq(stopband) ord,wn = sp_signal.filter_design.buttord(norm_passband, norm_stopband, max_passband_loss, min_stopband_attenuation) b, a = sp_signal.filter_design.butter(ord, wn, btype = 'bandpass') output_data = input_data for i,s in enumerate(output_data.signal): output_data.signal[i] = sp_signal.lfilter(b,a,s) return output_data ######################################### ## wrappers to numpy signal processing ## #########################################
def preemphasis(x, coef=0.97): """Pre-emphasis Args: x (1d-array): Input signal. coef (float): Pre-emphasis coefficient. Returns: array: Output filtered signal. See also: :func:`inv_preemphasis` Examples: >>> from nnmnkwii.util import example_audio_file >>> from scipy.io import wavfile >>> fs, x = wavfile.read(example_audio_file()) >>> x = x.astype(np.float64) >>> from nnmnkwii import preprocessing as P >>> y = P.preemphasis(x, coef=0.97) >>> assert x.shape == y.shape """ b = np.array([1., -coef], x.dtype) a = np.array([1.], x.dtype) return signal.lfilter(b, a, x)
def Gaussian2D(image, sigma, padding=0): n, m = image.shape[0], image.shape[1] tmp = np.zeros((n + padding, m + padding)) if tmp.shape[0] < 4: raise ValueError('Image and padding too small') if tmp.shape[1] < 4: raise ValueError('Image and padding too small') B, A = __gausscoeff(sigma) tmp[:n, :m] = image tmp = lfilter(B, A, tmp, axis=0) tmp = np.flipud(tmp) tmp = lfilter(B, A, tmp, axis=0) tmp = np.flipud(tmp) tmp = lfilter(B, A, tmp, axis=1) tmp = np.fliplr(tmp) tmp = lfilter(B, A, tmp, axis=1) tmp = np.fliplr(tmp) return tmp[:n, :m] #-----------------------------------------------------------------------------
def load_series(self, subject, series): min_freq = self.min_freq max_freq = self.max_freq key = (subject, series, min_freq, max_freq) if key not in self._series_cache: while len(self._series_cache) > self.MAX_CACHE_SIZE: # Randomly throw away an item self._series_cache.popitem() print("Loading", subject, series) data = read_csv(path(subject, series, "data")) # Filter here since it's slow and we don't want to filter multiple # times. `lfilter` is CAUSAL and thus doesn't violate the ban on future data. if (self.min_freq is None) or (self.min_freq == 0): print("Low pass filtering, f_h =", max_freq) b, a = butter_lowpass(max_freq, SAMPLE_RATE, FILTER_N) else: print("Band pass filtering, f_l =", min_freq, "f_h =", max_freq) b, a = butter_bandpass(min_freq, max_freq, SAMPLE_RATE, FILTER_N) self._series_cache[key] = lfilter(b, a, data, axis=0) return self._series_cache[key]
def _apply_filter(sys, dt, u): # "Correct" implementation of filt that has a single time-step delay # see Nengo issue #938 if dt is not None: num, den = cont2discrete(sys, dt).tf elif not sys.analog: num, den = sys.tf else: raise ValueError("system (%s) must be discrete if not given dt" % sys) # convert from the polynomial representation, and add back the leading # zeros that were dropped by poly1d, since lfilter will shift it the # wrong way (it will add the leading zeros back to the end, effectively # removing the delay) num, den = map(np.asarray, (num, den)) num = np.append([0]*(len(den) - len(num)), num) return lfilter(num, den, u, axis=-1)
def _process(self): # Receive input data. batch = self.input.receive() self._measure_time('start', frequency=100) # Process data. for i in xrange(self.nb_channels): batch[:, i], self.z[i] = signal.lfilter(self.b, self.a, batch[:, i], zi=self.z[i]) batch[:, i] -= np.median(batch[:, i]) if self.remove_median: global_median = np.median(batch, 1) for i in xrange(self.nb_channels): batch[:, i] -= global_median # Send output data. self.output.send(batch) self._measure_time('end', frequency=100) return
def _generate_noise(info, cov, iir_filter, random_state, n_samples, zi=None): """Helper to create spatially colored and temporally IIR-filtered noise""" from scipy.signal import lfilter noise_cov = pick_channels_cov(cov, include=info['ch_names'], exclude=[]) rng = check_random_state(random_state) c = np.diag(noise_cov.data) if noise_cov['diag'] else noise_cov.data mu_channels = np.zeros(len(c)) # we almost always get a positive semidefinite warning here, so squash it with warnings.catch_warnings(record=True): noise = rng.multivariate_normal(mu_channels, c, n_samples).T if iir_filter is not None: if zi is None: zi = np.zeros((len(c), len(iir_filter) - 1)) noise, zf = lfilter([1], iir_filter, noise, axis=-1, zi=zi) else: zf = None return noise, zf
def __call__(self, x): y, self.filter_state = lfilter(b=self.b, a=self.a, x = np.array(x)[None], zi = self.filter_state) return y[0]
def butter_bandpass_filter(data, highpass, fs, order=4): b, a = butter_bandpass(0, highpass, fs, order=order) y = lfilter(b, a, data) return y
def discount(x, gamma): """ Compute discounted sum of future values out[i] = in[i] + gamma * in[i+1] + gamma^2 * in[i+2] + ... """ return signal.lfilter([1],[1,-gamma],x[::-1], axis=0)[::-1]
def butter_bandpass_filter(data, lowcut, highcut, order=5): b, a = butter_bandpass(lowcut, highcut, order=order) filtered_data = lfilter(b, a, data) return filtered_data
def inv_preemphasis(x, coef=0.97): """Inverse operation of pre-emphasis Args: x (1d-array): Input signal. coef (float): Pre-emphasis coefficient. Returns: array: Output filtered signal. See also: :func:`preemphasis` Examples: >>> from nnmnkwii.util import example_audio_file >>> from scipy.io import wavfile >>> fs, x = wavfile.read(example_audio_file()) >>> x = x.astype(np.float64) >>> from nnmnkwii import preprocessing as P >>> x_hat = P.inv_preemphasis(P.preemphasis(x, coef=0.97), coef=0.97) >>> assert np.allclose(x, x_hat) """ b = np.array([1.], x.dtype) a = np.array([1., -coef], x.dtype) return signal.lfilter(b, a, x)
def Filter(self, LowCorner, HighCorner, Order=3): """ Butterworth bandpass filter """ FS = 1./self.HDR['TSMP'] if HighCorner >= FS/2.: print 'Warning: High corner must be < {0:.2f} Hz'.format(FS/2.) return if LowCorner < 0.: print 'Warning: Low corner must be > 0 Hz'.format(FS/2.) return # Corner frequencies Corners = [2.*LowCorner/FS, 2.*HighCorner/FS] # Butterworth filter b, a = _sig.butter(Order, Corners, btype='band') # Filtering records for I,S in enumerate(self.CHN): # self.CHN[I] = _sig.lfilter(b, a, S) zi = _sig.lfilter_zi(b, a); self.CHN[I],_ = _sig.lfilter(b, a, S, zi=zi*S[0]) #-----------------------------------------------------------------------------------------
def erb_frilter_bank(x, fcoefs): a0 = fcoefs[:, 0] a11 = fcoefs[:, 1] a12 = fcoefs[:, 2] a13 = fcoefs[:, 3] a14 = fcoefs[:, 4] a2 = fcoefs[:, 5] b0 = fcoefs[:, 6] b1 = fcoefs[:, 7] b2 = fcoefs[:, 8] gain = fcoefs[:, 9] output = np.zeros((np.size(gain, 0), np.size(x, 0))) for chan in range(np.size(gain, 0)): y1 = lfilter(np.array([a0[chan] / gain[chan], a11[chan] / gain[chan], a2[chan] / gain[chan]]), np.array([b0[chan], b1[chan], b2[chan]]), x) y2 = lfilter(np.array([a0[chan], a12[chan], a2[chan]]), np.array([b0[chan], b1[chan], b2[chan]]), y1) y3 = lfilter(np.array([a0[chan], a13[chan], a2[chan]]), np.array([b0[chan], b1[chan], b2[chan]]), y2) y4 = lfilter(np.array([a0[chan], a14[chan], a2[chan]]), np.array([b0[chan], b1[chan], b2[chan]]), y3) output[chan, :] = y4 return output
def rasta_filt(x): number = np.arange(-2., 3., 1.) number = -1. * number / np.sum(number*number) denom = np.array([1., -0.94]) zi = lfilter_zi(number, 1) zi = zi.reshape(1, len(zi)) zi = np.repeat(zi, np.size(x, 0), 0) y, zf = lfilter(number, 1, x[:,0:4], axis=1, zi=zi) y, zf = lfilter(number, denom, x, axis=1, zi=zf) return y
def _analog_filter(self, ftype, highpass=None, lowpass=None, order=3, zerophase=True): ' Use a classic analog filter on the data, currently butter or bessel' from scipy.signal import butter, bessel filter_types = {'butter': butter, 'bessel': bessel} afilter = filter_types[ftype] if highpass is None and lowpass is not None: b, a = afilter(order, lowpass / (self.sr / 2), btype='lowpass') elif highpass is not None and lowpass is None: b, a = afilter(order, highpass / (self.sr / 2), btype='highpass') elif highpass is not None and lowpass is not None: if highpass < lowpass: b, a = afilter(order, (highpass / (self.sr / 2), lowpass / (self.sr / 2)), btype='bandpass') else: b, a = afilter(order, (lowpass / (self.sr / 2), highpass / (self.sr / 2)), btype='bandstop') if zerophase: return self.filtfilt(b, a) else: return self.lfilter(b, a)
def lfilter(self, b, a): " Forward only filtering" from scipy.signal import lfilter def filter_func(x): return lfilter(b, a, x, axis=0) return self.new_stream(self.vector_map(filter_func))
def low_pass_filter(self, sig, low_cut, high_cut, nyq): """ Apply a low pass filter to the data to take the relevant frequencies and to smooth the drop-out regions No resample necessary because all files have the same sampling frequency """ b, a = signal.butter(5, np.array([low_cut, high_cut]) / nyq, btype='band') sig_filt = signal.lfilter(b, a, sig, axis=0) return(np.float32(sig_filt))
def butter_bandpass_filter(data, lowcut, highcut, fs, order): from scipy.signal import butter, lfilter nyq = 0.5 * fs low = lowcut / nyq high = highcut / nyq b, a = butter(order, [low, high], btype='bandpass') y = lfilter(b, a, data) return y
def butter_bandstop_filter(data, lowcut, highcut, fs, order): from scipy.signal import butter, lfilter nyq = 0.5 * fs low = lowcut / nyq high = highcut / nyq b, a = butter(order, [low, high], btype='bandstop') y = lfilter(b, a, data) return y
def butter_lowpass_filter(data, cut, fs, order, zero_phase=False): from scipy.signal import butter, lfilter, filtfilt nyq = 0.5 * fs cut = cut / nyq b, a = butter(order, cut, btype='low') y = (filtfilt if zero_phase else lfilter)(b, a, data) return y
def apply_ellip_filter(self, data, order, cutoff, btype): ripple = 3 attenuation = 50 cutoff /= self.nyq b, a = iirfilter(N=order, Wn=cutoff, rp=ripple, rs=attenuation, btype=btype, ftype='ellip') return lfilter(b, a, data, axis=0)
def apply(self, data): nyq = self.f / 2.0 cutoff = min(self.f, nyq-1) h = signal.firwin(numtaps=101, cutoff=cutoff, nyq=nyq) # data[i][ch][dim0] for i in range(len(data)): data_point = data[i] for j in range(len(data_point)): data_point[j] = signal.lfilter(h, 1.0, data_point[j]) return data
def _discount(self, x): a = np.asarray(x) return lfilter([1], [1, -self.discount_factor], a[::-1], axis=0)[::-1]
def discount(self, x, gamma): """ computes discounted sums along 0th dimension of x. inputs ------ x: ndarray gamma: float outputs ------- y: ndarray with same shape as x, satisfying y[t] = x[t] + gamma*x[t+1] + gamma^2*x[t+2] + ... + gamma^k x[t+k], where k = len(x) - t - 1 """ x = np.array(x) assert x.ndim >= 1 return lfilter([1], [1, -gamma], x[::-1], axis=0)[::-1]
def butter_bandpass_filter(data, *, lowcut, highcut, fs, order=5): """Band filter data using a butterworth filter.""" b, a = butter_bandpass(lowcut, highcut, fs, order=order) y = lfilter(b, a, data) return y
def test_lfilter_gen(self): x = np.random.normal(0, 1, 1000) hb = np.array([0, 0, 1, 0], dtype=np.float) f = signal.lfilter_gen(hb, 1) y = [f.send(v) for v in x] self.assertArrayEqual(np.append([0, 0], x[:-2]), y) hb, ha = sp.iirfilter(4, 0.01, btype='lowpass') y1 = sp.lfilter(hb, ha, x) f = signal.lfilter_gen(hb, ha) y2 = [f.send(v) for v in x] self.assertArrayEqual(y1, y2, precision=6)
def mfilter(s, x, axis=0, complex_output=False): """Matched filter recevied signal using a reference signal. :param s: reference signal :param x: recevied signal :param axis: axis of the signal, if multiple signals specified :param complex_output: True to return complex signal, False for absolute value of complex signal """ hb = _np.conj(_np.flipud(s)) x = _np.pad(x, (0, len(s)-1), 'constant') y = _sig.lfilter(hb, 1, x, axis)[len(s)-1:] if not complex_output: y = _np.abs(y) return y
def lfilter0(b, a, x, axis=0): """Filter data with an IIR or FIR filter with zero DC group delay. :func:`scipy.signal.lfilter` provides a way to filter a signal `x` using a FIR/IIR filter defined by `b` and `a`. The resulting output is delayed, as compared to the input by the group delay. This function corrects for the group delay, resulting in an output that is synchronized with the input signal. If the filter as an acausal impulse response, some precursor signal from the output will be lost. To avoid this, pad input signal `x` with sufficient zeros at the beginning to capture the precursor. Since both, :func:`scipy.signal.lfilter` and this function return a signal with the same length as the input, some signal tail is lost at the end. To avoid this, pad input signal `x` with sufficient zeros at the end. See documentation for :func:`scipy.signal.lfilter` for more details. >>> import arlpy >>> import numpy as np >>> fs = 250000 >>> b = arlpy.uwa.absorption_filter(fs, distance=500) >>> x = np.pad(arlpy.signal.sweep(20000, 40000, 0.5, fs), (127, 127), 'constant') >>> y = arlpy.signal.lfilter0(b, 1, x) """ w, g = _sig.group_delay((b, a)) ndx = _np.argmin(_np.abs(w)) d = int(round(g[ndx])) x = _np.pad(x, (0, d), 'constant') y = _sig.lfilter(b, a, x, axis)[d:] return y
def harvest_filter_f0_contour(f0, st, ed, b, a): smoothed_f0 = f0.copy() smoothed_f0[:st] = smoothed_f0[st] smoothed_f0[ed + 1:] = smoothed_f0[ed] aaa = sg.lfilter(b, a, smoothed_f0) bbb = sg.lfilter(b, a, aaa[::-1]) smoothed_f0 = bbb[::-1].copy() smoothed_f0[:st] = 0. smoothed_f0[ed + 1:] = 0. return smoothed_f0
def butter_lowpass_filter(data, cutoff, fs, order=2): """ This function designs a Butterworth lowpass filter and applies it to the data :param data: Audio data to be lowpass filtered :param cutoff: Cutoff frequency in Hz :param fs: Samplerate of audio :param order: Order of the filter, defaults to second order filter as required by timbral_metallic :return: Returns the filtered signal """ b, a = butter_lowpass(cutoff, fs, order=order) y = lfilter(b, a, data) return y
def discount_rewards(x, gamma): """ Given vector x, computes a vector y such that y[i] = x[i] + gamma * x[i+1] + gamma^2 x[i+2] + ... """ return signal.lfilter([1], [1, -gamma], x[::-1], axis=0)[::-1] # Source: http://stackoverflow.com/a/12201744/1735784
def SPBFILTER(df, n1 = 40, n2 = 60, n3 = 0, field = 'close'): if n3 == 0: n3 = int((n1 + n2)/2) a1 = 5.0/n1 a2 = 5.0/n2 B = [a1-a2, a2-a1] A = [1, (1-a1)+(1-a2), -(1-a1)*(1-a2)] PB = pd.Series(signal.lfilter(B, A, df[field]), name = 'SPB_%s_%s' % (n1, n2)) RMS = pd.Series(pd.rolling_mean(PB*PB, n3)**0.5, name = 'SPBRMS__%s_%s' % (n1, n2)) return pd.concat([PB, RMS], join='outer', axis=1)
def Fdeltas(x,w=9): ''' # D = deltas(X,W) Calculate the deltas (derivatives) of a sequence input feature*frame ''' # nr feature, nc frame nr,nc = x.shape # Define window shape hlen = int(np.floor(w/2)) w = 2*hlen + 1 win = np.arange(hlen,-hlen-1,-1) # normalisation win = win/np.sum(np.abs(win)) # pad data by repeating first and last columns xx = np.hstack((np.tile(x[:,0],(hlen,1)).transpose(),x,np.tile(x[:,-1],(hlen,1)).transpose())) # Apply the delta filter d = lfilter(win, 1, xx, 1) # filter along dim 1 (rows) # Trim edges d = d[:,2*hlen:2*hlen+nc] return d
def filter(self, s, axis=0): """Filter new data. """ if self.bandType == 'allpass': return s if self.bandType == 'allstop': return np.zeros_like(s) """ Should be very close to filtfilt, padding? XXX - idfah if self.zeroPhase: rev = [slice(None),]*s.ndim rev[axis] = slice(None,None,-1) #ziScaled = self.scaleZi(s[rev], axis) y, newZi = spsig.lfilter(self.numCoef, self.denomCoef, s[rev], axis=axis, zi=newZi) y = y[rev] """ # if zeroPhase and signal is shorter than padlen (default in filtfilt function) if self.zeroPhase and \ (3*max(len(self.numCoef), len(self.denomCoef))) < s.shape[axis]: # need astype below since filtfilt calls lfilter_zi, which does not preserve dtype XXX - idfah return spsig.filtfilt(self.numCoef, self.denomCoef, s, axis=axis, padtype='even').astype(self.dtype, copy=False) else: ziScaled = self.scaleZi(s, axis) # even padding to help reduce edge effects nPad = 3*max(len(self.numCoef), len(self.denomCoef)) sPad = np.apply_along_axis(np.pad, axis, s, pad_width=nPad, mode='reflect') # edge for constant padding slc = [slice(nPad,-nPad) if i == axis else slice(None) for i in range(s.ndim)] y, newZi = spsig.lfilter(self.numCoef, self.denomCoef, sPad, axis=axis, zi=ziScaled) return y[slc]
def filter(self, s, axis=0): if self.bandType == 'allpass': return s if self.bandType == 'allstop': return np.zeros_like(s) if not self.ziSaved: # use padding to find initial zi? XXX - idfah self.zi = self.scaleZi(s, axis) self.ziSaved = True y, self.zi = spsig.lfilter(nc, dc, s, axis=axis, zi=self.zi) return y
def butter_highpass_filter(data, cutoff, fs, order=5): b, a = butter_highpass(cutoff, fs, order=order) y = lfilter(b, a, data) return y
def butter_highpass_filter(data, cutoff, fs, order=5): b, a = AudioAnalyzer.butter_highpass(cutoff, fs, order=order) y = signal.lfilter(b, a, data) return y
def sosfilter(sos, zi_in, x): y = x zi_out = zi_in for i in range(len(sos)): y, zi_out[i] = lfilter(sos[i,:3], sos[i,3:], y, zi = zi_in[i]) return y, zi_out
def filter(self): path = os.getcwd()+'/trialGraspEventDetection_dataFiles' self.Fgr = np.sum(self.values[:, 9:15], axis=1) # SAI self.Fgl = np.sum(self.values[:, 0:7], axis=1) # SAI # can use this to plot in matlab graspeventdetection_plot.m np.savetxt(path+'/SAI_Fgr.txt', self.Fgr) # can use this to plot in matlab np.savetxt(path+'/SAI_Fgl.txt', self.Fgl) # 0.55*pi rad/samples b1, a1 = signal.butter(1, 0.55, 'high', analog=False) self.f_acc_x = signal.lfilter(b1, a1, self.acc_x, axis=-1, zi=None) self.f_acc_y = signal.lfilter(b1, a1, self.acc_y, axis=-1, zi=None) self.f_acc_z = signal.lfilter(b1, a1, self.acc_z, axis=-1, zi=None) # self.f_eff = signal.lfilter(b1, a1, self.eff, axis=-1, zi=None) # type(eff) self.FAII = np.sqrt(np.square(self.f_acc_x) + np.square(self.f_acc_y) + np.square(self.f_acc_z)) # can use this to plot in matlab np.savetxt(path+'/FAII.txt', self.FAII) # subtract base values from the values array self.values1 = self.values - self.values.min(axis=0) # pass the filter for each sensor self.fvalues1 = np.zeros(self.values1.shape) # 0.48*pi rad/samples b, a = signal.butter(1, 0.48, 'high', analog=False) for i in range(16): self.fvalues1[:, i] = signal.lfilter(b, a, self.values1[:, i], axis=-1, zi=None) self.FAI = np.sum(self.fvalues1, axis=1) # can use this to plot in matlab np.savetxt(path+'/FAI.txt', self.FAI)
def compute_fai(self): filtered_values = signal.lfilter(self.b, self.a, self.sensor_values, axis=0) # print shape() finger1 = filtered_values[-1][0:1].sum() finger2 = filtered_values[-1][1:2].sum() finger3 = filtered_values[-1][2:3].sum() msg = FingerSAI() msg.finger1 = float(finger1) msg.finger2 = float(finger2) msg.finger3 = float(finger3) self.fai_pub.publish(msg) if self.recording: t = time() self.data['fair'].append((t, right)) self.data['fail'].append((t, left)) self.data['faim'].append((t, middle))
def compute_faii(self): filtered_acc = signal.lfilter(self.b1, self.a1, self.acc, axis=0) self.faii = np.sqrt((filtered_acc**2).sum(axis=1)) # Publish just the last value val = self.faii[-1] self.faii_pub.publish(Float64(val)) if self.recording: t = time() self.data['faii'].append((t, val))
def __gausscoeff(s): """Python implementation of the algorithm by Young & Vliet, 1995.""" if s < .5: raise ValueError('Sigma for Gaussian filter must be >0.5 samples') q = 0.98711*s - 0.96330 if s > 0.5 else 3.97156 \ - 4.14554*np.sqrt(1.0 - 0.26891*s) b = np.zeros(4) b[0] = 1.57825 + 2.44413*q + 1.4281*q**2 + 0.422205*q**3 b[1] = 2.44413*q + 2.85619*q**2 + 1.26661*q**3 b[2] = -(1.4281*q**2 + 1.26661*q**3) b[3] = 0.422205*q**3 B = 1.0 - ((b[1] + b[2] + b[3])/b[0]) # convert to a format compatible with lfilter's # difference equation B = np.array([B]) A = np.ones(4) A[1:] = -b[1:]/b[0] return B, A
def Gaussian1D(signal, sigma, padding=0): n = signal.shape[0] tmp = np.zeros(n + padding) if tmp.shape[0] < 4: raise ValueError('Signal and padding too short') tmp[:n] = signal B, A = __gausscoeff(sigma) tmp = lfilter(B, A, tmp) tmp = tmp[::-1] tmp = lfilter(B, A, tmp) tmp = tmp[::-1] return tmp[:n]
def G_inv_mat(x,mode,NT,gs,gd_vec,bas_flag = True, c1_flag = True): """ Fast computation of G^{-1}*x and G^{-T}*x """ from scipy.signal import lfilter if mode == 1: b = lfilter(np.array([1]),np.concatenate([np.array([1.]),-gs]),x[:NT]) + bas_flag*x[NT-1+bas_flag] + c1_flag*gd_vec*x[-1] #b = np.dot(Emat,b) elif mode == 2: #x = np.dot(Emat.T,x) b = np.hstack((np.flipud(lfilter(np.array([1]),np.concatenate([np.array([1.]),-gs]),np.flipud(x))),np.ones(bas_flag)*np.sum(x),np.ones(c1_flag)*np.sum(gd_vec*x))) return b
