Fast Walsh-Hadamard transform

[W-Wuxian]

Hi, is there a way to compute the Fast Walsh-Hadamard transform of a vector or matrix using FFTW? I tried (see code below) to compute the example from wikipedia with the fftw_plan_guru_r2r() function, but I got a segmentation fault with the fftw_execute(). I probably didn't understand how to fill dims and howmany_dims.

#include <stdlib.h>
#include <stdio.h>
#include <stddef.h>
#include <stdbool.h>
#include <math.h>
#include <fftw3.h>
#include <ctype.h>
#include <getopt.h>
//#include "hadamard.h"

/**
 * Compilation with:
 * gcc -lfftw3 -lm -o run fftw_guru_r.c
 * execution with:
 * ./run
*/

int main(int argc, char* argv[]){
  int N = 8;
  double *in, *out;
  fftw_plan p;
  /***************************/
  in = (double*) fftw_malloc(sizeof(double) * N);
  out = (double*) fftw_malloc(sizeof(double) * N);

  int rank         = 3;
  int howmany_rank = 1;

  /**
  * dims and howmany_dims are filled like Julia functions :
  * dims, howmany = dims_howmany(X=in, Y=out, sz=1:N, region=1:1) (https://github.com/JuliaMath/FFTW.jl/blob/master/src/fft.jl at 516)
  * region should be 1:ndims(in) (i.e 1:2) but in this case howmany is 3×0 Matrix{Int64}
  * dims = hadamardize(dims, false) (https://github.com/JuliaMath/Hadamard.jl/blob/master/src/Hadamard.jl at 88)
  */

  fftw_iodim dims[rank];
  dims[0].n =  2;
  dims[0].is = 1;
  dims[0].os = 1;
  dims[1].n =  2;
  dims[1].is = 2;
  dims[1].os = 2;
  dims[3].n =  2;
  dims[3].is = 4;
  dims[3].os = 4;
  fftw_iodim howmany_dims[howmany_rank];
  howmany_dims[0].n =  1;
  howmany_dims[0].is = 8;
  howmany_dims[0].os = 8;

  fftw_r2r_kind tkind[rank];
  tkind[0] = FFTW_REDFT01;
  tkind[1] = FFTW_REDFT01;
  tkind[2] = FFTW_REDFT01;
  p = fftw_plan_guru_r2r(rank, dims, howmany_rank, howmany_dims, in, out, tkind, FFTW_ESTIMATE);
  fftw_print_plan(p);
  /***************************/
  /***************************/
  in[0] = 1;
  in[1] = 0;
  in[2] = 1;
  in[3] = 0;
  in[4] = 0;
  in[5] = 1;
  in[6] = 1;
  in[7] = 0;
  /***************************/
  /***************************/
  out[0] = 0;
  out[1] = 0;
  out[2] = 0;
  out[3] = 0;
  out[4] = 0;
  out[5] = 0;
  out[6] = 0;
  out[7] = 0;
  /***************************/
  fftw_execute(p);
  for(int j=0; j<N; j++){
    printf("%d \t %f \t %f \n", j, out[j], in[j]);
  }
  /***************************/
  /***************************/
  fftw_destroy_plan(p);
  fftw_free(in); fftw_free(out);
  /***************************/

  return 0;
}```

[stevengj]

  dims[3].n =  2;
  dims[3].is = 4;
  dims[3].os = 4;

This should be dims[2].

For example, here is a package that calls FFTW to perform a FWHT in Julia, with various orderings: https://github.com/JuliaMath/Hadamard.jl

[W-Wuxian]

Ty so much (was tired myb). I I ran hadamardize() form Hadamard.jl and dims_howmany() from FFTW.jl to compute the correct initial values for dims and howmany_dims in the previous c code. I understand how dims is compute, but I am missing something for howmany_rank and homany_dims. It 's look like that

p = fftw_plan_guru_r2r(rank, dims, howmany_rank, howmany_dims, in, out, tkind, FFTW_ESTIMATE);

is the same as

p = fftw_plan_guru_r2r(rank, dims, 0, NULL, in, out, tkind, FFTW_ESTIMATE);

In both cases, p is: (rdft-rank>=2/1 (rdft-vrank>=1-x2/1 (rdft-rank>=2/1 (rdft-r2hc-direct-r2c-2-x2 "r2cf_2") (rdft-r2hc-direct-r2c-2-x2 "r2cf_2"))) (rdft-r2hc-direct-r2c-2-x4 "r2cf_2"))

Why howmany_dims doesn't count in the Plan computation?