`irfft`: weird when odd

[Inkbottle007]

Notice that the rfft of odd and even length signals are of the same shape. (?) Hum, I'm not sure about that. However when composing irfft with rttf on odd size sample, below is what happens:

(* https://ocaml.xyz/book/signal.html
   rfft composed with irfft in odd case behaves unexpectedly. *)

open Owl
open Owl_plplot
open Core
module G = Dense.Ndarray.Generic

(* Size of sample below is 6. We probe it with 6 different frequencies.  One of them is
   the zero frequency; the others have values arranged in a Hermitian sort of way,
   because we are in the "real" case.  So we only need 5 / 2 = 3 values. `rfft` on that
   sample gives us 4 values: all is well. *)
let even_a = [| 1.; 2.; 1.; -1.; 1.5; 1.0 |] |> Fun.flip Arr.of_array [| 6 |]

(* When the size of the sample is 5, for the exact same reasons as above, `rfft` returns a
   size 1 + 4 / 2 = 3, array, just as expected. Now when we call `irfft` on that it
   doesn't know from which case it comes from. `irfft` assumes it is the even case, so 3
   must be 1 + 3 / 2... Anyway, you'll figure it out. *)
let odd_a = even_a |> G.get_slice [ [ 0; 4 ] ]
let correct_even_rfft = even_a |> Owl_fft.D.rfft |> Owl_fft.D.irfft
let faulty_odd_rfft = odd_a |> Owl_fft.D.rfft |> Owl_fft.D.irfft

let mimicking_faulty_odd_rfft =
  odd_a
  |> G.cast_d2z
  |> Owl_fft.D.fft
  |> G.get_fancy [ L [ 0; 1; 2; 4 ] ]
  |> Owl_fft.D.ifft
  |> Dense.Ndarray.Z.re
;;

And the output:

#use "bin/rfft_odd_is_weird.ml";;
module G = Owl.Dense.Ndarray.Generic
val even_a : Arr.arr = 
  C0 C1 C2 C3  C4 C5 
R  1  2  1 -1 1.5  1 

val odd_a : (float, Bigarray.float64_elt) G.t =
  
  C0 C1 C2 C3  C4 
R  1  2  1 -1 1.5 

val correct_even_rfft :
  (float, Bigarray.float64_elt) Owl_dense_ndarray_generic.t =
  
  C0 C1 C2 C3  C4 C5 
R  1  2  1 -1 1.5  1 

val faulty_odd_rfft : (float, Bigarray.float64_elt) Owl_dense_ndarray_generic.t =
  
       C0      C1       C2      C3 
R 1.70789 2.40844 -0.37367 0.75734 

val mimicking_faulty_odd_rfft : Dense.Ndarray.Z.cast_arr =
  
       C0      C1       C2      C3 
R 1.70789 2.40844 -0.37367 0.75734 

[jzstark]

Thanks a lot for helping to clarify this matter!