# Median of complex numbers different from Matlab output

**URL:** <https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352>\
**Category:** General Usage\
**Tags:** question, matlab, complex-numbers\
**Created:** [July 22, 2024, 5:06pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352 "2024-07-22T17:06:09Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![shindelr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shindelr/32/210829_2.png) [@shindelr](https://discourse.julialang.org/u/shindelr)\
**Post date:** [July 22, 2024, 5:06pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/1 "2024-07-22T17:06:10Z")

</div>

Hi there everyone! This is my first ever post in the Julia community and since I’m unfamiliar with the various tags and categories on here, I just put this question in General.

I’m currently in the process of porting a PIV Matlab script over to Julia for a professor at my university. However, neither Julia nor Matlab are in my CS curriculum so it’s slow going finding the ins and outs of both languages!

**The Issue I’m Having**  
The following code is part of a local median filter being iteratively used on small windows of a larger matrix composed of complex numbers and NaN values. I’ve created this little function which just filters out the NaN’s, then takes the median of the real and imaginary parts separately, and recombines them for the return value.

Julia code

```julia
using Statistics

function im_median(collection)
    i = filter(x -> !isnan(x), collection)

    if length(i) > 0
        real_part = median(real(i))
        im_part = median(imag(i))
        return real_part + im_part * im
    end

    # If all NaN
    return NaN
end

```

For your reference, here is the equivalent Matlab code I’m translating. Also, full transparency since this is a port, this function was created as part of the MatPIV package by John Peter Acklam:

```julia
function y = mnanmedian(x)
i = ~isnan(x);
if any(i)
  y = median(x(i));
else
  y = NaN;
end

end

```

I defined this matrix to try to pinpoint the differing outputs. A 3x3 matrix mimics the size of the filter window, basically a point in the grid and all of its neighbors:

```julia
nan_matrix = [
    0.5488135 + 0.71518937im NaN 0.42365480 + 0.64589411im;
    0.43758721 + 0.89177300im 0.96366276 + 0.38344152im NaN;
    0.56804456 + 0.92559664im NaN 0.02021840 + 0.83261985im
]

```

Which yields the following output:

```julia
# Julia output:
0.493200355 + 0.77390461im

# Matlab output:
0.4932 + 0.8035i

```

As you can see, the output for the imaginary end of things differs by enough to cause problems for the PIV algorithm down the road.

Is there something I’m missing with my own implementation of the complex median function?

---

<div class="post-metadata">

**Author:** ![adienes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adienes/32/37459_2.png) [@adienes](https://discourse.julialang.org/u/adienes)\
**Post date:** [July 22, 2024, 5:13pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/2 "2024-07-22T17:13:32Z")

</div>

how does MATLAB compute medians of complex numbers? I see you are doing this

> takes the median of the real and imaginary parts separately,

maybe MATLAB is taking the median via the absolute magnitudes of the numbers rather than separating the components?

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [July 22, 2024, 5:13pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/3 "2024-07-22T17:13:41Z")

</div>

Matlab understands the input as

```julia
>> mnanmedian

nan_matrix =

   0.5488 + 0.7152i NaN + 0.0000i 0.4237 + 0.6459i
   0.4376 + 0.8918i 0.9637 + 0.3834i NaN + 0.0000i
   0.5680 + 0.9256i NaN + 0.0000i 0.0202 + 0.8326i

```

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [July 22, 2024, 5:13pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/4 "2024-07-22T17:13:42Z")

</div>

Hi,

just sharing my thoughts:

```julia
julia> sort(collect(imag.(nan_matrix))[:])
6-element Vector{Float64}:
 0.38344152
 0.64589411
 0.71518937
 0.83261985
 0.891773
 0.92559664

```

So the median clearly is `0.77390461im` (as the mean of 0.71518937 and 0.83261985)

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [July 23, 2024, 1:51am UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/5 "2024-07-23T01:51:45Z")

</div>

Edit: As surmised by @adienes, above.

According to [this reply](https://www.mathworks.com/matlabcentral/answers/96926-why-do-i-receive-different-results-when-applying-median-or-nanmedian-to-my-matrix-vs-a-single-colu#answer_106277) from The Mathworks, when the Matlab `median` function is given a list of complex numbers, they are sorted according to magnitude, and the median of this sorted list is calculated. So the following Julia method for complex arrays should be equivalent to the Matlab function:

```julia
function im_median(collection::AbstractArray{Complex{T}}) where {T}
    i = filter(x -> !isnan(x), collection)
    isempty(i) && return NaN
    sort!(i; by=abs)
    n = length(i)
    no2 = n ÷ 2
    isodd(n) && return i[no2+1]
    return (i[no2] + i[no2+1]) / 2
end

```

It produces the following result on your sample matrix:

```julia
julia> im_median(nan_matrix)
0.493200355 + 0.8034811850000001im

```

which appears to agree with the Matlab result.

---

<div class="post-metadata">

**Author:** ![shindelr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shindelr/32/210829_2.png) [@shindelr](https://discourse.julialang.org/u/shindelr)\
**Post date:** [July 23, 2024, 5:03pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/6 "2024-07-23T17:03:05Z")

</div>

Thanks a ton everyone! This really smoothed out the results from program. Also, I had no idea you could use `&&` or `||` in the `<cond> <bool op> <statement>` format like that. I spend a lot of time while I’m working on this program wondering how true Julia writers might use the language differently.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [July 23, 2024, 5:24pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/7 "2024-07-23T17:24:16Z")

</div>

You also don’t need to sort the entire array — you just need to make sure that the median element(s) are in the right place. That’s precisely the reason for [`partialsort!`](https://docs.julialang.org/en/v1/base/sort/#Base.Sort.partialsort!). There’s surely more ways you could continue to speed this up, but doing less sorting work can be a quick (and huge!) win:

```julia
function im_median_faster(collection::AbstractArray{Complex{T}}) where {T}
    i = filter(x -> !isnan(x), collection)
    isempty(i) && return NaN
    n = length(i)
    v = partialsort!(i, div(n+1, 2, RoundDown):div(n+1, 2, RoundUp); by=abs)
    return sum(v)/length(v)
end

```

```julia-repl
julia> A = rand(Complex{Float64}, 10000);

julia> @btime im_median($A)
  1.717 ms (4 allocations: 312.59 KiB)
0.4918520015608561 + 0.5308805009885738im

julia> @btime im_median_faster($A)
  198.083 μs (4 allocations: 312.59 KiB)
0.4918520015608561 + 0.5308805009885738im

```

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [July 23, 2024, 5:53pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/8 "2024-07-23T17:53:59Z")

</div>

> [@mbauman](#):
>
> There’s surely more ways you could continue to speed this up

One easy way (that I should have used in my first post) is to use `abs2` rather than `abs`. This also gives a big speedup:

```julia
julia> function im_median_faster_still(collection::AbstractArray{Complex{T}}) where {T}
           i = filter(x -> !isnan(x), collection)
           isempty(i) && return NaN
           n = length(i)
           v = partialsort!(i, div(n+1, 2, RoundDown):div(n+1, 2, RoundUp); by=abs2)
           return sum(v)/length(v)
       end
im_median_faster_still (generic function with 1 method)

julia> using BenchmarkTools

julia> A = rand(ComplexF64, 10_000);

julia> @btime im_median($A)
  2.339 ms (4 allocations: 312.59 KiB)
0.49722004577323964 + 0.6006781196428463im

julia> @btime im_median_faster($A)
  345.400 μs (4 allocations: 312.59 KiB)
0.49722004577323964 + 0.6006781196428463im

julia> @btime im_median_faster_still($A)
  94.400 μs (4 allocations: 312.59 KiB)
0.49722004577323964 + 0.6006781196428463im

```

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [July 23, 2024, 6:27pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/9 "2024-07-23T18:27:17Z")

</div>

> [@shindelr](#):
>
> Also, I had no idea you could use `&&` or `||` in the `<cond> <bool op> <statement>` format like that

You probably already found this, but short-circuit evaluation is discussed [in this section](https://docs.julialang.org/en/v1/manual/control-flow/#Short-Circuit-Evaluation) of the Julia documentation.

---

<div class="post-metadata">

**Author:** ![shindelr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shindelr/32/210829_2.png) [@shindelr](https://discourse.julialang.org/u/shindelr)\
**Post date:** [August 1, 2024, 3:41pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/15 "2024-08-01T15:41:30Z")

</div>

Just wanted to report back in with my final solution! Sorting `by=abs2` worked great on almost all cases. In the end though there were a bunch of cases where Matlab was sorting things in a really non-sensical manner and so my median output wasn’t matching. Luckily a Matlab buddy of mine pointed out that Matlab sorts first by absolute magnitude and then by phase angle!

So in the end this worked like a charm:  
`by=x->(abs2(x), angle(x))`

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [August 1, 2024, 4:42pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/16 "2024-08-01T16:42:54Z")

</div>

> [@shindelr](#):
>
> by=x-\>(abs2(x), angle(x))

This may have some issues across the `angle()` function discontinuity at `+/- pi`?

 ![Complex_numbers_median_across_angle_discontinuity](https://global.discourse-cdn.com/julialang/original/3X/2/e/2e5595bfc7e869787733db74c39dca65a3034400.png)

> **Code example**
>
> ```julia
> function im_median(collection::AbstractArray{Complex{T}}) where {T}
> i = filter(x -> !isnan(x), collection)
> isempty(i) && return NaN
> n = length(i)
> v = partialsort!(i, div(n+1, 2, RoundDown):div(n+1, 2, RoundUp); by=x->(abs2(x), angle(x)))
> return sum(v)/length(v)
> end
> 
> Random.seed!(0)
> n = 99
> C = (1 .+ rand(n)) .* exp.(im .* (π .+ (rand(n) .- 0.5)/100))
> Cm = im_median(C)
> 
> using Plots
> scatter(C, c=:blues)
> scatter!((Cm.re, Cm.im), c=:red, ms=6)
> 
> ```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 1, 2024, 5:44pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/17 "2024-08-01T17:44:16Z")

</div>

> [@rafael.guerra](#):
>
> This may have some issues across the `angle()` function discontinuity at `+/- pi`?

The basic issue is that there isn’t a [good total ordering of the complex numbers](https://math.stackexchange.com/questions/487997/total-ordering-on-complex-numbers) that respects “common-sense” algebraic identities.

If your goal is “do whatever Matlab does”, then it indeed [documents that sorting (hence median)](https://www.mathworks.com/help/matlab/ref/sort.html) uses [phase angle](https://www.mathworks.com/help/matlab/ref/angle.html) (in [-\pi, \pi]):

> By default, the `sort` function sorts complex values by their magnitude, and breaks ties using phase angles.

which should match what @shindelr implemented in Julia despite being mathematically arbitrary.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [August 1, 2024, 6:17pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/18 "2024-08-01T18:17:26Z")

</div>

> [@stevengj](#):
>
> despite being mathematically arbitrary.

It’s so arbitrary that I have a hard time imagining any application where it makes sense at all.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [August 2, 2024, 9:21pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/20 "2024-08-02T21:21:41Z")

</div>

It doesn’t even match what you’d get if all the complex values happen to be real since it orders by absolute value. Truly whacky way of ordering complex numbers. Ordering lexicologically by (re, im) seems much better. In particular it would give `real(median(v)) == median(real.(v))` which is something at least.

---

<div class="post-metadata">

**Author:** ![Marc.Cox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marc.cox/32/7514_2.png) [@Marc.Cox](https://discourse.julialang.org/u/Marc.Cox)\
**Post date:** [October 5, 2024, 7:43pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/21 "2024-10-05T19:43:16Z")

</div>

(Brief History) After writing the code and testing I, the same as others in this thread, realized Matlab gives nonsense SORT ordering for REAL functions.  
So I searched Julia discourse again (4th time is a charm), found this thread, and so whipped up the code after Keying off StefanKarpinski comment / NOTE Here:  
StefanKarpinski Stefan\_Karpinski (Viking helmet man to the Rescue !!)  
Steward  
Aug 2nd (2024)  
It doesn’t even match what you’d get if all the complex values happen to be real since it orders by absolute value.  
Truly WHACKY way of ordering complex numbers.  
Ordering LEXICOlogically by (re, im) SEEMS much BETTER.

```julia-auto
function complex_min(x,y)
			  ## treat NaN as being GREATER than any other value and so are sorted to the end of the list.
			   ## So isnan IMMEDIATELY DECIDES which variable to return as LESSER
			  if isnan(x)
				  return y
			  end 
			  if isnan(y)
				  return x
			  end 

			  if real(y) < real(x)	## NOTE Strictly less than
				return y	## x is NOT LT y
			  end 

			  if real(y) > real(x)	## NOTE Strictly greater than
				return x 
              end

			  ## Here IFF real(y) = real(x), so sort by imag
			  if imag(y) < imag(x)	## NOTE Strictly less than
			    return y	## x is NOT LT y
			  else
				return x ## x IS LT y
			  end
		end ## end function

```

## !! TDD UNIT TEST - YAY PASSED !!

tx=NaN + 0.0im; ty=1.0 + 0.0im  
julia\> complex\_min(ty,tx)  
1.0 + 0.0im

## Reverse order SB same answer.

julia\> complex\_min(tx,ty)  
1.0 + 0.0im

```julia-auto
		function complex_isless(x,y)
			    ## treat NaN as being GREATER than any other value and so are sorted to the end of the list.
		         ## So isnan IMMEDIATELY DECIDES which variable to return as LESSER
		        if isnan(x)
					return false
		        end 
		        if isnan(y)
					return true
		        end 

			  if real(y) < real(x)	## NOTE Strictly less than
				return false	## x is NOT LT y
			  end 

			  if real(y) > real(x)	## NOTE Strictly greater than
				return true 
              end

			  ## Here IFF real(y) = real(x), so sort by imag
			  if imag(y) < imag(x)	## NOTE Strictly less than
			    return false	## x is NOT LT y
			  else
				return true ## x IS LT y
			  end
		end ## end function

```

## Marc Cox TDD QA TEST My Example code above Passes some TDD QA TESTS for Sort order of a 1Dim Vector of Complex numbers.

vec\_1 = [-3, NaN, NaN,1+im, -3, 1, 1-im, -3-im, -3+im]  
sort(vec\_1, lt=complex\_isless)

julia\> vec\_1 = [-3, NaN, NaN,1+im, -3, 1, 1-im, -3-im, -3+im]  
9-element Vector{ComplexF64}:  
-3.0 + 0.0im  
NaN + 0.0im  
NaN + 0.0im  
1.0 + 1.0im  
-3.0 + 0.0im  
1.0 + 0.0im  
1.0 - 1.0im  
-3.0 - 1.0im  
-3.0 + 1.0im

julia\> sort(vec\_1, lt=complex\_isless)  
9-element Vector{ComplexF64}:  
-3.0 - 1.0im  
-3.0 + 0.0im  
-3.0 + 0.0im  
-3.0 + 1.0im  
1.0 - 1.0im  
1.0 + 0.0im  
1.0 + 1.0im  
NaN + 0.0im  
NaN + 0.0im

# RE: Stefan Karpinski TDD Test: for a working example test case TRY vec\_2 WARNING try WITHOUT NaN

vec\_2 = [-3, 0.0, 0.5,1+im, -3, 1, 1-im, -3-im, -3+im]  
julia\> real(MEDIAN\_2(vec\_2)) == MEDIAN\_2(real.(vec\_2))  
true

## Marc Cox WARNING NOTE there are cases, e.g. NaN in the list,where it is critical that median always uses a SORTED list, NOT a middle FUNCTION

julia\> real(MEDIAN\_2(sort(vec\_2, lt=complex\_isless))) == MEDIAN\_2(real.(sort(vec\_2, lt=complex\_isless)))  
true

More TDD TESTS to do, to reverify intuition and correct operation of SORT ORDER in ALL CASES suggest try Visualizing Complex function SORT ordering with Julia DomainColoring.jl: Smooth Complex Plotting  
at [Home · DomainColoring.jl](https://eprovst.github.io/DomainColoring.jl/stable/)

---

<div class="post-metadata">

**Author:** ![Marc.Cox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marc.cox/32/7514_2.png) [@Marc.Cox](https://discourse.julialang.org/u/Marc.Cox)\
**Post date:** [December 2, 2024, 2:14pm UTC](https://discourse.julialang.org/t/median-of-complex-numbers-different-from-matlab-output/117352/22 "2024-12-02T14:14:36Z")

</div>

This post was temporarily hidden by the community for possibly being off-topic, unfocused, inappropriate, or spammy.