Possibly incorrect spherical harmonic transform

[jishnub]

Please see this discourse thread. I'm not entirely sure if this is a bug or if I'm doing something wrong. The dircourse thread did not generate a lot of discussion so I'm posting it here. I'm trying to compute the spherical harmonic transform of r.x, that is the x-component of the radial unit vector, but limited to only one face of the sphere, the other being set to zero. I'm following the example presented in sphere.jl.

julia> f_test(θ, ϕ) = (π/2 <= ϕ <= 3π/2) ? zero(Float64) : Float64(sin(θ)cos(ϕ))
f_test (generic function with 1 method)

julia> function threshold!(A::AbstractArray, ϵ)
           for i in eachindex(A)
               if abs(A[i]) < ϵ A[i] = 0 end
           end
           A
       end
threshold! (generic function with 1 method)

julia> lmax = 3
3

julia> N = lmax + 1
4

julia> θ = ((0.5:N-0.5)/N) .* pi
0.39269908169872414:0.7853981633974483:2.748893571891069

julia> M = 2N-1
7

julia> ϕ = ((0:M-1)/M) .* 2pi
0.0:0.8975979010256552:5.385587406153931

julia> F = [f_test(θ,ϕ) for θ in θ, ϕ in ϕ]
4×7 Array{Float64,2}:
 0.382683  0.238599  0.0  0.0  0.0  0.0  0.238599
 0.92388   0.576029  0.0  0.0  0.0  0.0  0.576029
 0.92388   0.576029  0.0  0.0  0.0  0.0  0.576029
 0.382683  0.238599  0.0  0.0  0.0  0.0  0.238599

julia> using FastTransforms

julia> P = plan_sph2fourier(F);

julia> PA = plan_sph_analysis(F);

julia> V = PA*F;

julia> threshold!(P\V, 400*eps())
4×7 Array{Float64,2}:
  0.888525  0.0  1.0394  0.0  0.417051   0.0  -0.0697631
  0.0       0.0  0.0     0.0  0.0        0.0   0.0
 -0.259657  0.0  0.0     0.0  0.0571638  0.0  -0.018645
  0.0       0.0  0.0     0.0  0.0        0.0   0.0

There are two things that stick out here -- firstly there's no conjugate symmetry despite the function being real. Secondly the last column has a non-zero entry in the third row, which is not what I would have expected given the storage pattern. It's possible that I'm not doing something correctly here, so I'll be glad if someone could point me in the right direction.

[MikaelSlevinsky]

Great question! There's a couple things to unpack here.

• Due to the angular conditional, f_test is not a spherical harmonic. The function is square integrable, so it does have an infinite spherical harmonic series. The spherical harmonic analysis in this library does recover coefficients in an infinite expansion subject to roundoff and aliasing. In this case, I think the aliasing is significant.

• In the discourse thread, one suggestion is to try a genuine spherical harmonic to help you see what's going on. Indeed, if we remove the branch in f_test, then the spherical harmonic is exactly recovered.

• *sph* plans and routines are real-to-real (so that means in phi it's a sine/cosine series) but I recently added complex-valued complex-to-complex transforms through spin-weighted spherical harmonics (which I think only make sense with complex Fourier series) via plan_spinsph2fourier with spin 0.

• You are right: U[3,7] shouldn't really be there. At one point in time, I had a post-processing function that zeroed-out spurious modes but I removed it because in some sense a non-zero mode out in the lower half (with respect to the array's second diagonal) is a symptom that the function has not been resolved to an appreciable precision.

Below, I show that f_test without branch is exactly recovered; I check it's 2-norm, and compare with a much larger sampling of the original f_test, whose 2-norm should be exactly half. I also show how to use the spin-0-weighted plans for complex-to-complex (in case you want to see (conjugate)-symmetry in signals).

julia> using FastTransforms, LinearAlgebra

julia> function threshold!(A::AbstractArray, ϵ)
           for i in eachindex(A)
               if abs(A[i]) < ϵ A[i] = 0 end
           end
           A
       end
threshold! (generic function with 1 method)

julia> f_test(θ, ϕ) = Float64(sin(θ)cos(ϕ))
f_test (generic function with 1 method)

julia> lmax = 3
3

julia> N = lmax+1
4

julia> θ = ((0.5:N-0.5)/N) .* pi
0.39269908169872414:0.7853981633974483:2.748893571891069

julia> M = 2N-1
7

julia> ϕ = ((0:M-1)/M) .* 2pi
0.0:0.8975979010256552:5.385587406153931

julia> F = [f_test(θ,ϕ) for θ in θ, ϕ in ϕ]
4×7 Array{Float64,2}:
 0.382683  0.238599  -0.0851551  -0.344786  -0.344786  -0.0851551  0.238599
 0.92388   0.576029  -0.205583   -0.832387  -0.832387  -0.205583   0.576029
 0.92388   0.576029  -0.205583   -0.832387  -0.832387  -0.205583   0.576029
 0.382683  0.238599  -0.0851551  -0.344786  -0.344786  -0.0851551  0.238599

julia> P = plan_sph2fourier(F)
FastTransforms Spherical harmonic--Fourier plan for 4×7-element array of Float64

julia> PA = plan_sph_analysis(F)
FastTransforms plan for FFTW Fourier analysis on the sphere for 4×7-element array of Float64

julia> V = PA*F
4×7 Array{Float64,2}:
 -9.93899e-17  1.95402e-16   1.77245      -9.08355e-18   6.85223e-17  -1.5761e-17   -1.03887e-16
  0.0          0.0          -7.85046e-17   0.0           2.27253e-17   0.0           0.0
  8.43351e-17  6.65986e-18   0.0           4.07983e-17  -2.23627e-17  -1.52999e-17   4.30315e-17
  0.0          0.0          -1.11022e-16   0.0           9.41315e-18   0.0           0.0

julia> U = P\V
4×7 Array{Float64,2}:
 -1.80314e-16  2.24093e-16   2.04665      -4.3329e-17   1.05784e-16  -1.02674e-17  -1.36464e-16
  0.0          0.0          -4.83183e-17   0.0          7.25671e-18   0.0           0.0
  7.11176e-17  6.57688e-18   0.0           1.83105e-17  2.66502e-17  -1.83837e-17   7.51536e-18
  0.0          0.0          -1.06997e-16   0.0          1.43728e-17   0.0           0.0

julia> threshold!(U, 400eps())
4×7 Array{Float64,2}:
 0.0  0.0  2.04665  0.0  0.0  0.0  0.0
 0.0  0.0  0.0      0.0  0.0  0.0  0.0
 0.0  0.0  0.0      0.0  0.0  0.0  0.0
 0.0  0.0  0.0      0.0  0.0  0.0  0.0

julia> norm(U)^2/π
1.3333333333333335

julia> # Now, we use 1024^2 modes to resolve the original f_test with a discontinuity in one component of its gradient along a great circle

julia> lmax = 1023
1023

julia> N = lmax+1
1024

julia> θ = ((0.5:N-0.5)/N) .* pi
0.0015339807878856412:0.0030679615757712823:3.1400586728019078

julia> M = 2N-1
2047

julia> ϕ = ((0:M-1)/M) .* 2pi
0.0:0.003069460335700824:6.280115846843886

julia> f_test(θ, ϕ) = (π/2 <= ϕ <= 3π/2) ? zero(Float64) : Float64(sin(θ)cos(ϕ))
f_test (generic function with 1 method)

julia> F = [f_test(θ,ϕ) for θ in θ, ϕ in ϕ];

julia> P = plan_sph2fourier(F)
FastTransforms Spherical harmonic--Fourier plan for 1024×2047-element array of Float64

julia> PA = plan_sph_analysis(F)
FastTransforms plan for FFTW Fourier analysis on the sphere for 1024×2047-element array of Float64

julia> V = PA*F;

julia> U = P\V;

julia> threshold!(U, 400eps());

julia> norm(U)^2/π # approximately half the 4/3 from the exact signal above (because the angular conditional cuts the signal in half)
0.6666666670497002

julia> F = ComplexF64[f_test(θ,ϕ) for θ in θ, ϕ in ϕ];

julia> P = plan_spinsph2fourier(F, 0)
FastTransforms Spin-weighted spherical harmonic--Fourier plan for 1024×2047-element array of Complex{Float64}

julia> PA = plan_spinsph_analysis(F, 0)
FastTransforms plan for FFTW Fourier analysis on the sphere (spin-weighted) for 1024×2047-element array of Complex{Float64}

julia> V = PA*F;

julia> U = P\V;

julia> threshold!(U, 400eps());

julia> norm(U)^2/π
0.6666666670496998

[jishnub]

Thanks, this is really helpful!