By “Nyquist Sample” you mean the highest frequency bin?
Yeah, Julia doesn’t have the symmetric flag. I guess you are referring to rfft instead?
Nyquist Sample only exist in the case the number of samples in Frequency Domain is even. Then it is the n / 2 sample (Before fftshift()) or the 1st (After fftshift()).
If you don’t have symmetric option, then just multiply it by 0.5 and take its real part.
1 Like
I think the using symmetric (rfft) is probably the better choice because by design it does the right thing, doesn’t it?
I was wondering about a complex array. Wouldn’t it be also better to use rfft for imaginary part and real part separately and add them afterwards back together?
Because for a purely real/complex array we know the symmetry property and we know that rfft does the right thing.
Puh, this one never ends.
Using rfft and irfft also induces some problems because the irfft expects that the input is basically hermitian (if you would unfold it to the original size). If you ensure this, you get entirely wrong results in some rows.
In this case the first row is mostly wrong.
julia> z
5×8 Array{Complex{Float64},2}:
0.809667-0.223829im -0.146846+0.483231im 0.126273-0.673676im -0.0995014-0.652765im 1.08065+0.179737im -0.0545377-1.29798im 0.291862+0.335662im 0.244291-0.939019im
0.0796073-0.551711im -0.261923-1.15531im -0.565023-0.699141im -0.374933-1.15981im -0.308387-1.51755im -0.194571-1.77467im -0.912418+0.388777im 0.0659049+0.528615im
-0.49447+0.298023im -0.905227+0.545565im -1.10754+0.15279im -0.776708+0.0774121im 0.647339-0.294794im 0.426044-0.803606im 0.393879-0.603508im -0.653537+0.483806im
0.521922-0.491875im -0.585717+1.07584im 1.85847-0.241757im -1.08196+0.141852im 0.990565-0.123956im 1.42857-0.675697im -0.0593938-0.513508im 0.353353+0.538698im
0.672973+0.311109im -0.134737-0.770579im 0.571894-1.15421im 1.0748-0.245317im -0.226026+0.132937im -0.424684+0.98478im -0.389344-0.941847im 0.866529-0.262219im
julia> rfft(irfft(z, 8))
5×8 Array{Complex{Float64},2}:
0.809667+0.0im 0.0487227+0.711125im 0.209068-0.504669im -0.0770195+0.322609im 1.08065+0.0im -0.0770195-0.322609im 0.209068+0.504669im 0.0487227-0.711125im
0.0796073-0.551711im -0.261923-1.15531im -0.565023-0.699141im -0.374933-1.15981im -0.308387-1.51755im -0.194571-1.77467im -0.912418+0.388777im 0.0659049+0.528615im
-0.49447+0.298023im -0.905227+0.545565im -1.10754+0.15279im -0.776708+0.0774121im 0.647339-0.294794im 0.426044-0.803606im 0.393879-0.603508im -0.653537+0.483806im
0.521922-0.491875im -0.585717+1.07584im 1.85847-0.241757im -1.08196+0.141852im 0.990565-0.123956im 1.42857-0.675697im -0.0593938-0.513508im 0.353353+0.538698im
0.672973+0.0im 0.365896-0.25418im 0.0912752-0.10618im 0.32506-0.615049im -0.226026+0.0im 0.32506+0.615049im 0.0912752+0.10618im 0.365896+0.25418im
So these functions are not the inverse of each other.
It basically ends at the same point, that you somehow need to make them “hermitian”. I already have this one, but for full fft(arr) arrays. For the rfft it seems just be more complex.
Maybe, I finally reached a working state which I’m satisfied with.
Sinc Interpolation with FFT
I’m always using fft for both real and complex signals. For real signals, in the even length signal case, I artificially change certain frequencies on the positive side, to create an array which has the hermitian symmetry. By doing so, the result from iffting is purely real.
Sinc Downsampling
I’m always using fft for both real and complex signals. In downsampling where my output size is even in an dimension, I artificially change the highest negative frequencies (where the positives partners are missing) to obey the symmetries so that iffting returns a purely real signal.
All in all, upsampling and downsampling are approximately inverse to each other and recover the initial signal pretty well.
I also registered that package and it’ll be merged in the next days
I still don’t understand why you bother. Just ignore the symmetries (they are bothersome to deal with, esp in higher dimensions), and if the input to the function was real, just discard the imaginary part after up/down sampling.
Hm, didn’t feel right to do so.
With the current routine, upsamling and downsampling are better inverse to each other in comparison to the routine where we take always the real part after iffting.
Nevertheless , thanks for your valuable input, I definitely learnt something!
If the resulting spurious imaginary parts are large, they probably indicate a bug in your algorithm. It’s pretty easy to get FFT-based interpolation (or differentiation etc.) wrong.
norm(imag(output)) / norm(real(output)) should be nearly on the order of machine precision.
2 Likes
Not in this very mparticular case (downsampling from odd to even size). Of course the resulting imaginary part should be small (or your downsampling will be bad anyway and you should have filtered it before) but not of order of the machine precision.
This is the wrong policy in my opinion.
But the code I supplied solves the case for downsmapling.
I really don’t understand what’s missing.
Not to say you can always do the downsmapling by upsampling and decimation.
@RoyiAvital It does maybe in 1D, but as far as I tested this, it doesn’t generalize to N dimensions. After iffting, the result contains a non-zero imaginary part. Using a symmetric FFT (as you did) provides real results, but at least with FFTW (via rfft) the real results are messed up a lot (low and high frequencies in the first dimension).
And downsampling by decimation introduces aliasing because you simply undersample a high frequency signal. Taking a frequency window from a FFT spectrum is not the same as decimation because it is alias free.
(by decimation you mean taking every N-th sample?)
What I mean:
rfft requires a special symmetry to produce a correct real result. If you don’t provide this symmetry, you really mess up low and high frequencies:
The symmetry can be seen here:
julia> rfft(randn((6,6)))
4×6 Array{Complex{Float64},2}:
4.36498+0.0im 5.59917+4.04078im 2.42256+1.53944im 0.583178+0.0im 2.42256-1.53944im 5.59917-4.04078im
0.658702+6.55594im 1.23211+10.0951im 3.90816+2.07009im -2.79509+2.73057im 6.50781-2.18167im -2.71237-4.55771im
-0.64458-4.56344im 4.39054-8.0552im -4.49956-3.66493im -1.37505-2.99915im -6.16341-6.68374im -2.62665-4.84961im
5.83408+0.0im -1.513+5.81849im -9.69012-0.0557974im -6.16479+0.0im -9.69012+0.0557974im -1.513-5.81849im
(1,1), (1,4), and (4,1) are purely real. Furthermore (2, : ) is has the hermitian symmetry.
However, if you would naively cut the highest frequency of a 7x7 array to get a 6x6 array, you don’t have the symmetry:
julia> rfft(randn((7,7)))
4×7 Array{Complex{Float64},2}:
3.66183+0.0im -1.02749+1.94808im 1.51224-14.965im -1.15603+6.66689im -1.15603-6.66689im 1.51224+14.965im -1.02749-1.94808im
-1.03278+4.97944im -0.93996+8.06737im -6.88103+3.68525im 8.17271-4.17101im 4.03284+5.643im -0.919963-2.10859im -6.10537-1.09231im
-1.03219+5.51279im 6.46411-11.5258im -0.71596-6.65423im -1.37923-10.1864im 0.996462+1.32339im -2.89172-1.68502im 4.21725-5.09607im
0.34954+4.8042im 4.59035-4.21084im 0.637454-1.86879im 11.3116+6.97052im 5.91684+7.88591im 7.2835-2.23067im 1.83935+3.89034im
You are not providing the right symmetry rfft expects.
Let’s look at this, to downsample an odd sized 2D array to an even sized
function ds_2D(arr)
a = rfft(arr)
a = a[:, 1:end-1]
a[:, end] .*= 0.5
a[end, :] .*= 0.5
a[end, end] *= 2
return irfft(a, size(arr)[1]-1)
end
ulia> rfft(x)
4×7 Array{Complex{Float64},2}:
5.19375+0.0im -1.22145-3.05277im 1.99423+3.27891im -5.63567+2.6869im -5.63567-2.6869im 1.99423-3.27891im -1.22145+3.05277im
0.468664+2.7739im -4.30717-2.67216im -1.11728+2.58736im -2.3358+3.52392im 0.419257-5.8995im -12.4438-5.26935im 0.413576+2.83563im
-5.94199+4.45657im 4.97952+2.69954im -5.49123+7.09337im -1.05639+5.62915im -11.0779-6.8679im 7.01256+5.8489im 4.35303-3.09863im
-4.61229-6.20257im -0.437476-8.4717im -3.10097-2.61191im 3.93253+5.40222im 1.05266-3.04878im 0.00995383-6.20056im 3.80061-8.75746im
julia> ds_2D(x) |> rfft
4×6 Array{Complex{Float64},2}:
5.19375+0.0im -0.112167-0.706655im -1.82072+2.98291im -5.63567+0.0im -1.82072-2.98291im -0.112167+0.706655im
0.468664+2.7739im -4.30717-2.67216im -1.11728+2.58736im -2.3358+3.52392im 0.419257-5.8995im -6.2219-2.63467im
-5.94199+4.45657im 4.97952+2.69954im -5.49123+7.09337im -1.05639+5.62915im -11.0779-6.8679im 3.50628+2.92445im
-2.30615+0.0im -0.106881-0.567786im -0.512076+0.109217im 1.96626+0.0im -0.512076-0.109217im -0.106881+0.567786im
You would think that the ds_2D(x) |> rfft is just a subset (maybe the highest frequency factor 0.5) of rfft(x) because we cut it in that way. But the first frequency line is completely different! And this will affect the 2D array because the artifacts are distributed over a whole frequency range in the first dimension.
Why do you insist on rfft()? Just do regular fft().
As I wrote, for real signal add to my code real() to the Nyquist Sample.
Another way is do downsample by upsample.
If your upsample works and validated, just do the trick I posted in Sinc Interpolation based on FFT - #27 by RoyiAvital.
Pay attention that while no LPF was applied this is the correct interpolation in the meaning those are the values of the “continuous” signal at those points.
To apply Downsamplign in the meaning of estimating the value as the signal was sampled at lower rate you need to follow what I did on the DFT.
I would still say you are doing it wrong if you are doing FFT-based downsampling and your algorithm doesn’t produce exactly conjugate-symmetric data in the frequency domain. Probably you have unbalanced handling of the Nyquist element(s).
By taking the real part of the output, you are throwing away the error you made, but as long as this error was in a small number of frequency components then it doesn’t affect the answer by much.
2 Likes
My point is that taking the real part of the output is exactly equivalent to a proper handling of the nyquist elements, while being much simpler to write Anda lot less error prone. (Of course the proper way of doing things is exploiting the real FFTs, but I can never be bothered.)
Yeah I agree on that.
For upsampling it’s reasonable, I agree on that.
But for downsampling: I still don’t understand what’s the proper way? In 1D it is clear (multiplication of the Nyquist sample by 0.5). But for 2 or N dimensions, it isn’t that simple as I pointed out above, the symmetry handling is tricky.
Am I wrong or did you manage to do the Upsampling for Multi Dimensional Arrays?
If so, the logic of handling the Nyquist Sample for downsampling is exactly the same.
I did, yes.
How can it be the same?
For upsampling you append slices with hermitian symmetry to the other boundary, right?
For downsampling you naively would take the frequencies you want to keep, but then you need to fix the symmetry in the “Nyquist slices”?
The rationale of the adjustment for upsampling / downsampling is the same. The only difference is the indices.
Sorry, I don’t really get what you mean by this 