thx for help. the problem was a state matrix (A). after transposition it's
ok (stupid error however difficult to find ... first good predictions made
me think that matrices were ok)
Reply by Michael Wirtzfeld●June 10, 20102010-06-10
grzegorz (Original Poster),
With respect to your Python code, I suggest you have a look at the line, "IS = R + dot(H, dot(P, H.T))", in your "update" function definition. This goes back to a previous idea of this thread.
As suggested by Michael (Plante), taking time to "Describe the problem (inputs/outputs/model) of the filter..." would help us understand the framework of your problem.
Your code fails, I believe, because the "IS" matrix becomes ill-conditioned as the observations are processed. From what I can tell it stops at the 5th iteration. I think the root-problem lies in the respective sizes of the principle process/model vectors and matrices you are using.
The "dot(H, dot(P, H.T))" portion of the above equation determines the innovation or residual covariance update. Looking at the sizes of the measurement/observation model vector (H) and the predicted (a priori) covariance estimate (P), might indicate the potential underlying issue.
"H", as you've defined in your code, is 1-by-3 and "P" is 3-by-3. With reference to, "http://tinyurl.com/Wiki-Kalman-Reference", the size of the innovation covariance update is then:
1-by-3 times 3-by-3 times 3-by-1 yielding a 1-by-1 or scalar result. Should it not be a 3-by-3 matrix result in light of the 3-by-3 measurement covariance matrix, "R"?
I am not 100% sure how Python handles the dimensions of the "array" and "mat" (matrix) types. In the event the measurement model vector's dimensions are reversed, then you have:
3-by-1 times 3-by-3 times 1-by-3 yielding a no-go invalid matrix result.
With the scalar result, the same value is added to all elements of the measurement covariance matrix, "R". With continued processing, your innovation covariance grows resulting in an ill-conditioned matrix "IS" that can not be inverted.
Perhaps you can flesh-out the vectors and matrices involved with your process and measurement models as suggested by Michael (Plante). From what I can see, you have:
----------
State Transition Matrix (A):
[[ 1. 1. 0.5]
[ 0. 1. 1. ]
[ 0. 0. 1. ]]
Process Noise Covariance Matrix (Q):
[[ 1.0 0. 0. ]
[ 0. 1.0 0. ]
[ 0. 0. 1.0]]
Measurement/Observation Model (H):
[1 0 0]
Measurement Noise Covariance Model (R):
[[ 1.0 0. 0. ]
[ 0. 1.0 0. ]
[ 0. 0. 1.0]]
Initial Covariance Estimate:
[[ 0.01 0. 0.]
[ 0. 0.01 0.]
[ 0. 0. 0.01]]
Initial State Estimate:
[0 0 0]
----------
Application details would be helpful... There are sections of your code that you may also what to clarify. For example:
# temporary freq. and modulation factor - not important if H= [1 0 0]
for n in range(1,n_iter):# dla n = 0 -> omega,alfa=0
omega_list[n]=(teta_list[n+1]-teta_list[n-1])/2
alfa_list[n]=teta_list[n+1]-2*teta_list[n]+teta_list[n-1]
for n in range(n_iter):
X[n][0]=teta_list[n]
# X[n][1]=omega_list[n]
# X[n][2]=alfa_list[n]
This is used to generate the vector "X" later on in your code to be used as your measurement or observations vector.
I hope all or part of this helps...
Michael.
---
frmsrcurl: http://compgroups.net/comp.dsp/kalman-filter-python-implementation
Reply by Tim Wescott●June 9, 20102010-06-09
On 06/09/2010 03:34 PM, Erik de Castro Lopo wrote:
> Tim Wescott wrote:
>
>> Matlab is _expensive_.
>
> Octave is a cheap alternative.
>
>> While it delivers considerable value, there are
>> free alternatives that do as well or better than Matlab. From a
>> standpoint of how much you gain vs. how much you pay, Matlab is _horrible_.
>
> I agree.
>
I started using Scilab when the Octave project was in a shambles. Now
the Scilab project is behind (still no debugger...) but I'm used to it.
--
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com
Reply by Erik de Castro Lopo●June 9, 20102010-06-09
Tim Wescott wrote:
> Matlab is _expensive_.
Octave is a cheap alternative.
> While it delivers considerable value, there are
> free alternatives that do as well or better than Matlab. From a
> standpoint of how much you gain vs. how much you pay, Matlab is _horrible_.
I agree.
For prototyping in a high level language I usually use Ocaml or
Haskell, but I would never expect anyone not familiar with those
languages to look at my code.
Erik
--
----------------------------------------------------------------------
Erik de Castro Lopo
http://www.mega-nerd.com/
Reply by Jerry Avins●June 9, 20102010-06-09
On 6/9/2010 11:24 AM, Michael wrote:
> Jerry,
>
> I did not insult the poster or the fact the he is using Python as his programming language of choice - "brickbats" is harsh to say the least. I am, however, keenly aware of the tone of his last message.
>
> It is okay to be frustrated with a problem you are trying to solve, but extending it into a public forum where you ask individuals for their time and their experiences to help is another.
>
>
> Michael.
> --
> Finding Nemo, Sharkbait - "He's probably American!"
>
> ---
> frmsrcurl: http://compgroups.net/comp.dsp/kalman-filter-python-implementation
Michael,
Yours was neither the first nor the more offensive letter critical of
his choice of language. Assume that his annoyance arose from the
accumulation, not from your response alone.
We tend to observe certain standards here, communicating in English and
Matlab chief among them. Nevertheless, we all should recognize that this
is an international forum and that there are many legitimate programming
languages.
All that said, nobody here is under any obligation to read a long
program or routine. When a question begins with a few pages of code in
any language, I usually ignore it.
Jerry
--
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������
Reply by Tim Wescott●June 9, 20102010-06-09
On 06/09/2010 12:07 PM, HardySpicer wrote:
> On Jun 9, 3:18 pm, "Nasser M. Abbasi"<n...@12000.org> wrote:
>> On 6/8/2010 6:55 PM, HardySpicer wrote:
>>
>>> On Jun 9, 7:41 am, "grzegorz g."
>>
>>> Why the Hell Python? Why not use Javascript while you're at it!
>>
>> fyi, Python is used a lot in mathematics software these days. Sage uses
>> Python, Numpy is in Python, as well as SciPy, and many more. Many
>> people write scientific software in Python.
>>
>> --Nasser
>
> What the hell is Sage and Numpty? Ever heard of Matlab? Far easier.
Until you go to distribute it, and find out that you have to pay enough
to buy a really nice car for each copy of Matlab that you give away.
Matlab is _expensive_. While it delivers considerable value, there are
free alternatives that do as well or better than Matlab. From a
standpoint of how much you gain vs. how much you pay, Matlab is _horrible_.
--
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com
Reply by grzegorz g.●June 9, 20102010-06-09
>Jerry,
>
>I did not insult the poster or the fact the he is using Python as his
programming language of choice - "brickbats" is harsh to say the least. I
am, however, keenly aware of the tone of his last message.
>
>It is okay to be frustrated with a problem you are trying to solve, but
extending it into a public forum where you ask individuals for their time
and their experiences to help is another.
>
>
>Michael.
Michael, I also didn't insult anybody and I didn't extend a frustration. I
only reacted on offtop (often agressive, without any arguments)
>[cut]
>> Can someone help me, cos I have no idea what's wrong ?
>
>
>Mmm... a bit too much code there :)
>Why not you try to the same computation using... err... octave?... and
>check numbers step-by-step to understand where the error is?
>
>
>I do python, so I could help... do you want to share directly .py
>files?
>
On Jun 9, 3:18�pm, "Nasser M. Abbasi" <n...@12000.org> wrote:
> On 6/8/2010 6:55 PM, HardySpicer wrote:
>
> > On Jun 9, 7:41 am, "grzegorz g."
>
> > Why the Hell Python? Why not use Javascript while you're at it!
>
> fyi, Python is used a lot in mathematics software these days. Sage uses
> Python, Numpy is in Python, as well as SciPy, and many more. �Many
> people write scientific software in Python.
>
> --Nasser
What the hell is Sage and Numpty? Ever heard of Matlab? Far easier.
Reply by claudegps●June 9, 20102010-06-09
[cut]
> Can someone help me, cos I have no idea what's wrong ?
Mmm... a bit too much code there :)
Why not you try to the same computation using... err... octave?... and
check numbers step-by-step to understand where the error is?
I do python, so I could help... do you want to share directly .py
files?
Reply by Tim Wescott●June 9, 20102010-06-09
On 06/09/2010 09:35 AM, Rob Gaddi wrote:
> On 6/8/2010 9:05 PM, Erik de Castro Lopo wrote:
>> Nasser M. Abbasi wrote:
>>
>>> On 6/8/2010 6:55 PM, HardySpicer wrote:
>>>
>>>> On Jun 9, 7:41 am, "grzegorz g."
>>>>
>>>> Why the Hell Python? Why not use Javascript while you're at it!
>>>
>>> fyi, Python is used a lot in mathematics software these days.
>>
>> But its not used much for DSP code. For DSP its usually C, C++
>> or matlab/octave.
>>
>>> Sage uses
>>> Python, Numpy is in Python, as well as SciPy, and many more. Many
>>> people write scientific software in Python.
>>
>> Many people drink too much alcohol and then drive motor vehicles.
>> That still doesn't make it a good idea.
>>
>> Erik
>
> There's no reason not to. It's handy in that it really allows you to
> bang out algorithms quickly rather than getting caught up in
> implementational details. Then if the ultimate target is some
> overpowered Core Duo, you might as well leave it Python (the speed hit's
> not as bad as you think), or if it's not then you rewrite it in C, C++,
> Ada, Forth, assembly, VHDL, etc depending on your target, but at least
> you know that if it doesn't work the error's in the porting, not the
> algorithm.
>
Which is exactly what you'd say about Scilab, Octave, Matlab, etc.
--
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com