Thanks for your code. However, it does not really help since it is restricted to 1D.
However, I’m currently working on a function make_hermitian to modify the array in such a way, that it is hermitian after padding and also preserves Parseval’s theorem.
The idea I pasted it for was to show what’s needed in order to do it right.
It is common knowledge how to do it. You can find it in many places over the net.
Yet few explains why the halving is done beside symmetry.
So the reason is when you take both constraints.
The actually reason comes form the understanding that the DFT is basically a DFS (Discrete Fourier Series).
In Fourier Series analysis we all remember the value of the rectangle window at the edges, it is half the value of the discontinuity jump. So in this case, it is 0.5 :-).
You have theory in place and you can do coding…
I finished coding of the make_hermitian function.
The name is not correct in strict mathematical notation, but it basically does, what you expect in 2D.
So we pad the array with extra rows, cols,… to achieve the symmetry property of real-valued FFTs.
If someone interested in the Math behind it, I just derived it as an answer for someone - \operatorname{sinc} \left( \cdot \right) Interpolation Based on DFT / FFT. It might be useful as a reference?
Thanks for your explanations.
However, I tried to implement your conclusions, but encountered a few problems.
Let’s discuss the case of having an array with odd length n and that we downsample it to a length n-1. lf I understand you procedure correctly, the highest positive frequency is simply removed:
julia> using FFTW
julia> function ds(y)
a = fftshift(fft(y))
b = a[1:end-1]
return ifft(ifftshift(b))
end
ds (generic function with 1 method)
julia> ds([1,2,3,4,5])
4-element Array{Complex{Float64},1}:
1.875 - 0.20307481014556647im
2.654522599411033 + 0.20307481014556647im
4.375 - 0.20307481014556647im
6.095477400588967 + 0.20307481014556647im
So the result is not purely real (as expected?).
However, if we add the highest positive frequency to the highest negative one:
julia> function ds_add(y)
a = fftshift(fft(y))
b = a[1:end-1]
b[1] += a[end]
return ifft(ifftshift(b))
end
ds_add (generic function with 1 method)
julia> ds_add([1,2,3,4,5])
4-element Array{Complex{Float64},1}:
1.25 + 0.0im
3.279522599411033 + 0.0im
3.75 + 0.0im
6.720477400588967 + 0.0im
The result is purely real.
I created a function downsample doing this procedure in N dimensions. So far it always returns a purely real result. In 2D (sampling from a odd N to N-1) it looks (simplified) like this:
julia> function ds_2D(arr)
a = fftshift(fft(arr))
a[:, 1] += a[:, end]
a[1, :] += a[end, :]
b = a[1:end-1, 1:end-1]
return ifft(ifftshift(b))
end
ds_2D (generic function with 1 method)
julia> x = randn((7,7))
7×7 Array{Float64,2}:
1.28566 0.347374 -0.327131 -0.811201 -0.0554991 -1.38513 -1.15124
0.342884 0.743379 -0.7558 -0.038343 -0.0961781 1.80734 1.07544
0.566522 0.129537 -1.06187 -0.280687 0.0809385 -0.934195 0.149161
1.63082 0.578416 -0.303113 -0.832909 0.0581663 1.80572 -0.29874
0.606241 -0.955447 0.820184 1.27619 -0.143522 -0.716512 0.949464
-0.552843 1.11347 2.24364 -0.711459 0.0939156 0.443312 -0.214561
-0.789011 1.18366 -0.399333 -0.279353 0.588281 -1.69753 -1.83028
julia> ds_2D(x)
6×6 Array{Complex{Float64},2}:
1.74993+0.0im 0.217544+0.0im -0.72932+0.0im -0.646827+0.0im -0.919558+0.0im -2.11846+0.0im
0.153758+0.0im 0.443005+0.0im -1.08628+0.0im 0.413973+0.0im 0.820712+0.0im 2.06099+0.0im
1.6024+0.0im 0.236729+0.0im -1.22807+0.0im -0.73624+0.0im 0.190207+0.0im -0.768012+0.0im
1.60426+0.0im -0.723955+0.0im 0.6624+0.0im 0.0153613+0.0im 1.13091+0.0im 0.700818+0.0im
-0.151597+0.0im 0.871951+0.0im 2.37298+0.0im -0.516958+0.0im -0.309972+0.0im 0.233701+0.0im
-1.35343+0.0im 2.15146+0.0im -0.918603+0.0im 0.552748+0.0im -0.521778+0.0im -2.15894+0.0im
julia> ds_2D(x) |> fft |> fftshift
6×6 Array{Complex{Float64},2}:
8.76707+0.0im 2.16954-1.31279im 4.13384+6.00899im -4.59698+0.0im 4.13384-6.00899im 2.16954+1.31279im
7.66118-0.173246im 10.8922-5.08887im 4.52256-1.21773im -0.234256+7.15146im -1.03317+0.540592im -2.43672-2.49465im
-5.51836+0.708964im 2.0933+7.65555im 8.69814+7.08368im -6.45624+1.60355im -8.06179+1.61042im 2.16754-1.71655im
2.84007+0.0im 2.38215+5.68463im 5.36487+3.4028im 3.29783+0.0im 5.36487-3.4028im 2.38215-5.68463im
-5.51836-0.708964im 2.16754+1.71655im -8.06179-1.61042im -6.45624-1.60355im 8.69814-7.08368im 2.0933-7.65555im
7.66118+0.173246im -2.43672+2.49465im -1.03317-0.540592im -0.234256-7.15146im 4.52256+1.21773im 10.8922+5.08887im
At the moment, it’s just a observation without any proof. If we look at the FFT of a random real matrix, we also see the same symmetry:
julia> fft(fftshift(randn((6,6))))
6×6 Array{Complex{Float64},2}:
-5.68469+0.0im -3.65677+0.370793im -1.80369+8.03579im -3.75849+0.0im -1.80369-8.03579im -3.65677-0.370793im
3.46143-5.19503im 1.6603+2.00749im 6.96219+3.45326im 2.1068-2.61896im -2.25073-0.424612im -2.34791-9.26965im
-0.272129-6.99491im 0.547095+3.00501im 12.6262+3.82233im 2.85939+5.94222im 6.08718-3.48567im -4.73672-6.57986im
6.47331+0.0im 5.59877+3.68285im 2.38269+3.51618im -0.935849+0.0im 2.38269-3.51618im 5.59877-3.68285im
-0.272129+6.99491im -4.73672+6.57986im 6.08718+3.48567im 2.85939-5.94222im 12.6262-3.82233im 0.547095-3.00501im
3.46143+5.19503im -2.34791+9.26965im -2.25073+0.424612im 2.1068+2.61896im 6.96219-3.45326im 1.6603-2.00749im
We can see that (1,1), (4, 1), (4,4) and (1,4) are purely real. This is achieved if we add the highest negative and highest positive frequency, and store it at the position of the highest negative one.
Where do these disagreements to your conclusions come from?
@roflmaostc thanks a lot for FFTInterpolations.jl! Now I can delete my ad-hoc code specialized on dimension and evenness. Please register it!
I am not sure I got your message correctly.
But for downsampling you might do one of 3 things:
If I get you right, you tried 3 (Which is what’s posted as a sample code in my answer).
I will add to my answer the code for the other 2.
But to answer your questions, let’s say the input has n samples and you want to interpolate it to m < n samples.
What you can do, as in (1) above is:
subSampleFctr = floor(numSamples / numSamplesO) + 1;
numSamplesO = subSampleFctr * numSamplesO;
The upsample to numSamplesO and then decimate by vX = vX(1:subSampleFctr:numSamplesO);.
Regarding your code.
My code supports Complex and Real signals.
When I know the output should be Real I use ifft(SincInterpolationDft(fft(vX), numSamplesOutput), 'symmetric');. Since for real signals the Nyquist sample must be real you can either force specifically that sample or just on the inverse DFT use the symmetric property.
For the Nyquist sample (The one which I multiply by 0.5) taking only the real part is equivalent to what you did (You added a conjugate term).
The Fourier Series equation, for real signals, allows using half the data with cos() instead of complex numbers. Yet it doesn’t have all the optimizations of the fft() / ifft() so I’m not sure it will be faster.
I plan to do so, but I want to be sure that it’s doing the high frequency content correctly.
I’m confused by what your problem is. Is it that you get a complex signal as output for a real input? If so, you can just discard the imaginary part, which is equivalent to imposing hermiticity of the FFT array (so, for an odd array, symmetrizing the high/low frequencies, and for an even array, zeroing out the lone high frequency). Of course you lose information this way (eg parseval is no longer valid) but the loss should be minimal if the underlying data is smooth and sampled sufficiently.
Thanks for your answer!
Yes, I want to do downsampling in frequency domain because that doesn’t introduce any aliasing (as in comparison to taking every subSampleFctr-th sample).
@antoine-levitt
Yes this what I’m struggling with.
For doing the downsampling does the FFT of an even length really adds the high positive and high negative frequencies together? Otherwise I don’t understand where the 0.0im come from.
In my downsampling (see ds_add) procedure I’m currently doing this and it seems to be fine and delivers correct results.
So maybe as summary:
I artificially reconstruct hermitian property by halving the negative highest frequency boundary, conjugate&mirror it and glue it to the positive highest frequency. Result is purely real.
You imply that this is the same as taking the real part after performing iffting. Is there a proof somewhere?
In downsampling to even length, I take the highest positive frequency entry and add it to the highest negative frequency entry (without conjugate&mirror).
This result is (after iffting) purely real but is also different from that where taking the real part only.
So for upsampling I guess, I understood what’s needed.
For downsampling I’m not sure.
the real part of the IFFT is the IFFT of the even part, and the imaginary part of the IFFT is the IFFT of the odd part, where even and odd are understood with respect to the “time-reversal” symmetry x_{-k} = x_k^*. So by linearity if you take the real part of the IFFT, you take the IFFT of the even part of the input.
But my signal is not guaranteed to be even, is it?
are you referring to the even-ness of the length of the signal? It doesn’t matter. The even-ness I mean is even-ness wrt the symmetry I mentioned: an array is even if x_{-k} = x_k^* and odd if x_{-k} = -x_k^*.
I’m reffering to the symmetry eveness. My real signal (even/odd length) does not have a prior symmetry, right?
So therefore I can’t really follow you 
julia> a = randn(ComplexF64, 3)
3-element Array{Complex{Float64},1}:
-0.2883118512559681 - 0.8504598493836352im
-0.29301253729827786 - 0.5871569599861151im
-0.23443071019678133 + 0.5328734013167514im
julia> e = [a[1]+a[1]' a[2]+a[3]' a[2]'+a[3]]/2 # even part
1×3 Array{Complex{Float64},2}:
-0.288312+0.0im -0.263722-0.560015im -0.263722+0.560015im
julia> ifft(e)
1×3 Array{Complex{Float64},2}:
-0.271918+0.0im 0.315128+0.0im -0.331522+0.0im
julia> real(ifft(a))
3-element Array{Float64,1}:
-0.27191836625034244
0.315128172796569
-0.3315216578021947
Hm, I guess it’s clear.
Do you have any clue regarding my latest question in the Downsampling?
To get the same result as the real part, you must symmetrize/hermitianize the array, which in the case of even-sized arrays means zeroing out the highest frequency. That’s what I wrote above:
(so, for an odd array, symmetrizing the high/low frequencies, and for an even array, zeroing out the lone high frequency)
Hm, I observe something different 
The operations in ds_2D seems to produce the same. zeroing the highest frequency produces something different.
julia> using FFTInterpolations
julia> function ds_2D(arr)
a = fftshift(fft(arr))
a_new = copy(a)
a_new[:, 1] += a[:, end]
a_new[1, :] += a[end, :]
b = a_new[1:end-1, 1:end-1]
b[:, 1] .*= 0.5
b[1, :] .*= 0.5
b[1, 1] = 0.5 .* (a[1,1] .+ a[end, end])
return ifft(ifftshift(b))
end
ds_2D (generic function with 1 method)
julia> function ds_2D_r(arr)
a = fftshift(fft(arr))
b = a[1:end-1, 1:end-1]
return real(ifft(ifftshift(b)))
end
ds_2D_r (generic function with 1 method)
julia> x = randn((7,7));
julia> x_ft = fftshift(fft(x))
7×7 Array{Complex{Float64},2}:
4.78581-2.04962im 3.53796+3.86416im 7.97549+2.14692im 2.33653-0.85103im 7.12309+0.863984im 1.91516-0.566075im -3.83162-3.00291im
2.2941+6.30089im 1.76841+5.39204im 5.67055-13.9238im 1.84566-6.2981im -9.94878+3.93367im 0.746254+1.35699im 0.38076+0.818838im
0.131401-3.36596im 2.09302+1.08541im -1.29329-7.21876im 10.6937-4.72555im -3.27202-0.188586im 2.89425+0.0890002im 7.62254-2.30664im
4.03318-12.5314im 0.0435514+2.99394im -2.02828-4.23962im 3.86994+0.0im -2.02828+4.23962im 0.0435514-2.99394im 4.03318+12.5314im
7.62254+2.30664im 2.89425-0.0890002im -3.27202+0.188586im 10.6937+4.72555im -1.29329+7.21876im 2.09302-1.08541im 0.131401+3.36596im
0.38076-0.818838im 0.746254-1.35699im -9.94878-3.93367im 1.84566+6.2981im 5.67055+13.9238im 1.76841-5.39204im 2.2941-6.30089im
-3.83162+3.00291im 1.91516+0.566075im 7.12309-0.863984im 2.33653+0.85103im 7.97549-2.14692im 3.53796-3.86416im 4.78581+2.04962im
julia> a = ds_2D_r(x);
julia> a_ft = fftshift(fft(a))
6×6 Array{Complex{Float64},2}:
4.78581+0.0im 2.72656+2.21512im 7.54929+0.641469im 2.33653+0.0im 7.54929-0.641469im 2.72656-2.21512im
1.33743+3.55986im 1.76841+5.39204im 5.67055-13.9238im 1.84566-6.2981im -9.94878+3.93367im 0.746254+1.35699im
3.87697-2.8363im 2.09302+1.08541im -1.29329-7.21876im 10.6937-4.72555im -3.27202-0.188586im 2.89425+0.0890002im
4.03318+0.0im 0.0435514+2.99394im -2.02828-4.23962im 3.86994+0.0im -2.02828+4.23962im 0.0435514-2.99394im
3.87697+2.8363im 2.89425-0.0890002im -3.27202+0.188586im 10.6937+4.72555im -1.29329+7.21876im 2.09302-1.08541im
1.33743-3.55986im 0.746254-1.35699im -9.94878-3.93367im 1.84566+6.2981im 5.67055+13.9238im 1.76841-5.39204im
julia> b = ds_2D(x);
julia> b_ft = fftshift(fft(b))
6×6 Array{Complex{Float64},2}:
4.78581+0.0im 2.72656+2.21512im 7.54929+0.641469im 2.33653+0.0im 7.54929-0.641469im 2.72656-2.21512im
1.33743+3.55986im 1.76841+5.39204im 5.67055-13.9238im 1.84566-6.2981im -9.94878+3.93367im 0.746254+1.35699im
3.87697-2.8363im 2.09302+1.08541im -1.29329-7.21876im 10.6937-4.72555im -3.27202-0.188586im 2.89425+0.0890002im
4.03318+0.0im 0.0435514+2.99394im -2.02828-4.23962im 3.86994+0.0im -2.02828+4.23962im 0.0435514-2.99394im
3.87697+2.8363im 2.89425-0.0890002im -3.27202+0.188586im 10.6937+4.72555im -1.29329+7.21876im 2.09302-1.08541im
1.33743-3.55986im 0.746254-1.35699im -9.94878-3.93367im 1.84566+6.2981im 5.67055+13.9238im 1.76841-5.39204im
julia> a_ft ≈ b_ft
true
Woops sorry, hermitianizing in the even case doesn’t mean it has to be zero, it means it has to be real.
using FFTW
data_in = randn(5)
fft_in = fft(data_in)
fft_down = [fft_in[1]; fft_in[2]] # downsampling from 5 to 2
fft_down_hermitianized = [fft_in[1]; real(fft_in[2])]
display(real(ifft(fft_down)))
display(ifft(fft_down_hermitianized))
For downasmpaling my code is perfect.
Just if you have knolwedge the input signal is real you may:
symmetric flag for the fft().Either of those and you’re set.