The Hadamard transform works similarly to the Fourier transform: it takes a vector and maps it to its frequency components, which are the Walsh functions. Instead of sine waves in the Fourier transform, the Walsh functions are discrete "square waves" so to speak. This makes it easier to compute.

The size of a Hadamard matrix is a power of two, 2^n x 2^n. We can build a Hadamard matrix for a given n recursively. $$H_1 = beginpmatrix1&1\1 & -1endpmatrix$$ And then build block matrices in the same pattern to get $$H_n = beginpmatrixH_n-1&H_n-1\H_n-1 & -H_n-1endpmatrix$$

For instance,

$$H_3 = beginpmatrix1&1&1&1&1&1&1&1\1&-1&1&-1&1&-1&1&-1\1&1&-1&-1&1&1&-1&-1\1&-1&-1&1&1&-1&-1&1\1&1&1&1&-1&-1&-1&-1\1&-1&1&-1&-1&1&-1&1\1&1&-1&-1&-1&-1&1&1\1&-1&-1&1&-1&1&1&-1endpmatrix = beginpmatrixH_2&H_2\H_2&-H_2endpmatrix$$

Note how the entries are always 1 or -1 and symmetric! Easy! Then the Hadamard transform is simply multiplying H_n by a column vector and dividing by 2^n.

In order to apply the Hadamard transform to 2D images, we have to apply the 1D HT across every column and then every row. In short, for an input matrix X, the output will be a new matrix Y given by two matrix multiplications:

$$Y = frac12^2nH_nXH_n$$

Another way to implement this is by first applying the HT to every column in the image separately. Then apply it to every row, pretending they were column vectors. For example, to compute the 2D HT of $$beginpmatrix2 & 3\ 2 & 5endpmatrix$$

we first do

$$frac12beginpmatrix1&1\1 & -1endpmatrixbeginpmatrix2 \ 2endpmatrix = beginpmatrix2 \ 0endpmatrix quadquadquad frac12beginpmatrix1&1\1 & -1endpmatrixbeginpmatrix3 \ 5endpmatrix = beginpmatrix4 \ -1endpmatrix quadrightarrowquad beginpmatrix2&4\0 & -1endpmatrix $$

which takes care of the columns, and then we do the rows,

$$frac12beginpmatrix1&1\1 & -1endpmatrixbeginpmatrix2 \ 4endpmatrix = beginpmatrix3 \ -1endpmatrix quadquadquad frac12beginpmatrix1&1\1 & -1endpmatrixbeginpmatrix0 \ -1endpmatrix = beginpmatrix-0.5 \ 0.5endpmatrix quadrightarrowquad beginpmatrix3&-1\-0.5 & 0.5endpmatrix $$

Finally, one can also compute this recursively using the fast Hadamard transform. It's up to you to find out if this is shorter or not.

Challenge

Your challenge is to perform the 2D Hadamard transform on a 2D array of arbitrary size 2^n x 2^n, where n > 0. Your code must accept the array as input (2D array or a flattened 1D array). Then perform the 2D HT and output the result.

In the spirit of flexibility, you have a choice for the output. Before multiplying by the prefactor 1/(2^(2n)), the output array is all integers. You may either multiply by the prefactor and output the answer as floats (no formatting necessary), or you may output a string "1/X*" followed by the array of integers, where X is whatever the factor should be. The last test case shows this as an example.

You may not import the Hadamard matrices from some convenient file that already contains them; you must build them inside the program.

Test Cases

Input:

2,3
2,5

Output:

3,-1
-0.5,0.5

Input:

1,1,1,1
1,1,1,1
1,1,1,1
1,1,1,1

Output:

1,0,0,0
0,0,0,0
0,0,0,0
0,0,0,0

Input:

1,0,0,0
1,1,0,0
1,1,1,0
1,1,1,1

Output:

0.6250,0.1250,0.2500,0
-0.1250,0.1250,0,0
-0.2500,0,0.1250,0.1250
0,0,-0.1250,0.1250

Input:

0,0,0,0,0,0,0,1
0,0,0,0,0,0,2,0
0,0,0,0,0,3,0,0
0,0,0,0,4,0,0,0
0,0,0,5,0,0,0,0
0,0,6,0,0,0,0,0
0,7,0,0,0,0,0,0
8,0,0,0,0,0,0,0

Output:

1/64*
36,4,8,0,16,0,0,0
-4,-36,0,-8,0,-16,0,0
-8,0,-36,-4,0,0,-16,0
0,8,4,36,0,0,0,16
-16,0,0,0,-36,-4,-8,0
0,16,0,0,4,36,0,8
0,0,16,0,8,0,36,4
0,0,0,-16,0,-8,-4,-36

Scoring

This is code golf, so the smallest program in bytes wins.

Haskell, 111 102 97 92 bytes

import Data.Matrix
t x|n<-nrows x=h n*x*h n
h 1=1
h i|m<-h$div i 2=(/2)<$>m<->m<|>(m<->(-m))

This uses the Data.Matrix module. The input matrix has to be of type Matrix Float or Matrix Double.

TIO doesn't provide Data.Matrix, so no link.

Edit: -9 bytes thanks to @ceased to turn counterclockwis, -5 bytes thanks to @H.PWiz, -5 bytes thanks to @xnor

JavaScript (ES6), 114 bytes

This version is golfed a bit further by embedding the function $h$ within the callback of the deepest loop.

M=>(g=z=>M=M.map((r,y)=>r.map((_,x)=>r.map(v=>s+=(h=m=>m?-h(m&~-m):z?v:M[i-1][x])(i++&(z?x:y)),i=s=0)&&s/i)))(g())

Try it online!

JavaScript (ES6),  161 ... 117  116 bytes

M=>(h=m=>m?-h(m&~-m):1,g=z=>M=M.map((r,y)=>r.map((_,x)=>r.map(v=>s+=(z?v:M[i][x])*h(i++&(z?x:y)),s=i=0)&&s/i)))(g())

Try it online!

How?

Each 0-indexed cell $(x,y)$ of a Hadamard matrix $H_n$ is equal to:

$$H_n(x,y)=cases
1,&textif $(xAnd y)$ is evil\
-1,&textotherwise
$$

or:

$$H_n(x,y)=(-1)^operatornamepopcount(x And y)$$

where $And$ is a bitwise AND and $operatornamepopcount(m)$ is the number of bits set in $m$.

The function $h$ computes $(-1)^operatornamepopcount(m)$ recursively:

h = m => m ? -h(m & ~-m) : 1

(This is similar to my answer to Is this number evil?)

The function $g$ computes either $frac12^nH_nM$ or $frac12^nMH_n$, depending on the flag $z$:

g = z => // z = flag cleared on the 1st pass, set on the 2nd
 M = // update M[]:
 M.map((r, y) => // for each row r[] at position y in M[]:
 r.map((_, x) => // for each cell at position x in r[]:
 r.map(v => // for each value v in r[]:
 s += // add to s:
 (z ? v // v = M(i, y) if z is set,
 : M[i][x]) * // or M(x, i) otherwise
 h(i++ & (z ? x // multiplied by H(x, i) if z is set,
 : y)), // or H(i, y) otherwise
 // and increment i afterwards
 s = i = 0 // start with s = i = 0
 ) // end of inner map() over the columns
 && s / i // yield s divided by the width of the matrix
 ) // end of outer map() over the columns
 ) // end of map() over the rows

We invoke it twice to compute $frac12^2nH_nMH_n$.

MATL, 12 bytes

&nKYLw/Gy,Y*

Try it online! Or verify all test cases.

Explanation

Consider input [2,3; 2,5] as an example.

&n % Input (implicit). Push number of rows and number of columns
 % STACK: 2, 2
KYL % Hadamard matrix of specified size
 % STACK: 2, [1 1; 1 -1]
w/ % Swap, divide element-wise
 % STACK: [0.5 0.5; 0.5 -0.5]
G % Push input again
 % STACK: [0.5 0.5; 0.5 -0.5], [2,3; 2,5]
y % Duplicate from below
 % STACK: [0.5 0.5; 0.5 -0.5], [2,3; 2,5], [0.5 0.5; 0.5 -0.5]
, % Do twice
 Y* % Matrix multiply
 % End (implicit)
 % STACK: [3 -1; -0.5 0.5]
 % Display (implicit)