Inverse Fourier transforms
in [Matlab]
|
From: Walter Roberson on 19 Jun 2010 11:53 Mat Hunt wrote: >> k = -100:.1:100; >> t = 16.2668*exp(-1/4*k.^2)/(144+2*k.^4-12*k.^2+9*abs(k)); >> ti = ifft([0 t]); >> plot(ti) > I used the above code and didn't get anything remotely sensible. The > code I use is: > > h=1; > E_b=0.1; %0.2; > B=0.2; > rho=1; > F=0.2; > mu=0; > %Start the integration > x=-8:0.01:8; > n=length(x); > y=zeros(n,1); > dk=0.001; > k=-20:dk:20; > l=length(k)-1; > for j=1:n > z=f_1(k,F,B,E_b,h,rho,mu).*exp(i*k*x(j)); > y(j)=real(trapz(k,z)); > end > > plot(x,y) > xlabel('x'); > ylabel('\eta_{1}(x)'); > > Along with the function: > > function eta=f_1(k,F,B,E_b,h,rho,mu) > m=length(k); > for j=1:m > eta(j)=-(sqrt(pi)*exp(-0.25*k(j)^2)/(2*rho*9.8065))./(1-F+(h*k(j))^4/90+0.5*(B-1/3)*(h*k(j))^2-0.5*E_b*h*abs(k(j))+i*mu); > > end > If I don't get anything like the above function then the answer is > incorrect. The code that you gave me didn't give anything like the > above M file. I missed copying in a "." t = 16.2668*exp(-1/4*k.^2)/(144+2*k.^4-12*k.^2+9*abs(k)); should have been t = 16.2668*exp(-1/4*k.^2)./(144+2*k.^4-12*k.^2+9*abs(k)); Anyhow, if you take the fft of your y, you will see that the result is complex whereas the function you provided earlier is not complex and thus cannot be the discrete fourier transform of your y. Is your f_1 intended to be a more detailed version of the function you gave earlier? If so then note that your f_1 function has an imaginary component whereas the function you gave earlier does not. My answer was not incorrect for the question you posted earlier, and I don't appreciate being expected to read your mind about what you are really trying to do.
From: Greg Heath on 19 Jun 2010 13:10 On Jun 17, 6:29 pm, "Mat Hunt" <hunt_...(a)hotmail.com> wrote: > Walter Roberson <rober...(a)hushmail.com> wrote in message <hve5un$mj...(a)canopus.cc.umanitoba.ca>... > > Mat Hunt wrote: > > > I want to compute an inverse Fourier transform and plot the solution. I > > > have the analytical expression of the inverse Fourier transform. I can > > > either supply it as an analytical function (and use ifourier) or give it > > > as a numerical array and use ifft. > > > > Which would be best? > > > I would try ifourier first as it would be more accurate -- using the discrete > > fourier is going to give you round-off problems. > > > If your function is sufficiently complex then there might not be a nice > > symbolic inverse transform; in that case you would have to resort to ifft. > > The function in question is: > > 16.2668*exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k)) Are you sure about abs(:) ?. The function is not analytic. Hope this helps. Greg
From: Greg Heath on 19 Jun 2010 15:12 On Jun 18, 6:04 am, "Mat Hunt" <hunt_...(a)hotmail.com> wrote: > No, there are no problems with it, no complex numbers. The integrand is an even function and real, as a result, the complex part of the solution is zero. No, the function has an imaginary part. I'll include a demo in my next post. Greg
From: Greg Heath on 19 Jun 2010 15:59 On Jun 18, 9:37 am, Walter Roberson <rober...(a)hushmail.com> wrote: > Mat Hunt wrote: > > The function in question is: > > > 16.2668*exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k)) > > > I tried ifourier and I got the result: > > > 16.2668*ifourier(exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k)),k,x) > > > So I think I will have to use ifft. I would ideally like to plot the > > original function but I don't know how to plot it. > > I get something vaguely plausible if I use > > k = -100:.1:100; > t = 16.2668*exp(-1/4*k.^2)/(144+2*k.^4-12*k.^2+9*abs(k)); > ti = ifft([0 t]); > plot(ti) I assume the zero is there to deal with the wild oscillations that occur in the real and imaginary parts otherwise. I have no idea why that happens. It also happens when abs(k) is replaced by k or -k. The following works well: close all, clear all, clc, p = 0 kb = -100:.1:100; % bipolar k dk = kb(2)-kb(1) % 0.1 N = length(kb) % 2001; N is odd k = dk*(0:N-1); % unipolar k Yb = 16.2668*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2+9*abs(kb)); % Yb = 16.2668*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2+9*kb); % Yb = 16.2668*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2-9*kb); Y = fftshift(Yb); p = p+1, figure(p) % figure 1 subplot(211) plot(kb,Yb) subplot(212) plot(k,Y) dx = 1/(N*dk) x = dx*(0:N-1); xb = x-dx*(N-1)/2; y = ifft(Y); ry = real(y); iy = imag(y); ay = abs(y); py = angle(y); p = p+1, figure(p) % figure 2 subplot(221) plot(x,ry) subplot(222) plot(x,iy) subplot(223) plot(x,ay) subplot(224) plot(x,py) % yb = ifft(Yb); yb = ifftshift(y); ryb = real(yb); iyb = imag(yb); ayb = abs(yb); pyb = angle(yb); p=p+1,figure(p) % figure 3 subplot(221) plot(xb,ryb) subplot(222) plot(xb,iyb) subplot(223) plot(xb,ayb) subplot(224) plot(xb,pyb) sortfigs(1,p)
From: Mat Hunt on 20 Jun 2010 06:33
Hi Greg, I guess without knowing the full background to what I am doing. I am working on electrified flow down a channel with a moving pressure distribution. I am working on the weakly nonlinear theory which is a KdV like equation. I ignored the nonlinear terms and solve the equation via Fourier transforms (the equation contained a Benjamin-Ono term (a Hilbert transform term)) and that is where the |k| term comes in as I am taking the derivative of a Hilbert transform of e^{ikx} which comes out as sgn(k)*k, which is |k|, that is the origin of the |k|. The complex term is nothing but Rayleigh dispersion, which I have set to zero for the parameters I am using . I solved the equation numerically and I want to compare my answers for the linear case, and currently they don't match at all which is very worrying. Your solution was close to what I get for my code which is very encouraging. Mat |