Compare ifft in Matlab with Julia for complex input

Specific Domains Signal and Image Processing
1

My use case is frequency demodulation using the method of Coats and Randall, with the link being one of many from a web search. They use the Hilbert transform for this purpose.

When I use ifft in MATLAB on a complex vector, the resulting vector is ~2x the length of the input vector.

In the MATLAB case the complex vector does not have the symmetry of having the positive and negative frequencies present. So what I think is happening is that MATLAB adds in the negative frequencies and then performs the ifft from the call to the FFTW library. Using this hypothesis I wrote the following code to do what I think MATLAB is doing:

# sf is the complex frequency vector obtained via rfft
# In this case it is assumed that the amplitude has been normalized by multiplying by 2 (except the DC component and the Nyquist frequency if the length is even.
nsf = length(sf)
if iseven(nsf) 
    sfd = Vector{ComplexF64}(undef, 2nsf-2)
    sfd[1] = sf[1]
    sfd[2:nsf-1] = sf[2:nsf-1] ./ 2
    sdf[nsf] = sf[nsf]
    sfd[nsf+1:end] = sfd[nsf-1:-1:2]
else
    sfd = Vector{ComplexF64}(undef, 2nsf-1)
    sfd[1] = sf[1]
    sfd[2:nsf] = sf[2:nsf] ./ 2
    sfd[nsf+1:end] = sfd[nsf:-1:2]
end
    
sc = ifft(sfd).*n                              # complex time from hilbert transform

Doing this I get a result that is different from Matlab. I am wondering if anyone has insight into this.

1 Like
2

Just guessing: If your input is computed with rfft, you will likely want to use irfft to compute the inverse. API · AbstractFFTs.jl

3

I don’t think Matlab’s and Julia’s fft/ifft behave differently, and I’m skeptical that the following is correct:

It seems to contradict the documentation of ifft here: Inverse fast Fourier transform - MATLAB ifft - MathWorks Nordic

Can you provide a minimal working example to show the difference you are seeing?

(And, as @fatteneder says, rfft is not the same as fft, that is probably a crucial point here.)

1 Like
4

Thanks for the challenge of making a MWE. It turns out that some of my thinking was wrong. This is what I should have been doing:

  • Transform to frequency domain
  • Set negative frequency components to zero
  • Multiply positive frequency components by 2 (do not change zero frequency component and take care with the Nyquist frequency component)
  • Transform back to time domain
  • Take real part of the analytic signal

Ergo, you have the original signal.

As an example

using FFTW

# Create signal
dt = 0.2
t = 0:dt:1
y = sin.(2*pi*t)
df = 1/t[end]
yf = fft(y)

# Make signal one sided like rfft
n = length(t)
np = div(n, 2) + 1
yf[np+1:end] = zeros(ComplexF64, n-np)

# Put energy lost from zeroing -ve frequency back in 
if iseven(n)
    yf[2:np-1] = 2*yf[2:np-1]
else
    yf[2:np] = 2*yf[2:np]
end

# Inverse FFT
ynew = ifft(yf)
[y ynew]   # compare y with real part of ynew

MATLAB does the same thing

1 Like