On Wed, May 6, 2009 at 10:21 AM, <josef.pktd@gmail.com> wrote:
On Wed, May 6, 2009 at 10:00 AM, Talbot, Gerry <Gerry.Talbot@amd.com> wrote:
Sorry, I guess I wasn't clear, I meant:
for n in xrange(1,N): y[n] = A*x[n] + B*y[n-1]
So y[n-1] is the result from the previous loop iteration.
I was using scipy.signal for this but I have to look up what I did exactly. I think either signal.correlate or using signal.lti.
No, its signal.lfilter, below is a part of a script I used to simulate and estimate an AR(1) process, which is similar to your example. I haven't looked at it in a while but it might give you the general idea. Josef # Simulate AR(1) #-------------- # ar * y = ma * eta ar = [1, -0.8] ma = [1.0] # generate AR data eta = 0.1 * np.random.randn(1000) yar1 = signal.lfilter(ar, ma, eta) etahat = signal.lfilter(ma, ar, y) np.all(etahat == eta) # find error for given filter on data print 'AR(2)' for rho in [0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.79, 0.8, 0.81, 0.9]: etahatr = signal.lfilter(ma, [1, --rho], yar1) print rho,np.sum(etahatr*etahatr) print 'AR(2)' for rho2 in np.linspace(-0.4,0.4,9): etahatr = signal.lfilter(ma, [1, -0.8, -rho2], yar1) print rho2,np.sum(etahatr*etahatr) def errfn(rho): etahatr = signal.lfilter(ma, [1, -rho], yar1) #print rho,np.sum(etahatr*etahatr) return etahatr def errssfn(rho): etahatr = signal.lfilter(ma, [1, -rho], yar1) return np.sum(etahatr*etahatr) resultls = optimize.leastsq(errfn,[0.5]) print 'LS ARMA(1,0)', resultls resultfmin = optimize.fmin(errssfn, 0.5) print 'fminLS ARMA(1,0)', resultfmin