Inverse Fourier transforms
in [Matlab]
|
From: Greg Heath on 20 Jun 2010 15:16 On Jun 20, 6:33 am, "Mat Hunt" <hunt_...(a)hotmail.com> wrote: > 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. I am intrigued by the rapid oscillations in the real and imaginary part of ifft(Yb) opposed to the nice and smooth plot of abs(ifft(Yb)). Hope this helps. Greg
From: Mat Hunt on 21 Jun 2010 04:08 Greg Heath <heath(a)alumni.brown.edu> wrote in message <4cbcef89-1c16-4d58-8eab-7da8e0e7eff2(a)r27g2000yqb.googlegroups.com>... > On Jun 20, 6:33 am, "Mat Hunt" <hunt_...(a)hotmail.com> wrote: > > 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. > > I am intrigued by the rapid oscillations in the real and imaginary > part of ifft(Yb) opposed to the nice and smooth plot of abs(ifft(Yb)). > > Hope this helps. > > Greg So What do I need to do if I have the Fourier transform F(k), of a function, f as a function of the Fourier transform variable (say k, for example) back to the original variable x, say? I can send you my other numerical code which I used to compute the integral if you like. Mat
From: Greg Heath on 21 Jun 2010 07:36 On Jun 20, 3:16 pm, Greg Heath <he...(a)alumni.brown.edu> wrote: > On Jun 20, 6:33 am, "Mat Hunt" <hunt_...(a)hotmail.com> wrote: > -----SNIP > > > Your solution was close to what I get for my code which is very encouraging. > > I am intrigued by the rapid oscillations in the real and imaginary > part of ifft(Yb) opposed to the nice and smooth plot of abs(ifft(Yb)). % 1. Used trial and error to determine whether to use % fftshift or ifftshift (different for N odd) for % shifting from bipolar k to unipolar k: % Y = ifftshift(Yb) % and from unipolar x to bipolar x % y = fftshift(ifft(Y)) % 2. Size is reduced from 2001 to 81 to facilitate % observing plot details % 3. Compared yb = fftshift(ifft(Yb)) with y % a. yb is complex but y is real % b. They differ by a phase shift that is linear % in x % 5. Changing the use of fftshift and ifftshift results in % both yb and y being complex. The corresponding plots % for N = 2001 contain relatively wild oscillations (Left % as an excercise for the reader). close all, clear all, clc, p = 0 %Original Formulation kb = -100:.1:100; % bipolar k N0 = length(kb) % N0 = 2001; N0 is odd Yb = 16.2668*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2+9*abs(kb)); minmaxYb = minmax(Yb) % [ 0 0.11296 ] [Ybmax ibmax] = max(Yb) % [ 0.11296 1001] [Ymax imax] = max(fftshift(Yb)) % [ 0.11296 2001] [Ymax imax] = max(ifftshift(Yb)) % [ 0.11296 1 ] % Convenient Reformulation kb = -4 : 0.1 : 4; % bipolar k N = length(kb) % N = 81; N is odd dk = kb(2)-kb(1) % 0.1 k = dk*(0:N-1); % unipolar k Yb = 144*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2+9*abs(kb)); minmaxYb = minmax(Yb) % [ 0.0052749 1 ] Y = ifftshift(Yb); [Ymax imax] = max(Y) % [ 1 1 ] %Similar Transforms Ybp = 144*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2+9*kb); Ybm = 144*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2-9*kb); p = p+1, figure(p), hold on % figure 1 plot(kb,Ybp,'r') plot(kb,Ybm,'g') plot(kb,Yb) % Pole placements for analytic inversion % % % Four poles of Yb are at +/-kb0 and +/-conj(kb0) % No symbolic solution because of abs(k) % Consider numerical search starting at symbolic % solution for poles of Ybp and Ybm: kbp0 = double(solve('144 + 2*kb^4 - 12*kb^2 + 9*kb = 0','kb')) % kbp0 = 2.3977 - 1.794i % 2.3977 + 1.794i % -2.3977 - 1.5099i % -2.3977 + 1.5099i kpm0 = -kbp0; % function denom = denomfunc(kb) % % denom = 144+2*kb.^4-12*kb.^2+9*abs(kb) % % return kb0 = fsolve(@denomfunc,kbp0(1),optimset('fsolve')) % kb0 = 2.4756 - 1.7688i denomfunc(kb0) % 9.3344e-011 +1.3068e-010i denomfunc(-kb0) % 9.3344e-011 +1.3068e-010i denomfunc(conj(kb0)) % 9.3344e-011 -1.3068e-010i denomfunc(-conj(kb0)) % 9.3344e-011 -1.3068e-010i % Return to numerical inversion dx = 1/(N*dk) % 0.12346 x = dx*(0:N-1); minmax(x) % [ 0 9.8765 ] xb = x-dx*(N-1)/2; minmax(xb) % [ -4.9383 4.9383 ] yb = fftshift(ifft(Yb)); y = fftshift(ifft(Y)); absyb = abs(yb); absy = abs(y); maxabs(absy-absyb) % 5.0849e-017 [ymax imax] = max(absy) % [ 0.42304 41] ryb = real(yb); ry = real(y); iyb = imag(yb); iy = imag(y); p = p+1, figure(p) % figure 2 subplot(221) plot(xb,ryb,'b') title('Real(yb)') subplot(222) plot(xb,iyb,'b') title('Imag(yb)') subplot(223) plot(xb,ry,'r') title('Real(y)') subplot(224) plot(xb,iy,'r') title('Imag(y)') phaseyb = angle(yb); phasey = angle(y); p = p+1, figure(p) % figure 3 subplot(221) plot(xb,absyb,'b') title('Mag(yb)') subplot(222) plot(xb,phaseyb,'b') title('Phase(yb)') subplot(223) plot(xb,absy,'r') title('Mag(y)') subplot(224) plot(xb,phasey,'r') title('Phase(y)') sortfigs(1,p) Hope this helps. Greg
From: Mat Hunt on 21 Jun 2010 17:42 Hi Greg, So I need to do an fftshift after the ifft to get the answer? do I also need to do this with the variable? I can evaluate the integral using residue calculus if I ignore the |k| terms and that should provide me with another check on my working. I had a look at you code and there were some commands which didn't run and that matlab. Are you using an old version? Mat
From: Greg Heath on 21 Jun 2010 19:58
On Jun 21, 5:42 pm, "Mat Hunt" <hunt_...(a)hotmail.com> wrote: > Hi Greg, > > So I need to do an fftshift after the ifft to get the answer? Yes, fft and ifft assume the nonnegative intervals x = dx*(0:N-1) and k = dk*(0:N-1). fftshift and it's inverse ifftshift allow the translation to and from the bipolar interval with indexing (0:N-1)-ceil((N-1)/ 2). When N is even, fftshift and ifftshift yield the same results. As in my code, when N is odd, I have to use trial and error to figure out which one to use. For example, for N =6: fftshift([-3:2]) = [ 0 1 2 -3 -2 -1 ] ifftshift([-3:2]) = [ 0 1 2 -3 -2 -1 ] whereas for N = 5 fftshift([-2:2]) = [ 1 2 -2 -1 0 ] ifftshift([-2:2]) = [ 0 1 2 -2 -1 ] It may appear trivial, but a one bin translation can cause an imaginary part with rapid ocscillations to appear when N is large. >do I also need to do this with the variable? No, I have just given you the variable indexing. >I can evaluate the integral >using residue calculus if I ignore the |k| terms and that should provide > me with another check on my working. ??? Reread my post. I have already found the poles without ignoring the abs(k). > I had a look at you code and there were some commands which didn't run > and that matlab. Are you using an old version? Yes. It is old as dirt. >> ver ------------------------------------------------------------------------------------- MATLAB Version 6.5.1.199709 (R13) Service Pack 1 MATLAB License Number: XXXXXX Operating System: Microsoft Windows XP Version 5.1 (Build 2600: Service Pack 3) Java VM Version: Java 1.3.1_01 with Sun Microsystems Inc. Java HotSpot(TM) Client VM ------------------------------------------------------------------------------------- sortfigs is my own 3 line figure reordering function. Hope this helps. Greg |