Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Reverse FFT of a ratio of two FFTs

  1. Dec 14, 2015 #1
    Dear all,
    I don't know if this is the correct place for this question, but I did a little search on the forum and I saw that most FFT related questions have been posted here.
    My question is this: I need to deconvolve two real signals (in my case they are two probability density functions), so I go in Fourier space, I divide them and then I go back to real space.
    Now, before apply FFT I need to sample the signals, i.e. I have to define a sampling interval and frequency. Now, if the sampling interval is on the positive axis (the signal s is sampled at discrete time points s[t], t>=0) then everything is fine, but if the sampling interval includes negative points, then when I perform the inverse FFT on the FFTs of the two signals, the final signal results swapped, i.e. the
    signal corresponding to the negative sampling points is put to the right, after the positive sampling points.
    Can you please help me understand why this happens? Notice that if I perform iFFT(FFT(s)), where iFFT is the inverse FFT and s is my sampled signal, everything is fine even when the sampling interval comprises negative time points (i.e. s[t], t<0).

    Here you can see an example of the code (R language):

    Code (Text):
    #Define sampling interval
    fs <- 0.05
    n0 <- -30
    nMax <- 100
    N=(nMax-n0)/fs
    mids <- seq(0, (N-1))*fs+n0

    #Sample two continuous probability density functions (fO and fB) that must be deconvolved
    fO = density(rnorm(1000000, 40, 6), bw=1.5, from=n0, to=nMax, n=N)
    fB = hist(rnorm(1000000, 14, 2), breaks=c((mids-fs/2), (tail(mids,1)+fs/2)), plot=F)$density

    #Perform the FFT
    fO.fft <- fft(fO$y)
    fB.fft <- fft(fB)

    #Divide the FFTs (deconvolution) and go back to real space
    fT <- abs(fft(fO.fft/(fs*fB.fft), inverse=T))/length(fO$x)

    #Plot the results to check the deconvolution
    plot(fO$x, fB, type='l', col=1)  #First signal
    lines(fO, type='l', col=2) #Second signal
    shift <- which(mids==min(abs(mids)))  #Now I need to swap the deconvolved signal... why????
    fT.shifted <- c(tail(fT, shift), fT[shift : length(fT)-shift]) #Swapped deconvolved signal
    lines(mids, abs(fT.shifted), type='l', col=3)
    lines(mids, dnorm(mids, 26, sqrt(36-4)), col=4) #Theoretical deconvolved signal
    legend(x='topright', legend=c('B', 'O', 'Estimated T', 'True T'), col=1:4, lwd=1)

     
     
  2. jcsd
  3. Dec 14, 2015 #2

    mfb

    User Avatar
    2016 Award

    Staff: Mentor

    I don't know R, but it probably treats the distributions as periodic (at least the underlying mathematics does that), having an entry at position x or at x+N is the same for the algorithm.

    I moved the thread to the programming section.
     
  4. Dec 14, 2015 #3

    Ibix

    User Avatar
    Science Advisor

    Note: I've been taught multiple ways to define the Fourier transform, with factors of ##2\pi## in different places. This post uses ##F(u)=\int f(x)e^{-i2\pi ux}dx##. Depending on your convention, you may find stray factors of ##2\pi## popping up below.

    Your problem is a consequence of the combination of your FFT algorithm assuming that the first point of your input array lies at x=0 and the implicit periodicity that mfb mentioned.

    It's fairly easy to show that, if the Fourier transform of ##f(x)## is ##F(u)##, then the Fourier transform of ##f(x+a)## is ##e^{i2\pi au}F(u)##. Trivially, then, the inverse Fourier transform of ##e^{i2\pi au}F(u)## returns ##f(x+a)##.

    If you take two functions, ##f(x+a)## and ##g(x+a)##, Fourier transform them and divide them, the complex exponential cancels out leaving just ##F(u)/G(u)##. When you inverse Fourier transform you get back ##f(x)## with ##g(x)## deconvolved out - in other words, the answer you want offset by ##a## compared to what you were expecting. So actually your negative time points have dropped off the beginning of your array. Fortunately, discrete Fourier transforms assume that their input is periodic - so the negative time points of the next iteration have appeared at the end of your array (if that makes sense).

    Assuming you had ##n## negative time points, all you need to do is chop off the last ##n## points (or is it the last ##n+1##? Not sure where the zero goes any more - experiment with something you know the answer to!) of the final array and stick them on again at the beginning as, I think, you have realised. This was usually called a block-swapping operation back when I last did it (10+ years ago, so beware of language drift).
     
  5. Dec 14, 2015 #4

    mfb

    User Avatar
    2016 Award

    Staff: Mentor

    Alternatively, multiply the fft ratio by the factor ##e^{i 2 \pi a u}## (up to factors of ##\pm 2 \pi##).
     
  6. Dec 15, 2015 #5
    Thank you all for your helpful replies!
     
  7. Dec 15, 2015 #6

    Daz

    User Avatar
    Gold Member

    Just a practical comment: deconvolution using the FFT is prone to numerical instability if the input data contains noise. Solutions that look good on paper sometimes don’t work terribly well in a real-life digital computer.

    There is a fantastic book that I’d highly recommend: Numerical Recipes - The Art of Scientific Computing by Press, Teukolsky, Vettering and Flannery.

    Quoting directly from the section on deconvolution: “The process is generally quite sensitive to noise in the input data and the accuracy to which the response function is known. Perfectly reasonable attempts at deconvolution can sometimes produce nonsense for these reasons.” (Page 542, Second Edition, C-language version.)

    The book goes on to discuss solutions to these problems and also presents alternative algorithms.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Reverse FFT of a ratio of two FFTs
  1. FFT library (Replies: 2)

  2. Getting FFTs to behave (Replies: 0)

  3. FFT code help please (Replies: 16)

  4. [C#] FFT in read world (Replies: 33)

Loading...