plz help.... ifft of 1 / (1 - fft(g(t)) )
2 次查看(过去 30 天)
显示 更早的评论
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
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?
回答(10 个)
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
0 个评论
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:
- Scale by dt for the FFT, and by Fs for the IFFT
- Scale by 1/M for the FFT, and by M for the IFFT
- Scale by 1 for the FFT, and by 1 for the IFFT
- Scale by 1/sqrt(M) for the FFT, and by sqrt(M) for the IFFT.
- 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 个评论
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
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
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
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
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
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
0 个评论
J J
2011-9-14
6 个评论
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.
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 个评论
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.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!