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

4 次查看（过去 30 天）
J J
J J 2011-9-10
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?
  2 个评论
Rick Rosson
Rick Rosson 2011-9-10
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
J J 2011-9-10
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（10 个）

Rick Rosson
Rick Rosson 2011-9-12
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
Rick Rosson
Rick Rosson 2011-9-13
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
  2 个评论
J J
J J 2011-9-13
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
Rick Rosson 2011-9-13
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.

请先登录，再进行评论。

Wayne King
Wayne King 2011-9-10
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
  1 个评论
J J
J J 2011-9-10
well i dont get it. and btw, why did you drop dt in you recirc = ifft(...) ? dt~=1 , it is 0.2

请先登录，再进行评论。

Rick Rosson
Rick Rosson 2011-9-12
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
  1 个评论
Walter Roberson
Walter Roberson 2011-9-12
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".

请先登录，再进行评论。

Wayne King
Wayne King 2011-9-12
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
  2 个评论
J J
J J 2011-9-12
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
J J 2011-9-12
and we need to include df in our ifft term as well.

请先登录，再进行评论。

J J
J J 2011-9-13
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).
Rick Rosson
Rick Rosson 2011-9-14
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
J J
J J 2011-9-14
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????
  6 个评论
显示 4更早的评论
Rick Rosson
Rick Rosson 2011-9-15
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
J J 2011-9-15
yes you are correct - i was not using the standard notation f instead of omega. i corrected it now.

请先登录，再进行评论。

Rick Rosson
Rick Rosson 2011-9-15
Hi J J,
I don't think I can help you. Maybe someone else can.
Best of luck,
Rick
  1 个评论
J J
J J 2011-9-15
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.

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2011-9-15
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.
  2 个评论
J J
J J 2011-9-15
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
Walter Roberson 2011-9-15
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by