plz help.... ifft of 1 / (1 - fft(g(t)) )

I do not quite understand how fft and ifft work in MATLAB. I have two normalized functions g1 and g2. I need to calculate ifft(1 / (1 - fft(g1).fft(g2) ) ). I don't get what I should get, so I know I am missing something in my code.

MATLAB code

alpha2 = 8; beta2 = .8;
alpha1 = 3; beta1 = 2;
td = 8;
dt = .2; t = [0:dt:359]; T = length(t);
g1 = (t>t(1)) .* ((t-t(1)) .^alpha1) .* exp(-(t-t(1)) /beta1) / ...
  ( (beta1 ^(alpha1+1)) * gamma(alpha1+1) );
g2 = (t>td) .* ((t-td).^alpha2) .* exp(-(t-td)/beta2) / ...
  ( (beta2 ^(alpha2+1)) * gamma(alpha2+1) );
ft_g1 = fft(g1);
ft_g2 = fft(g2);  
ft_recirc = ones(size(ft_g1))./(ones(size(ft_g1))-ft_g1.*ft_g2);
recirc = dt.*ifft(ft_recirc);

What am I missing?

[Rick Rosson]

When you say "I don't get what I should get...", what specifically were you expecting to get, and how are your actual results different from your expected results?

[J J]

Hi Rick, if I use maclaurin series expansion for x = ft_g1.*ft_g2 , I get:

1 / (1-ft_g1.ft_g2) = 1 + ft_g1.ft_g2 + (ft_g1.ft_g2)^2 + (ft_g1.ft_g2)^3 + (ft_g1.ft_g2)^4 +...

taking ifft of this expansion, I get what I expect biologically speaking. This is a biological problem where I know the answer. I am just struggling with ifft of the division above. I want to avoid using the maclaurin expansion. thx.

Answer 1

[Rick Rosson]

Hi J J,

I don't know what the issue is just yet, but I have a few suspicions that I am trying to verify.

Meanwhile, I have modified your code a bit to make it easier to read and added some plots to help visualize what is going on. Maybe you or someone else can figure it out from this code:

   alpha2 = 8; beta2 = .8;
   alpha1 = 3; beta1 = 2;

   td = 8;

   dt = 0.2; 
   t = (0:dt:359)'; 
   M = size(t,1);

   Fs = 1/dt;
   dF = Fs/M;
   f = -Fs/2:dF:Fs/2-dF;

   g1 = (t>t(1)) .* ((t-t(1)) .^alpha1) .* exp(-(t-t(1)) /beta1) /( (beta1 ^(alpha1+1)) * gamma(alpha1+1) );
   g2 = (t>td)  .* ((t-td).^alpha2) .* exp(-(t-td)/beta2) / ( (beta2 ^(alpha2+1))* gamma(alpha2+1));

   G1 = dt*fftshift(fft(g1));
   G2 = dt*fftshift(fft(g2));

   G = G1.*G2;

   H = 1./(1-G);

   h = Fs*ifft(ifftshift(H));

   figure;
   plot(t,g1,t,g2);

   figure;
   plot(f,abs(G1),f,abs(G2),f,abs(G));
   grid on;

   figure;
   plot(f,abs(H));
   xlim([-0.5,0.5]);
   ylim([0.5,2.0]);
   grid on;

   figure;
   plot(t,h);
   xlim([0,95]);
   ylim([-2.00924,-2.00922]*1e4);

HTH.

Best,

Rick

Answer 2

[Rick Rosson]

In most situations, scaling is really not all that important. The overall shape of the spectrum matters much more than the absolute scale.

But if you really are worried about it, there are several different conventions from which you can choose:

1. Scale by dt for the FFT, and by Fs for the IFFT
2. Scale by 1/M for the FFT, and by M for the IFFT
3. Scale by 1 for the FFT, and by 1 for the IFFT
4. Scale by 1/sqrt(M) for the FFT, and by sqrt(M) for the IFFT.
5. and so on.

I generally use either option #1 or option #2 depending on my mood and whether it's raining outside.

All of these conventions have one thing in common: The product of the two scaling factors is always 1. Please note that the ifft function in MATLAB includes a scaling factor of 1/M as part of the computation, so that the overall round-trip scaling is 1/M.

HTH.

Rick

[J J]

Rick, today is a sunny day - so I ca use option 1! ;-p

but shouldnt we use Fs instead of dF, as matlab's ifft has that 1/length(t) ?

I will have more questions shortly. Thanks for your responses.

[Rick Rosson]

Yes, that is correct. I did not realize that IFFT includes the 1/M scaling, but I checked the documentation, and it does.

I fixed my answers to account for this factor.

Answer 3

[Wayne King]

Hi, If you have formulated what you want to find the inverse Fourier transform of correctly, then:

ft_g1 = fft(g1);
ft_g2 = fft(g2);
recirc = ifft(1./(1-ft_g1.*ft_g2));

Wayne

[J J]

well i dont get it. and btw, why did you drop dt in you recirc = ifft(...) ? dt~=1 , it is 0.2

Answer 4

[Rick Rosson]

Hi J J,

Can you explain a bit of what this particular math operation is supposed to accomplish? I see that you are taking the FFT of two signals, combining them as a more complex transfer function, and then taking the IFFT back to the time domain. The signals are exponential pulses with various time delays and coefficients. What is the goal of this technique?

Thanks.

Rick

[Walter Roberson]

I was thinking that perhaps there was some theorem about this transformation that I hadn't happened to come across, so I experimented with some functions. With the first pair of functions I choose, after bashing it with a bunch of mathematical transformations, I got two suggestive solutions: one of them a piecewise that was undefined for w=1 or w=-1 and otherwise was Dirac(t). The second suggestive solution did was plain Dirac(t).

Thinking that maybe I had found a generalizable solution, I tried a different set of functions, simple functions. After a lot of hammering, the solutions I got were all pretty obtuse, but two of the solutions clearly simplified to "undefined".

Thus, if there is some general theorem that is being explored here, either it is too sophisticated for the kind of bashing I did, or else the theorem is that "the result is undefined".

Answer 5

[Wayne King]

I have the same question that Rick has in that I'm not sure what you are trying to accomplish with that operation, but you do not need to scale the inverse DFT by dt. In fact, if you were to scale the inverse DFT by dt, you would not get perfect inversion.

For example:

        % Sampling rate is 1 kHz
        t = 0:0.001:1-0.001;    
        x = cos(2*pi*100*t);
        xdft = fft(x);
        % Invert the DFT
        y = ifft(xdft);
        % Get the Linf norm of the difference
        norm(x-y,Inf)
        % Get the L2 norm of the difference
        norm(x-y,2)

Now repeat the above except scale the inverse DFT by 0.001. You'll see that is not what you want to do.

Wayne

[J J]

I dont understand - we need dt to scale the amplitude of fft( x(t) ) equal to its analytical Fourier transform x(f). dont we - cause analytically we have integral wrt to dt? i was wrong in including dt in my ifft - as Rick pionted out I should have included it in my fft term. am i correct?

[J J]

and we need to include df in our ifft term as well.

Answer 6

[J J]

ok another question: implement the code:

       alpha2 = 8; beta2 = .8;
       alpha1 = 3; beta1 = 2;
       td = 8;
       dt = 0.2; 
       t = (0:dt:359)'; 
       M = size(t,1);
       Fs = 1/dt;
       dF = Fs/M;
       f = -Fs/2:dF:Fs/2-dF;
       g1 = (t>t(1)) .* ((t-t(1)) .^alpha1) .* exp(-(t-t(1)) /beta1) /( (beta1 ^(alpha1+1)) * gamma(alpha1+1) );
       g2 = (t>td)  .* ((t-td).^alpha2) .* exp(-(t-td)/beta2) / ( (beta2 ^(alpha2+1))* gamma(alpha2+1));
       G1 = dt*fftshift(fft(g1));
       G2 = dt*fftshift(fft(g2));
       G = G1.*G2;
       H = 1./(1-G);  
       h = Fs*ifft(ifftshift(H));
       plot(t,h); grid on;

now,

    h_corrected = h; 
    h_corrected(1)=h_corrected(2); 
    plot(t,h_corrected); grid on;     % this is close to the one i was looking for!  this is closer to the biologically correct answer.  but why do we have that mistake at h(1)??  and still -  why do we have that dip h(2:td)?  you see that dip?  it should not be there.

Also, I dont understand why the amplitude of h is so large O(10^4).

Answer 7

[Rick Rosson]

Hi J J,

I still do not understand what this operation represents either mathematically or biologically. Can you please give us a brief explanation of what this code is supposed to do and why you are trying to do it? Alternatively, can you provide a link to a web-based reference that explains what it is?

Thanks!

Rick

Answer 8

[J J]

ok. before I get into dirty details of my application, lets assume a much simpler case - the g=exp(-t) case:

Theory:

g(t) = exp(-t) -> G(f) = fft(g(t)) = 1/(1+2i.pi.f)

H(f) = 1/[1-G(f)] = (1+2i.pi.f) / (2i.pi.f)

h(t) = ifft(H(f)) = ifft(1/(2i.pi.f)) + ifft(1) = sgn(t)/2 + dirac(t)

THUS,

plot(t,h(t)) = plot( t , sgn(t)/2+dirac(t) )

Matlab:

dt = .04;
t = [0:dt:31];
T = length(t);
fs = 1/dt;
df = fs/T;
f = -fs/2:df:fs/2-df; 
g = exp(-t);
G = dt.*fftshift(fft(g));
H = 1./(1-G);
h = abs(fs.*ifft(ifftshift(H)));
figure; plot(t,h); grid on

Now, why dont we get what we expect????

[J J]

why do we have a sloped line?

Why is it much larger than 0.5?

why does it have imaginary parts with a large magnitude ? the result should be real, amplitude at .5, no slope.

if we decrease dt, things work better - but still.

[Walter Roberson]

It appears to me that you have used Fourier Transforms for continuous functions; I am not certain that they apply in all the details to discrete fourier transforms.

[Rick Rosson]

Walter is correct. The theory that you are showing is for the Continuous Time Fourier Transform (CTFT), whereas the FFT is an implementation of the Discrete Fourier Transform (DFT). There is absolutely a relationship between the CTFT and the DFT, but they are by no means the same thing, and you should not expect the DFT to give you the same result as the CTFT Theory predicts.

[Rick Rosson]

Furthermore, you are comparing the CTFT as a function of angular frequency (symbol omega in units of radians per second) whereas the MATLAB code is plotting the FFT versus frequency (symbol f in units of cylces per second or hertz). Frequency and angular frequency are not the same thing.

[Rick Rosson]

J J,

Please tell us what you are trying to do, what you expect to see, and why your results are different from expectation. I am not a mind reader. I have been trying to help you for 4 days, but at this point I feel like I am just shooting in the dark. If you would have provided more information 4 days ago (which I asked for), we might have been able to answer your question much more quickly instead of wasting all of this time.

[J J]

yes you are correct - i was not using the standard notation f instead of omega. i corrected it now.

Answer 9

[Rick Rosson]

Hi J J,

I don't think I can help you. Maybe someone else can.

Best of luck,

Rick

[J J]

Rick, you have helped me a great deal - at least you corrected me for dt, fftshift, iffshift, and f, df and fs! no time was wasted - at least not for me. further, I found out that in the simple exponential case, the discrepancies are exactly due to the discrete nature of fft (as dt decreases we get closer and closer to the continuous FT result). And I observed the same pattern of discrepancy (the sloped line and larger magnitude) in my lil more complex gamma distribution case. I also realized that the spike at t=0 is correct - according to theory it corresponds to dirac(t). and yes , so you and walter are correct.

but i was wondering whether we know the exact nature of this discrepancy, the mechanisms governing it, and then whether this mechanism can be formulated and then discrepancies somehow corrected mathematically - maybe i could somehow correct it in my code. you know?

and explaining my problem - it is very long and hairy. basically, i am modeling a biological process: injection of a an agent into the blood pool, its different patterns of circulation, absorption, elimination, diffusion into various tissue types, retention in some tissues,.... I could send you my paper, but i think it d be a waste of time for u - as i said it is hairy. specially now that we've seen the same pattern of discrepancy in the simple exponential case, its results can be expanded to my gamma variate case.

Answer 10

[Walter Roberson]

J J, you wrote in your theory section that

g(t) = exp(-t) -> G(w) = fft(g(t)) = 1/(1+2i.pi.w)

As mentioned above, you have attempted to use continuous fourier transform theory in your equations. Unfortunately, the transform you have used is not correct. It is not possible to take the continuous fourier transform of exp(-t) . This is explained in more detail in Peter Conlon's posting near the bottom over here

The result you show, 1/(1+2i*pi*w), is the continuous fourier transform for a different function,

1/(2*Pi) * exp(-(1/2)*t/Pi)*Heaviside(t)

and for comparison, the continuous fourier transform of the cleaner exp(-t)*Heaviside(t) goes to 1/(1+1i*w) (without the 2 or the Pi in the denominator)

As your initial transform is incorrect and cannot be corrected, the remainder of your calculation fails.

If you tested experimentally with fft (the discrete one), then you could have been mislead about the transform of exp(-t) if you had happened to only test for nonnegative t.

[J J]

Hi walter, my t does NOT go from -inf to +inf; rather it starts at 0, thus my integration starts from 0 - that does take care of the heaviside function.

and yes, that integration is correct with the 2 and with the pi according to the conventional formula that i use for my continuous fourier calculations. I think cause earlier i used omega instead of f (my bad - but i meant f), you used the conventional fourier integration for omega - and yeah there is no pi and no 2! but i 'd corrected my post and had changed omega to f. I should stick to standard notation. so it s all fine.

and thx for your time.

i just gotta find a way to formulate the error as a function of dt and Fs now, because this error gets larger and larger in my long long experimental time scale - days!

[Walter Roberson]

fourier transforms are *defined* in terms of integrating from -infinity to +infinity . You cannot just choose to start from 0. You can multiply by the Heaviside function, but you need to include it in your theory.

(IMHO) you have done too much retcon in this topic: you already know what result you want and change the question as you go along when your result becomes untenable based on the question.