Venturing further afield, it looks like normalizing means instead of using
kp =: *&$ ($,) 0 2 1 3 |: */
use
Kp =: (%%:2) * *&$ ($,) 0 2 1 3 |: */
Or, put differently, multiply the unnormalized result r by
2^(1-2^.#r)%2 or (same numbers) by -@-:@<:&.(2&^.) #r or (again, same
numbers) 0.5 _0.5&p.&.(2&^.) #r
Or, in other words, values from the sequence %:1 2 4 8 16 32 ... with
1 being used for 2 by 2 arrays. This suggests that the order 0 Harr
matrix should be ,.%:0.5
People mostly ignore this, it seems, because it slightly obscures the
rest of the structure, and because matrix inverse works just fine with
or without it.
Ordering... I'm going to guess that Mike Day's H2N represents the
canonical order for the Haar wavelet. Mostly, because I can't think of
any reason to pick any other order and because the formulation "feels
right".
In other words:
HaarM=: (Kp&(,:1 1), Kp&(,:1 _1)@(=@i.@#))@]^:[&(,.%:0.5)
uHaarM=: (kp&(,:1 1), kp&(,:1 _1)@(=@i.@#))@]^:[&(,.1)
HaarM 2
0.353553 0.353553 0.353553 0.353553
0.5 0.5 _0.5 _0.5
0.707107 _0.707107 0 0
0 0 0.707107 _0.707107
uHaarM 2
1 1 1 1
1 1 _1 _1
1 _1 0 0
0 0 1 _1
Wrapped up as transforms:
uhaar=: +/ .*~ uHaarM@(2 ^. #)
iuhaar=: +/ .*~ %.@uHaarM@(2 ^. #)
And the normalized versions follow the same pattern (except no 'u' in
the names nor definitions)
Revisiting my (not normalized) original definitions, to get things in
the proper order, they should have been:
Haar=: _2&(Haar^:(1<#)@(+/\), -/\)
and
iHaar -:@(+/,@,.-/)@(,:~iHaar^:(1<#))/@($~ 2,-:@#)
I need to think a bit more about whether these versions should be
normalized - normalizing makes them significantly more complex:
(haar % uhaar) p: i.2
0.5 0.707107
(haar % uhaar) p: i.4
0.353553 0.5 0.707107 0.707107
(haar % uhaar) p: i.8
0.25 0.353553 0.5 0.5 0.707107 0.707107 0.707107 0.707107
Anyways, I've probably gone way past the saturation point on this
particular wavelet, and I should probably be looking at other
wavelets, instead.
Thanks,
--
Raul
On Mon, Jan 8, 2018 at 10:57 AM, Raul Miller <[email protected]> wrote:
> Rather than focus on performance (fun, yes, but probably premature),
> I'd like to think a bit about correctness.
>
> Trying a random example from that wikipedia page, I notice two things:
>
> [1] Different ordering (which you also mentioned, briefly), and
>
> [2] Scaling seems to be an issue (the wikipedia article suggests
> "normalizing" the transform).
>
> But I do not know what criteria are relevant in judging correctness of
> an approach for either of those issues, and I've yet to spot anything
> in that wikipedia writeup which seems to help here.
>
> Thoughts?
>
> Thanks,
>
> --
> Raul
>
>
> On Mon, Jan 8, 2018 at 5:01 AM, 'Mike Day' via Programming
> <[email protected]> wrote:
>> However, on this Lenovo Windows 10 laptop,
>>
>> also in a brief trial for q =: 128?10000, I find
>>
>> Haar < hw ~= hw3 << applyHaar for time,
>> and
>>
>> Haar < hw < hw3 << applyHaar for space - as on the iPad
>>
>> (details below)
>>
>> fwiw, I realised that the definition for the H2n matrix
>> required for Raul's ordering is similar in spirit to that
>> in the Wiki article:
>>
>> H2 = [1 _1]
>> [1 1]
>>
>> H2n = [In X [1 _1]]
>> [Hn X [1 1]]
>>
>> where X is the (matrix) Kronecker Product.
>> So now I have:
>>
>> H2n =: 3 : 0 NB. H2n as redefined above for Raul's ordering
>>
>> n =. y
>>
>> h =. 1 _1,:1 1
>>
>> n2=. 2
>>
>> for_i. i. <: 2 <.@^. n do.
>>
>> h =. (h kp ,:1 1) ,~ (=i. n2) kp ,: 1 _1
>>
>> n2=. +:n2
>>
>> end.
>>
>> h
>>
>> )
>>
>> It is fairly easy to define iH2n, the inverse of H2n to
>>
>> allow a matrix-based inverse of applyHaar. I won't bore
>>
>> you with it here.
>>
>> Mike
>>
>> brief timing details:
>> ts'Haar qk'
>>
>> 0.000101511 33920
>>
>> ts'hw qk'
>>
>> 0.000178184 75776
>>
>> ts'hw3 qk'
>>
>> 0.000169545 90880
>>
>> ts'applyHaar qk'
>>
>> 0.104641 3.83121e7
>>
>>
>> On 07/01/2018 23:34, 'Mike Day' via Programming wrote:
>>>
>>> In a brief trial for q =: 128?10000, on this iPad, I find
>>> hw3 < hw < Haar << applyHaar for time,
>>> and
>>> Haar < hw < hw3 << applyHaar for space.
>>> (details below)
>>>
>>> So the matrix method wins no prizes for performance, although examining
>>> the matrices H2n and H2na might throw more light on the matter.
>>>
>>> Cheers,
>>>
>>> Mike
>>>
>>> Brief details
>>> #~.q=:1024?100000
>>> 1024
>>> ts'hw3 q'
>>> 0.000214 94720
>>> ts'hw q'
>>> 0.000255 79616
>>> ts'Haar q'
>>> 0.000309 34432
>>> ts'applyHaar q'. NB. !!!!!
>>> 0.084739 3.83503e7
>>>
>>>
>>> Sent from my iPad
>>>
>>>> On 7 Jan 2018, at 22:44, Jimmy Gauvin <[email protected]> wrote:
>>>>
>>>> I was playing with alternative ways to get to the result:
>>>>
>>>> hw=: 3 : '; _2 -/\ each (-2^i.1+2^.#y) +/\ each < y'
>>>>
>>>> hw3=: 3 : '_2 -/\ ; (-2^i.1+2^.#y) +/\ each < y'
>>>>
>>>>
>>>> and I noticed that hw is faster even though it has an extra each.
>>>>
>>>> Shouldn't hw3 be faster with less boxing and unboxing ?
>>>>
>>>>
>>>>
>>>> As usual Raul' Haar is faster and consumes less memory.
>>>>
>>>>
>>>>
>>>> On Sun, Jan 7, 2018 at 4:07 PM, 'Mike Day' via Programming <
>>>> [email protected]> wrote:
>>>>
>>>>> Not at all familiar with this, but I played with the Matrices H2n in
>>>>> the
>>>>> wiki
>>>>> article, and found they reproduce these results, subject to
>>>>> reordering.
>>>>>
>>>>> So, fwiw, here are a few verbs:
>>>>> (I assume the argument size is a power (> 0) of 2; no checking!)
>>>>>
>>>>> kp =: *&$ ($,) 0 2 1 3 |: */ NB. Kronecker product from J Wiki
>>>>>
>>>>> H2n =: 3 : 0 NB. H2n as defined in Wiki
>>>>> n =. y
>>>>> h =. 1 1,:1 _1
>>>>> n2=. 2
>>>>> for_i. i. <: 2 <.@^. n do.
>>>>> h =. (h kp ,:1 1) , (=i. n2) kp ,: 1 _1
>>>>> n2=. +:n2
>>>>> end.
>>>>> h
>>>>> )
>>>>>
>>>>> H2na =: 3 : 0 NB. H2n reordered to match Raul's Haar verb
>>>>> i =. i. n =. y
>>>>> I =. ''
>>>>> while. 1 < n =. -: n do.NB. build map from Wiki H2n to required order
>>>>> i =. (-n) <\ i
>>>>> I =. I, ;{: i
>>>>> i =. ;}: i
>>>>> end.
>>>>> (I, |. i) { H2n y
>>>>> )
>>>>>
>>>>> applyHaar =: H2na@# +/ . * ] NB. do matrix product
>>>>>
>>>>> Haar 3 2 3 2 2 2 1 1 NB. Raul's verb
>>>>> 1 1 0 0 0 2 4 16
>>>>> applyHaar 3 2 3 2 2 2 1 1 NB. verb using H2n matrix
>>>>> 1 1 0 0 0 2 4 16
>>>>>
>>>>> Might help a bit!
>>>>>
>>>>> Mike
>>>>>
>>>>>
>>>>>
>>>>>> On 07/01/2018 17:41, Raul Miller wrote:
>>>>>>
>>>>>> https://en.wikipedia.org/wiki/Haar_wavelet
>>>>>>
>>>>>> https://pbs.twimg.com/media/DS7mmaLWAAE6-US.jpg
>>>>>>
>>>>>> Haar=: _2&(-/\ , Haar^:(1<#)@(+/\))
>>>>>> Haar 3 2 3 2 2 2 1 1
>>>>>> 1 1 0 0 0 2 4 16
>>>>>>
>>>>>> Like the fast fourier transform, this is only defined on arguments
>>>>>> whose length is a power of 2 (which might be enforced by requiring an
>>>>>> array whose dimensions are all 2 and then raveling it).
>>>>>>
>>>>>> Unlike the FFT, however, the Haar transform is not self inverting.
>>>>>> Still, since every non-lossey transform deserves an inverse transform:
>>>>>>
>>>>>> iHaar=: -:@(+/,@,.-~/)@(,:$:^:(1<#))/@($~ 2,-:@#)
>>>>>> iHaar 1 1 0 0 0 2 4 16
>>>>>> 3 2 3 2 2 2 1 1
>>>>>>
>>>>>> As an aside, though: I found the wikipedia article nearly useless for
>>>>>> this implementation - I used the example contained in the image to
>>>>>> write this code and only briefly skimmed the wikipedia text for rough
>>>>>> agreement. (And since I've encountered plenty of errors in wikipedia
>>>>>> in the past, that's all I'm inclined to do with that page at the
>>>>>> moment.) Still, if someone with more familiarity with wavelets than I
>>>>>> could tell me if I've screwed up royally (and, if so, how and where),
>>>>>> please let me know.
>>>>>>
>>>>>> Thanks,
>>>>>>
>>>>>>
>>>>> ---
>>>>> This email has been checked for viruses by Avast antivirus software.
>>>>> https://www.avast.com/antivirus
>>>>>
>>>>>
>>>>> ----------------------------------------------------------------------
>>>>> 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
