CODEWITHSUNDEEP
mathimage-processingsignal-processinghaar-wavelet
user2372976
user2372976

Reputation: 725

Create NxN Haar Matrix

I can't find a definition for generating an un-normalized NxN Haar matrix. So what is the equation?

see: http://en.wikipedia.org/wiki/Haar_wavelet

Thx, Chris

Upvotes: 0

Views: 8027

Answers (5)

cmarkle
cmarkle

Reputation: 1

The answer by Nikolai Janakiev is close but its not properly normalized -- for n > 4, HH^T != 1. I brute forced the proper normalization below.

def haarMatrix(n, normalized=False):
    # Allow only size n of power 2
    l = np.log2(n)
    assert l.is_integer(), "Only powers for 2 are permitted."
    n = 2**np.ceil(l)
    if n > 2:
        h = haarMatrix(n // 2)
    else:
        return np.array([[1, 1], [1, -1]])

    # calculate upper haar part
    h_n = np.kron(h, [1, 1])
    # calculate lower haar part
    h_i = np.kron(np.eye(len(h)), [1, -1])
    # combine parts
    h = np.vstack((h_n, h_i))

    if normalized:
        h /= np.sqrt(np.linalg.norm(h, ord=1, axis=-1, keepdims=True))

    return h

Upvotes: 0

Hakan Demir
Hakan Demir

Reputation: 327

def HT(N):
    if N == 1: return np.array([[1]])
    return 1/np.sqrt(2)*np.concatenate(
        (
            np.kron(HT(N//2),[1,1])
        ,
            np.kron(np.identity(N//2),[1,-1])
        ),axis = 0
        )

I used the equation in Wikipedia adding a normalization factor 1/sqrt(2). The matrix becomes orthonormal. Does it work?

Upvotes: 0

Nikolai Janakiev
Nikolai Janakiev

Reputation: 366

Here is the algorithm for both the normalized and un-normalized Haar matrix based on the recursive formula for the Haar matrix:

from the paper "Discrete wavelets and perturbation theory" by W-H Steeb, et al. Here is the implementation in Python

def haarMatrix(n, normalized=False):
    # Allow only size n of power 2
    n = 2**np.ceil(np.log2(n))
    if n > 2:
        h = haarMatrix(n / 2)
    else:
        return np.array([[1, 1], [1, -1]])

    # calculate upper haar part
    h_n = np.kron(h, [1, 1])
    # calculate lower haar part 
    if normalized:
        h_i = np.sqrt(n/2)*np.kron(np.eye(len(h)), [1, -1])
    else:
        h_i = np.kron(np.eye(len(h)), [1, -1])
    # combine parts
    h = np.vstack((h_n, h_i))
    return h

Upvotes: 4

user2372976
user2372976

Reputation: 725

Thx all. Wikipedia gives the 'equation':

I've wrote a recursive solution for generating an un-normalized NxN Haar matrix in octave.

function [h] = haar(n)
h = [1];
if n > 2
    h = haar(n/2);
endif
% calculate upper haar part
h_n = kron(h,[1,1]); 
% calculate lower haar part 
h_i = kron(eye(length(h)),[1,-1]);
% combine parts
h = [h_n; h_i];
endfunction

disp(haar(8));

Upvotes: 1

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

It depends on what exactly you want to achieve. The Haar matrix is the 2x2 DCT matrix, so inversly, you can treat the NxN DCT(II) matrix as the Haar matrix for that block size.

Or if the N is dyadic, N=2^n, then you might be asking for the transform matrix for n stages of the Haar transform. Which might be a problem because of the sampling rate decimation in each step.

Upvotes: 1

Related Questions