I would argue that the correct symmetry handling is exactly that which results in zeroing out the imaginary part of the result, without changing its real part. The geometrical way of seeing this is that you’ve got the subspace X of purely real real-space signals, which is dual by the FFT to the subspace Y of frequency-space sequences satisfying a particular symmetry. If you’ve got a frequency-space signal y that doesn’t belong to Y, you can project it in that subpsace, ie find the ytilde in Y that is closest to y in the least-squares. By Parseval, that operation corresponds in real space to projecting in the least-squares sense the ifft of y on X, which is exactly taking the real part. Again, the nice thing with this is that 1) you don’t have to make a choice for how you fix the nyquist components, that choice is given by the least-squares projection 2) you don’t even have to make explicit what Y is, or what the projection on Y is. Of course you could write it down explicitly, but why would you bother doing that?
I’m bothering because taking the real part of the iffting result produces worse results than my current routine:
begin
using Revise, FFTInterpolations, FFTW
end
begin
x = range(-10.0, 10.0, length=129)[1:end-1]
y = x'
arr = abs.(x) .+ abs.(y) .+ sinc.(sqrt.(x .^2 .+ y .^2))
end
function ds_r(y)
a = fftshift(fft(y))
b = a[1:end-1, 1:end-1]
b[1:end, 1] .*= 2
b[1, 1:end] .*= 2
b[1, 1] /= 2
return real(ifft(ifftshift(b))) .* length(b) ./ length(y)
end
# 19167.48880005109
sum(abs2.(arr))
sum(abs2.(arr .- ds_r(sinc_interpolate(arr, (129, 129)))))
# 1.1928437627612847e-5
sum(abs2.(arr .- downsample(sinc_interpolate(arr, (129, 129)), (128, 128))))
# 5.592278121264315e-27
However, what is more in accordance, if I modify ds_r:
function ds_r(y)
a = fftshift(fft(y))
b = a[1:end-1, 1:end-1]
b[1:end, 1] .*= 2
b[1, 1:end] .*= 2
b[1, 1] /= 2
return real(ifft(ifftshift(b))) .* length(b) ./ length(y)
end
sum(abs2.(arr .- ds_r(sinc_interpolate(arr, (129, 129)))))
# 5.483458710517112e-27
sinc_interpolation makes the array hermitian and divides the boundaries by 2 (as @RoyiAvital did on StackExchange.
OK so your point is that upsampling then downsampling a real signal has to produce the same result? Indeed the naive scheme of ignoring symmetry and then taking real parts does not preserve this property: if eg upsampling from 2 to 3, the FFT [a b] will get upsampled to [a b 0], but then taking the real part will get [a b/2 b*/2], which will get downsampled to [a b/2]. I don’t think it’s a big deal: the error is of the order of the amplitude of the Nyquist component, which should be small - note that in your experiment you interpolated a function with discontinuous derivatives, hence the relatively large error even at high N. If you really really want to preserve that property indeed it makes sense to look for fancy schemes, but I would be surprised if you didn’t lose other properties in the process.
Yeah after some time I thought I would be a nice property. Otherwise it is hard to test the functions.
The discontinuity is unfortunate, but the result was similar for the pure sinc.
I’m satisfied with the current state of FFTInterpolations.jl
And I leave it like this for now.
Thanks for your input!
If you’re downsmapling from N to M<N samples and M is an even number all needed is to keep the first (In the fftshift() grid) M/2 element (At the index -M/2), multiply it by 0.5 and take only its real part. You remove out the other new nyquist sample (At the M/2 index) and all other smaples with indes greater than M/2 then do regular fft() and the result should be real (Up to floating point errors).
Just as in my code, just the real part of the negative nyquist sample.
@Arrigo_Benedetti, thanks for sharing the Candocia and Principe (1998) reference. It addresses 1D-sinc interpolation of periodic band-limited complex series but it might be extended to more dimensions. It may complement @roflmaostc’s FFTResampling.jl package, by allowing interpolating a coarsely regularly sampled array at a set of arbitrary points within the support of the input.
Some examples below for 1D-sinc interpolation of real and complex series at arbitrary points:
The example below of 2D-sinc interpolation uses two-passes, along x and y, of the above Candocia’s 1D-sinc interpolation. It includes also a nice resampling result from FFTResampling.jl:
The basic Julia code with examples and annotations takes some 150 lines and for that reason was not posted it here. Feel free to ping me if it might be of interest.
I was also considering to rename it and you gave men the final push 
FFTResampling.jl
Very nice work! I love the 3D complex plot …
A really neat package that I now use to interpolate my periodic functions (not in the field of signal processing though). I’d say it is a perfect fit for my package which also uses FFTW internally. Thank you for your great work!
To my use-case, I will do either interpolation or downsampling in a single optimization run. I notice that sinc_interpolate and downsample methods have similar (or same?) interfaces. I thus wonder is it possible to have a, for example, resample method that combines the two? Currently, I use
resample = N_new > N_old ? sinc_interpolate : downsample
A_new = resample(A_old, N_new)
Yes, thanks for your suggestion that’s definitely possible. It’s just happened to be so without much considerations 
I’ll suggest you leave the workaround and I will release a breaking 0.3.0 soon!
I thought about your wish.
We can even generalize resample to handle both down and upsampling at the same time!
julia> x = [1 2 3; 4 5 6; 7 8 9]
3×3 Array{Int64,2}:
1 2 3
4 5 6
7 8 9
julia> size(x)
(3, 3)
julia> resample([1 2 3; 4 5 6; 7 8 9], (4, 2))
4×2 Array{Float64,2}:
1.0 3.0
2.26795 4.26795
7.0 9.0
5.73205 7.73205
The new code is here
It is not very elegant, but I basically do first a upsampling of all dimensions which are larger, and then I do the final downsampling step. In that way we can mix the operations in one step.
I’ll try to clean everything (examples etc.) and then tag a new version by the end of the week.
We now implemented FourierTools.jl a first version of which will appear soon. I want to use this discussion to outline how a proper padding and cutting in Fourier space (also mixed versions) needs to be done.
- We want: Reversibility for cases where we do not reduce sizes.
- We need to not change the integrated power spectrum
The latter part is crucial when it comes to dealing with even sized dimensions in FFT and RFFT of real valued data alike. The lowest (biggest negative) frequency stores in these cases the sum of amplitude of the negative highest frequency PLUS the corresponding value at the aliased positive highest frequency (e.g. [-kxmax, ky, kz] + [+kxmax, ky, kz]). Note that this value is generally NOT real except for special choices of ky and kz.
When changing the region of interest in FFT- (or RFFT-) space the algorithm consists of two steps:
fix_dim_before(Array, dim)-
fix_dim_after(Array, dim)
Both of these routines need to be run for each dimension consecutively. The first before the size change and the second after it.
fix_dim_before(Array, dim) is only active for a dimension if it is even-sized and enlarged
In this case the highest frequency subslice(dim) needs to be halfed in value and this half also copied (without conjugation) to the equivalent positive frequency (e.g. [+kxmax, ky, kz]).
fix_dim_after(Array, dim) is only active for a dimension if it is even-sized and reduced
Here we need to copy the lowest frequency which is just about dropped (e.g. [+kxmax, ky, kz]) and sum it to the other side (here [-kxmax, ky, kz]))
We implemented this for array views (overloading getindex, such that no extra memory allocation is needed.
A remark on just taking the real part and ignoring such (rather complicated) actions:
An even lowest frequency in FFTs (and RFFTs) store both, the positive and negative highest frequencies. Going back to real-space after padding in Fourier-space without copying will then drop some of the energy at the moment when the real part is taken. When reverting this situation (cutting away the previously introduced zeros) it is clear to see that this real-space energy is still lost.
An interesting consequence is that also some other functions need modifications on which we are currently working: For example: multiplication with a complex exponential needs a modification to the lowest frequency along an even dim.
It seems to be that also calculating the Energy or Power spectrum needs a modification to such frequencies: sum((a+a)^2) !== 2*sum(a^2)
The should be fairly easy to check via a comparison to the real-space integral of the square of the signal.
If array sizes are fairly large, all of these effects are small, but wouldn’t it be nice to get this correct?
This is what FourierTools.jl will be trying to achieve.
I have to correct myself: The lowest frequency in an FFT does NOT represent the sum of this frequency content and the corresponding (aliased) frequency from the positive max frequency, but it corresponds to the AVERAGE. This is easily seem by performing an FFT of a delta-shaped distribution at zero coordinate: using FFTW, IndexFunArrays; fft(delta((4,4),offset=(1,1))).
This does make a significant difference for the necessary operations and lends validity to the real-casting in real space. Yet, I think this only works in 1D.
Discussing with @roflmaostc we concluded that the essential problem with even-sized arrays is the following ambiguity: The highest frequency may represent a cosine-shaped wave with phase 0 or any phase shifted wave with a much higher amplitude but a different phase. As opposed to odd-sized arrays: Not only can the FFT highest frequency not distinguish between these cases but also the real-space array itself is ambiguous in this respect. In Fourier space, the lowest (highest negative) frequency is really the sum of two (unknown) phasors with unknown relative phase angle and corresponding magnitude such that the average yields the stated real part.
This means that there is really no unambiguous “correct” way to disentangle this, but one can argue about different ways to do this.
- Assuming the (fundamentally unknown) components to have real-only phase.
The consequence is that on average (randomly) distributed phases, we are overestimating the power of this frequency. Yet the advantage is that centrosymmetric distributions (known to be real in Fourier space) are fine. Splitting is (opposed to what I wrote before) not required to halve the value, but just to copy it upon padding in Fourier space. - We keep real-valuedness but mulitply this frequencies by sqrt(2) upon splitting. This assumes that these frequency did realy have an unkown randomized phase before they were averaged together.
Since we have no chance to correct the (unkown) phase, we can at least try to correct the statistical effect on the magnitude the averaging was having. Advantage: The power density after padding (or even before padding) has ON AVERAGE a better representation. Disadvantage: For centrosymmetric structures we get a too strong powerspectrum. Disadvantage: Upon clipping the array (back) what do we do then? To be reversible in the operation it sounds like we would need to devide be sqrt(2) after averaging the corresponding pixels. But here we actually do know the phases and simple averaging is the best result. - We set the highest frequency to zero before padding. This is a loss of data which probably should be avoided if possible.
Taken together it looks like option one is the way to go.
@stevengj has expertise on this topic.
Yes, the Nyquist component is ambiguous in this sense, but there are ways to disambiguate. For example, for a band-limited signal there is a unique choice that minimizes the mean-square slope of the interpolated function, and the same choice also has the nice property of giving a real interpolant from real inputs. I have some notes on this: https://math.mit.edu/~stevenj/fft-deriv.pdf
I am not sure, I understand. I think this ambiguity is pretty fundamental. Let me give an example:
A cosine curve at the limit with even array size yields maximum contrast, but a sine curve is entirely invisible as it yields only zeros. If we now sample a function with a phase in between the two extreme cases, we would always interprete this as a cosine-only. Maybe it is the best we can do, but it is unfortunately not correct. Of course this is NOT a problem of the FFT but a problem of just about incorrect sampling as we are AT the Nyquist limit rather than sampling above it. @stevengj do you agree?
This is just an example of aliasing. Yes, looking only at the samples n it is impossible to distinguish between 2a\cos(\pi n) = ae^{i\pi n} + ae^{-i\pi n} and 2a\cos(\pi n) + 2b\sin(\pi n) = (a-ib)e^{i\pi n} + (a+ib)e^{-i\pi n} for any b, because the sine term vanishes on all the samples (equivalently, e^{\pm i \pi n} are aliased). You must assume something about what happens in between the samples to determine b. For example, if you want the interpolant that minimizes the mean-square slope you get b=0, since the mean-square slope is:
which is obviously minimized for b=0, whereas a is fixed by the sample observations.
Interpolating inherently involves choices, because there are infinitely many functions that interpolate between a discrete set of samples. Trigonometric interpolation is no different. Fortunately, there are reasonable criteria (like bandlimited interpolants with minimal mean-square slope) that lead to a simple choice of interpolant.
In a sense, yes — if you sample a DFT exactly at the Nyquist limit then the bandlimited assumption is not sufficient to determine the interpolant. You need an additional criterion like the minimal-slope condition above. But I don’t see this as a big problem; after all, the band-limited assumption itself is usually just an approximation.
(For the classic form of the Nyquist–Shannon sampling theorem, you have infinitely many samples — a DTFT, not a DFT — giving a continuous set of frequencies in [-\pi,\pi] and corresponding amplitudes. In this case, the Nyquist frequency \pm \pi is a set of measure zero and doesn’t contribute anything to the interpolated signal as long as the Fourier transform is regular (no delta functions). Hence the bandlimited assumption is sufficient by itself to determine a unique interpolant/reconstruction.)
Thanks for the detailed explanation. As I am in the field of optics, we are in the lucky position that for the ideal (i.e. noise-free) signal the band limit is indeed fulfilled.
On a more practical side though: When it comes to implementing the Fourier-padding as we are attempting in our FourierTools package right now, we did implement it under the assumption of the highest (1D) frequency being real, which corresponds to the case you mention (minimal curvature) as far as I can see. Yet, if I understand it correctly this frequency in a DFFT is NOT just the aliasing, which would be a sum, but it is an average. But if we have very informative signals, they should be close to random and this would probably mean a more or less random phase. And this would mean that “on average” we would be slightly better off to assume this value to be representing a by sqrt(2) reduced version. Agreed?
I am not insisting to implement it like this as many practical signal are symmetric (i.e. real for this frequency). On the other hand, if we have a symmetric signal that is randomly offset by subpixel positions: There you go. This frequency may be on average better represented by the sqrt(2) assumption.
Aliasing means you can’t distinguish e^{\pm i\pi x}, i.e. frequencies \pm \pi, and can only obtain the sum of the two coefficients, which is exactly what is going on here as I explained above (you get 2a = (a-ib) + (a+ib)).
I don’t know what this means. What is an “informative signal”? Are your signals actually periodic, or are you using the DFT as a windowed approximation for the DTFT? Is this measured data that is sampling a physical continuous-time signal that you are trying to reconstruct, or…
RainerHeintzmann:
Yet, if I understand it correctly this frequency in a DFFT is NOT just the aliasing, which would be a sum, but it is an average.
Aliasing means you can’t distinguish e^{\pm i\pi x}, i.e. frequencies \pm \pi, and can only obtain the sum of the two coefficients, which is exactly what is going on here as I explained above (you get 2a = (a-ib) + (a+ib)).
I know what aliasing means. The point I am making is the instead 2|a|, which is what I would have expected at first, the DFFT as implemented yield only |a|. You can check for yourself by Fouriertransforming a delta function on an even sized grid.
But if we have very informative signals, they should be close to random and this would probably mean a more or less random phase. And this would mean that “on average” we would be slightly better off to assume this value to be representing a by sqrt(2) reduced version. Agreed?
I don’t know what this means. What is an “informative signal”? Are your signals actually periodic, or are you using the DFT as a windowed approximation for the DTFT? Is this measured data that is sampling a physical continuous-time signal that you are trying to reconstruct, or…
I agree, “informative” is not really an established scientific term. What I meant is that the signal looks close do noise (== cannot be compressed == informative) as opposed to something that is for example known to be real and symmetric. Let’s assume that the signal is band-limited including only the highest frequency in question but nothing beyond. Maybe I can do a simulation of 100s of noisy signals which are filtered, to demonstrate how you probably get get a better estimate by the sqrt(2) scaling at the Nyquist sampling limit than by the mean (or sum as you mentioned).