[SciPy-User] python and filter design: calculating optimal "S" transform

[josef.pktd at gmail.com]

On Tue, Jun 1, 2010 at 9:46 PM,  <josef.pktd at gmail.com> wrote:
> On Tue, Jun 1, 2010 at 6:55 PM, robert somerville
> <rbrt.somerville at gmail.com> wrote:
>> Hi;
>> this is an airy question.
>>
>> does anybody have some code or ideas on how to calculate the optimal "S"
>> transform of user specified order (wanting the coefficients)  for a
>> published filter response curve, ie.
>>
>> f(s) = (b1*S^2 + b2*S) / (a1*S^2 + a2*S + a3)
>>
>> I am parameterizing the response of a linear device (Out = Response*In). I
>> have the measured frequency response for the device (amplitude, phase) for a
>> range of frequencies.
>>
>> I wish to model that measured response via a ratio of polynomials in the
>> s-domain (or Laplace domain), where I define the polynomial orders (for
>> numerator and denominator).
>>
>> Something like the "Yule-Walker" method is what I'm after except, to my
>> knowledge, Yule-Walker approach is strictly for responses involving a
>> denominator polynomial (i.e. strictly autoregressive) only. I need some
>> thing to discover both numerator and denominator coefficients.
>>
>
> I have quite a bit of difficulty translating the terminology and to
> understand how "optimal" would be defined in your case.
>
> I have coded up quite a bit for working with ARMA models, especially
> for going in between ma and ar representation (infinite denominator or
> infinite numerator). If I understand the question correctly, then the
> closest I have is to use numerical optimization to solve for the
> ARMA(p,q) for fixed p,q that mimimizes the squared integrated distance
> between the theoretical impulse response function of this process and
> one that is given as target. I don't remember if I ever added a
> function to convert the ma term to an invertible form, since the
> optimization is unconstraint.
>
> If nobody else a more signal theoretic solution, I can look up my code
> which is somewhere in scikits.statsmodels.sandbox.tsa
>
> Josef
>

On Wed, Jun 2, 2010 at 11:28 AM, robert somerville
<rbrt.somerville at gmail.com> wrote:
> My electrical guy says this looks very interesting, and he would like to see
> if it is what we are are looking for, although he says it looks very close .
>
> we are trying to model the impulse response of some geophysical
> instrumentation from published response curves. I believe we are trying to
> determine the best coefficients in the modeling polynomial
>
> Thanks;
> Robert Somerville
>

(Can you try to reply to the thread to keep the information together?)

some documentation is here:
http://statsmodels.sourceforge.net/sandbox.html#time-series-analysis-tsa

for most parts I have both sample and theoretical properties of arma processes

for example
theoretical impulse response function of ARMA process:

http://bazaar.launchpad.net/~scipystats/statsmodels/trunk/annotate/head%3A/scikits/statsmodels/sandbox/tsa/arima.py#L278

I think the conversion code and the theoretical acf, pac, impulse
response functions are fully tested,
What I have written in frequency domain is more eclectic and less
tested, because I have less experience working in frequency than time
domain. I think theoretical spectrum is easy, but I don't remember how
far I got with theoretical impulse response function in frequency
domain.

The closest I have to what you might need,  is ar2arma

http://bazaar.launchpad.net/~scipystats/statsmodels/trunk/annotate/head%3A/scikits/statsmodels/sandbox/tsa/try_fi.py#L79

def ar2arma(ar_des, p, q, n=20, mse='ar', start=None):
    '''find arma approximation to ar process

    This finds the ARMA(p,q) coefficients that minimize the integrated
    squared difference between the impulse_response functions
    (MA representation) of the AR and the ARMA process. This does
    currently not check whether the MA lagpolynomial of the ARMA
    process is invertible, neither does it check the roots of the AR
    lagpolynomial.

This was written to approximate an infinite order AR process
(fractionally integrated ARMA), but I think it should be possible to
change the function to match a given impulse response function,
instead of the one implied by ar_des. The function has some examples
but is not really tested, because didn't finish the missing pieces to
get the invertible process.

some functions where written to help me with the translation of the
terminology between signal analysis, scipy.signal and time series
analysis,   essentially ma, ar are the same as num, denom (but I never
remember which is which). The lag polynomials sometimes have the
leading 1 sometimes not. general for is ar(L) y_t = ma(L) u_t , with
u_t input and y_t output, L lag-operator
In the optimization and estimation, ar(0)=ma(0)=1 is usually assumed.

Note: I worked on this heavily more than half a year ago, and have
barely looked at it since.

Hope that helps and I would appreciate any feedback.

Josef


>
>
>
>
>
>
>>
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
>