I've just read what it says about Discrete Hartley Transform in: http://en.wikipedia.org/wiki/Fast_fourier_transform#FFT_algorithms_specialized_for_real_and.2For_symmetric_data But this is such a weaselly assertion, I've a good mind to try it out, when I get round to it.
On Wed, Jun 27, 2012 at 8:11 AM, Ian Clark <earthspo...@gmail.com> wrote: >> Why would a Hartley transform be faster than the Fast Fourier Transform? > > That should be "might" not "would", because it's an off-the-cuff > conjecture, until I've tried it. The Hartley transform works entirely > with real numbers, avoiding the need for symmetrisation, which nearly > doubles the size of the processed vector. It avoids the use of (o.) to > ensure complex datatype doesn't creep in. Also multiplying a complex > vector (by a floating scalar) takes twice as long as multiplying a > same-length real vector ...at least it does in my rough timer tests. > > J wiki says the addon "math/fftw" works only for Windows, so I assumed > FFT was not available in J for me as a Mac user. I am wrong. I see it > uses a dylib (libfftw3.3.dylib) when loaded on my iMac, and it appears > to work. (I have added a note to this effect to: > http://www.jsoftware.com/jwiki/Addons/math/fftw .) > > I'm glad I've discovered that. I last played with FFT way back in 1974 > and was impressed then with what it could do. > > So my earlier thoughts about adapting Cooley-Tukey to the Hartley > transform may be wasted effort since FFTW may outperform it anyway. > Even with symmetrization. > > However, provided it runs in minutes, not hours or days, performance > is low priority for me on this task. And I rather like the fact that > ABC computes C first, because then I have the option of telling it the > "true" C if I know it. > > Ian > > On Tue, Jun 26, 2012 at 8:13 AM, Bo Jacoby <bojac...@yahoo.dk> wrote: >> Hi Ian. >> >> >> I'm glad that your test was successful, albeit slow. The important >> optimization is to replace my slow fourier transform by a fast fourier >> transform. Your example has 100 items in the time series, and the real >> transform RT asks the fourier transform FT to perform 98 vector >> multiplications by a 198*198 complex matrix, (98*198*198)=3841992, which is >> embarassingly inefficient. The FFT will speed it up by a factor >> ((%2&^.)198)=26. >> >> >> Another optimization is to incorporate the computation of A and B into the >> computation of C, rather than starting from scratch after computing C. >> >> The choice of index origin for C is a matter of taste. >> >> >> Why would a Hartley transform be faster than the Fast Fourier Transform? >> >> >> It remains to estimate the inaccuracies of A and B. It is the >> root-mean-square value of the noise vector X-A+B*C<#X. Note that the noise >> has only (#X)-2 degrees of freedom, because A and B have been estimated, so >> the root-mean-square is >> >> >> RMS=:(+/%#-2:)&.:*: >> >> - Bo >> >> >>>________________________________ >>> Fra: Ian Clark <earthspo...@gmail.com> >>>Til: Programming forum <programming@jsoftware.com> >>>Sendt: 2:41 tirsdag den 26. juni 2012 >>>Emne: Re: [Jprogramming] Kalman filter in J? >>> >>>...I should add that my CC returns 1+(the value of your 'C'). This is >>>for compatibility with how I generate the step fn, using (the >>>constant) C. >>> >>>On Tue, Jun 26, 2012 at 1:37 AM, Ian Clark <earthspo...@gmail.com> wrote: >>>> @Bo -Here is your solution running in my test rig... >>>> >>>> A= 100 bias constant (popn mean of X before the step) >>>> B= 20 scaling factor for height of step >>>> C= 90 epoch of actual occurrence of the Heaviside step >>>> D= 10 scaling factor for Gaussian noise component (sigma) >>>> N= 100 number of epochs in time series X >>>> step detected at: (wait...) >>>> 90 >>>> ABC X >>>> 99.8246 20.0353 90 >>>> >>>> ABC takes approx 10 s to run on my iMac. >>>> CC (--your 'C') takes approx 5 s. >>>> >>>> I'm actually quite impressed with how accurately it estimates A, B, C, >>>> compared with other methods (which don't do A or B). I'm writing up my >>>> experiments and will post a wiki link. >>>> >>>> Whilst I like your symmetrisation idea, might instead using the >>>> discrete Hartley Transform (which works entirely in the real domain) >>>> speed things up? >>>> >>>> Ian >>>> >>>> On Mon, Jun 25, 2012 at 12:00 PM, Bo Jacoby <bojac...@yahoo.dk> wrote: >>>>> Hi Ian. Here is your solution. >>>>> >>>>> ABC=: 3 : 0 >>>>> c=.C y >>>>> y=.(0,>:c){RT"1&(,H&#)y >>>>> y=.y,(%/1{"1 y)*1{y >>>>> y=.y,(0{y)-2{y >>>>> y=.y,(#|:y){.{.{:y >>>>> y=.4 2{{:"1 RT"1 y >>>>> y,c >>>>> ) >>>>> ABC 101 103 100 210 230 200 >>>>> 101.15 115.101 2 >>>>> >>>>> >>>>> - Bo >>>>> >>>>> >>>>> >>>>> >>>>>>________________________________ >>>>>> Fra: Ian Clark <earthspo...@gmail.com> >>>>>>Til: Programming forum <programming@jsoftware.com> >>>>>>Sendt: 6:46 mandag den 25. juni 2012 >>>>>>Emne: Re: [Jprogramming] Kalman filter in J? >>>>>> >>>>>>Interesting approach, Bo. >>>>>> >>>>>>Turned my attention to CUSUM (an adaptive filter) which is rather easy >>>>>>to compute and gives good results, but needs parameters setting "by >>>>>>eye". I've yet to try a Bayesian approach, which I suspect will be the >>>>>>most sensitive of all. >>>>>> >>>>>>But I did play with your FT, using *:@(10&o.) to plot a power >>>>>>spectrum. I'm using rather noisier data than your example: my step is >>>>>>only roughly (sigma) high. Nevertheless with the onset of the step I >>>>>>can clearly see a high-frequency spike appear at the far end of the >>>>>>transformed time-series, due to the transformed Heaviside fn. Same >>>>>>problem as with the CUSUM statistic however: designing a detector >>>>>>which can be adjusted for false negatives and positives, and doesn't >>>>>>take too many subsequent samples to detect the step. >>>>>> >>>>>>Will try out your C function in my context, and try to tweak it. Your >>>>>>solution detects when the step happens, solving 1 in my stated >>>>>>objectives: 1-4. I fear however that I am more interested in 2-4. >>>>>> >>>>>> >>>>>> >>>>>>On Sun, Jun 24, 2012 at 6:06 PM, Bo Jacoby <bojac...@yahoo.dk> wrote: >>>>>>> Hi Ian >>>>>>> The Fourier Transform FT transforms a complex n-vector X into a complex >>>>>>> n-vector Y.If X is real, then Y is symmetric. And if X is symmetric >>>>>>> then Y is real. So I think it is a good idea to symmetrize X before >>>>>>> taking the FT. (This makes a slow program even slower. Optimize later!) >>>>>>> >>>>>>> >>>>>>> NB.sm=: symmetrize >>>>>>> sm=:[:,}:"1&(,:|.) >>>>>>> sm i.6 NB. see what I mean >>>>>>> >>>>>>> 0 1 2 3 4 5 4 3 2 1 >>>>>>> NB.RT=: Real Transform. (9 o. removes every trace of complexity) >>>>>>> RT=:9 o.#{.FT&sm >>>>>>> NB.H=:Heaviside steps >>>>>>> H=:i.&<:</i. >>>>>>> >>>>>>> NB.C=:Index of the last item before the change >>>>>>> C=:[:(i.<./)[:+/"1[:*:[:}.[:(-"1{.)[:}."1[:(%{."1)[:}."1[:RT"1],[:H# >>>>>>> NB.test >>>>>>> C 101 204 202 200 203 >>>>>>> 0 >>>>>>> C 101 104 202 200 203 >>>>>>> 1 >>>>>>> C 101 104 102 200 203 >>>>>>> 2 >>>>>>> C 101 104 102 100 203 >>>>>>> 3 >>>>>>> >>>>>>> This is not a Kálmán filter, but it solves the problem. >>>>>>> - Bo >>>>>>> >>>>>>> >>>>>>> >>>>>>>>________________________________ >>>>>>>> Fra: Ian Clark <earthspo...@gmail.com> >>>>>>>>Til: Programming forum <programming@jsoftware.com> >>>>>>>>Sendt: 3:20 fredag den 22. juni 2012 >>>>>>>>Emne: Re: [Jprogramming] Kalman filter in J? >>>>>>>> >>>>>>>>Thanks, Bo. I'll look at it when I get back, probably over the weekend. >>>>>>>> >>>>>>>>Ian >>>>>>>> >>>>>>>>On Fri, Jun 22, 2012 at 12:44 AM, Bo Jacoby <bojac...@yahoo.dk> wrote: >>>>>>>>> Hi Ian. >>>>>>>>> The following is a tool, but not yet a solution. >>>>>>>>> >>>>>>>>> FT is a slow Fourier Transform. When the problem has been solved it >>>>>>>>> is time to optimize. This FT has the nice property that FT^:_1 is >>>>>>>>> identical to FT. >>>>>>>>> >>>>>>>>> NB.FT =: Fourier Transformation >>>>>>>>> FT =: +/&(+*((%:%~_1&^&+:&(%~(*/])&i.))&#)) >>>>>>>>> NB.rd =: round complex number. >>>>>>>>> rd =: (**|)&.+. >>>>>>>>> NB.H =: Heaviside step function >>>>>>>>> H =: (]#0:),-#1: >>>>>>>>> 10 H 3 >>>>>>>>> 0 0 0 1 1 1 1 1 1 1 >>>>>>>>> NB. FT = FT^:_1 >>>>>>>>> rd FT FT 10 H 3 >>>>>>>>> 0 0 0 1 1 1 1 1 1 1 >>>>>>>>> >>>>>>>>> >>>>>>>>> NB. For fixed C the values of A and B to minimize +/*:X-A+B*(#X)H C >>>>>>>>> are found by solving (2{.FT X) = 2{.FT A+B*(#X)H C >>>>>>>>> >>>>>>>>> - Bo >>>>>>>>>>________________________________ >>>>>>>>>> Fra: Ian Clark <earthspo...@gmail.com> >>>>>>>>>>Til: Programming forum <programming@jsoftware.com> >>>>>>>>>>Sendt: 18:15 torsdag den 21. juni 2012 >>>>>>>>>>Emne: Re: [Jprogramming] Kalman filter in J? >>>>>>>>>> >>>>>>>>>>Yes, Bo, that looks like a fair statement of the problem, in terms of >>>>>>>>>>least-squares minimization. >>>>>>>>>> >>>>>>>>>> >>>>>>>>>>On Thu, Jun 21, 2012 at 4:46 PM, Bo Jacoby <bojac...@yahoo.dk> wrote: >>>>>>>>>>> Consider the step functions >>>>>>>>>>> >>>>>>>>>>> H=:(]#0:),-#1: >>>>>>>>>>> 10 H 0 >>>>>>>>>>> 1 1 1 1 1 1 1 1 1 1 >>>>>>>>>>> 10 H 3 >>>>>>>>>>> 0 0 0 1 1 1 1 1 1 1 >>>>>>>>>>> 10 H 10 >>>>>>>>>>> 0 0 0 0 0 0 0 0 0 0 >>>>>>>>>>> The problem is to find constants A,B,C to minimize the expression >>>>>>>>>>> >>>>>>>>>>> +/*:X-A+B*(#X)H C >>>>>>>>>>> Right? >>>>>>>>>>> -Bo >>>>>>>>>>>>________________________________ >>>>>>>>>>>> Fra: Ian Clark <earthspo...@gmail.com> >>>>>>>>>>>>Til: Programming forum <programming@jsoftware.com> >>>>>>>>>>>>Sendt: 14:34 torsdag den 21. juni 2012 >>>>>>>>>>>>Emne: Re: [Jprogramming] Kalman filter in J? >>>>>>>>>>>> >>>>>>>>>>>>Thanks, Bo and Ric. Yes, I'd tried searching jwiki on "Kalman" and >>>>>>>>>>>>Pieter's page was the only instance. >>>>>>>>>>>> >>>>>>>>>>>>I neglected to mention that in the most general case I can't really >>>>>>>>>>>>be >>>>>>>>>>>>confident that X1 and X2 have the same distribution function, let >>>>>>>>>>>>alone the same variance. But looking at it again, I see that under >>>>>>>>>>>>the >>>>>>>>>>>>restrictions I've placed the problem simplifies immensely to >>>>>>>>>>>>fitting a >>>>>>>>>>>>step-function H to: X=: K+G+H. If I can just do that, I'll be happy >>>>>>>>>>>>for now. >>>>>>>>>>>> >>>>>>>>>>>>Repeated application of FFT should allow me to subtract the noise >>>>>>>>>>>>spectrum F(G), or at least see a significant change in the overall >>>>>>>>>>>>spectrum emerge after point T, and that might handle the more >>>>>>>>>>>>general >>>>>>>>>>>>cases as well. >>>>>>>>>>>> >>>>>>>>>>>>Anyway it's simple-minded enough for me, and worth a try. FFT is, >>>>>>>>>>>>after all, "fast" :-) >>>>>>>>>>>> >>>>>>>>>>>>But won't an even faster transform do the trick, such as (+/X)? On >>>>>>>>>>>>the >>>>>>>>>>>>above model, X performs a drunkard's walk around a value M1 until >>>>>>>>>>>>some >>>>>>>>>>>>point T, after which it walks around M2. Solution: simply estimate >>>>>>>>>>>>M1 >>>>>>>>>>>>and M2 on an ongoing basis. >>>>>>>>>>>> >>>>>>>>>>>>I get the feeling I ought to be searching on terms like "edge >>>>>>>>>>>>detection", "step detection" and CUSUM. >>>>>>>>>>>> >>>>>>>>>>>>Anyway, there's enough here to try. >>>>>>>>>>>> >>>>>>>>>>>>On Thu, Jun 21, 2012 at 10:28 AM, Ric Sherlock <tikk...@gmail.com> >>>>>>>>>>>>wrote: >>>>>>>>>>>>> Hi Ian, >>>>>>>>>>>>> >>>>>>>>>>>>> A quick search of the J wiki finds this: >>>>>>>>>>>>> http://www.jsoftware.com/jwiki/Stories/PietdeJong >>>>>>>>>>>>> Sounds like he might have what you're after? >>>>>>>>>>>>> >>>>>>>>>>>>> Cheers, >>>>>>>>>>>>> Ric >>>>>>>>>>>>> >>>>>>>>>>>>> On Thu, Jun 21, 2012 at 4:33 PM, Ian Clark >>>>>>>>>>>>> <earthspo...@gmail.com> wrote: >>>>>>>>>>>>>> Can anyone help? Has anyone written a Kalman filter in J? >>>>>>>>>>>>>> >>>>>>>>>>>>>> I'm not a specialist in either statistics or control theory, so >>>>>>>>>>>>>> I'm >>>>>>>>>>>>>> only guessing a Kalman filter is what I need. Though I do have a >>>>>>>>>>>>>> passing acquaintance with the terms: stochastic control and >>>>>>>>>>>>>> linear >>>>>>>>>>>>>> quadratic Gaussian (LQG) control. I am aware that a "Kalman >>>>>>>>>>>>>> filter" >>>>>>>>>>>>>> (like ANOVA) is more a topic than a black-box. >>>>>>>>>>>>>> >>>>>>>>>>>>>> So let me explain what I want it for. >>>>>>>>>>>>>> >>>>>>>>>>>>>> I have a time series X which I am assuming can be modelled like >>>>>>>>>>>>>> this: >>>>>>>>>>>>>> >>>>>>>>>>>>>> X=: K + G + (X1,X2) >>>>>>>>>>>>>> >>>>>>>>>>>>>> where >>>>>>>>>>>>>> >>>>>>>>>>>>>> K is constant >>>>>>>>>>>>>> G is Gaussian noise >>>>>>>>>>>>>> X1 is a random variable with mean: M1 and variance: V1 >>>>>>>>>>>>>> X2 is a random variable with mean: M2 and variance: V2 >>>>>>>>>>>>>> >>>>>>>>>>>>>> Typically X is a sequence of sensor readings, but may also be >>>>>>>>>>>>>> measurements from a series of user trials conducted on a working >>>>>>>>>>>>>> prototype, which suffers a design-change at a given point T. >>>>>>>>>>>>>> >>>>>>>>>>>>>> Simplifying assumptions (which unfortunately I may need to relax >>>>>>>>>>>>>> in due course): >>>>>>>>>>>>>> >>>>>>>>>>>>>> (a) X is not multivariate >>>>>>>>>>>>>> (b) X1 and X2 are Gaussian >>>>>>>>>>>>>> (c) V1=V2 (only the mean value changes, not the variance). >>>>>>>>>>>>>> >>>>>>>>>>>>>> The problem: >>>>>>>>>>>>>> >>>>>>>>>>>>>> 1. Estimate T=: 1+#X1 -- the point at which X1 gives way to X2. >>>>>>>>>>>>>> >>>>>>>>>>>>>> 2. Given T, estimate (M2-M1) -- the "underlying improvement", if >>>>>>>>>>>>>> any, >>>>>>>>>>>>>> of the change to the prototype. >>>>>>>>>>>>>> >>>>>>>>>>>>>> 3. (subcase of 2.) Given T, test the null hypothesis: M1=M2, viz >>>>>>>>>>>>>> that >>>>>>>>>>>>>> there has been no underlying improvement. >>>>>>>>>>>>>> >>>>>>>>>>>>>> 4. Estimate U=: #X2 -- the minimum number of samples needed >>>>>>>>>>>>>> after T in >>>>>>>>>>>>>> order to achieve 1-3 above with 95% confidence. >>>>>>>>>>>>>> >>>>>>>>>>>>>> In other words, detect the signal-in-noise: M1-->M2, and do so >>>>>>>>>>>>>> in real-time. >>>>>>>>>>>>>> >>>>>>>>>>>>>> Because of 4, the need to estimate T and (M2-M1) on an ongoing >>>>>>>>>>>>>> basis, >>>>>>>>>>>>>> I can't do a randomised block design. I gather that a Kalman >>>>>>>>>>>>>> filter, >>>>>>>>>>>>>> or some sort of adaptive filter, will handle this problem. >>>>>>>>>>>>>> >>>>>>>>>>>>>> But maybe something simpler will turn out good enough? >>>>>>>>>>>>>> >>>>>>>>>>>>>> Supposing I can get hold of a "black box" Kalman filter, I >>>>>>>>>>>>>> propose to >>>>>>>>>>>>>> test it out on generated data and compare its performance to some >>>>>>>>>>>>>> simple-minded approach, like estimating M1 / M2 from a simple >>>>>>>>>>>>>> moving >>>>>>>>>>>>>> average of the last U samples, or applying the F-test to 2 sets >>>>>>>>>>>>>> of U >>>>>>>>>>>>>> samples taken either side of T. >>>>>>>>>>>>>> >>>>>>>>>>>>>> But since the technique aims to be published, or at least >>>>>>>>>>>>>> critically >>>>>>>>>>>>>> scrutinised (and maybe incorporated in a software product), I'd >>>>>>>>>>>>>> rather >>>>>>>>>>>>>> depend on a state-of-art packaged solution than reinvent the >>>>>>>>>>>>>> wheel: a >>>>>>>>>>>>>> large and very well-turned wheel it appears to me. >>>>>>>>>>>>>> >>>>>>>>>>>>>> Ian Clark >>>>>>>>>>>>>> ---------------------------------------------------------------------- >>>>>>>>>>>>>> For information about J forums see >>>>>>>>>>>>>> http://www.jsoftware.com/forums.htm >>>>>>>>>>>>---------------------------------------------------------------------- >>>>>>>>>>>>For information about J forums see >>>>>>>>>>>>http://www.jsoftware.com/forums.htm >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>> ---------------------------------------------------------------------- >>>>>>>>>>> For information about J forums see >>>>>>>>>>> http://www.jsoftware.com/forums.htm >>>>>>>>>>---------------------------------------------------------------------- >>>>>>>>>>For information about J forums see http://www.jsoftware.com/forums.htm >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>> ---------------------------------------------------------------------- >>>>>>>>> For information about J forums see http://www.jsoftware.com/forums.htm >>>>>>>>---------------------------------------------------------------------- >>>>>>>>For information about J forums see http://www.jsoftware.com/forums.htm >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>> ---------------------------------------------------------------------- >>>>>>> For information about J forums see http://www.jsoftware.com/forums.htm >>>>>>---------------------------------------------------------------------- >>>>>>For information about J forums see http://www.jsoftware.com/forums.htm >>>>>> >>>>>> >>>>>> >>>>> ---------------------------------------------------------------------- >>>>> For information about J forums see http://www.jsoftware.com/forums.htm >>>---------------------------------------------------------------------- >>>For information about J forums see http://www.jsoftware.com/forums.htm >>> >>> >>> >> ---------------------------------------------------------------------- >> For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
