MATLAB Central Newsreader - integrating with FFT

integrating with FFT

Paul — Thu, 27 Dec 2007 01:56:07 -0500

You can integrate in the frequency domain by dividing the
complex FFT by iw. I've been trying this but to no avail
since I seem to be getting hung up on the symmetry of the
FFT around the nyquist.

For example:

x = 1:10;

and the integral of x is trivial.

But I try to do this via fft, division by iw and then
inverse fft and end up with junk. I recall seeing this
done someplace before but I checked the archives and found
nothing.

Can anybody show how this is done!

Re: integrating with FFT

Greg Heath — Thu, 27 Dec 2007 04:09:26 -0500

On Dec 26, 8:56=A0pm, "Paul " <p...@ceri.memphis.edu> wrote:
> You can integrate in the frequency domain by dividing the
> complex FFT by iw. =A0I've been trying this but to no avail
> since I seem to be getting hung up on the symmetry of the
> FFT around the nyquist.
>
> For example:
>
> x =3D 1:10;
>
> and the integral of x is trivial.
>
> But I try to do this via fft, division by iw and then
> inverse fft and end up with junk. =A0 I recall seeing this
> done someplace before but I checked the archives and found
> nothing.
>
> Can anybody show how this is done!

The answer is probably in the details of your code. Post it
and useful advice will probably follow.

Hope this helps.

Greg

Re: integrating with FFT

Greg Heath — Thu, 27 Dec 2007 08:49:58 -0500

On Dec 26, 8:56 pm, "Paul " <p...@ceri.memphis.edu> wrote:
> You can integrate in the frequency domain by dividing the
> complex FFT by iw.

With caveats.

To be consistent with the periodic DFT, make the following
assumptions

F(iw) = INT(0,T){ dt f(t)*exp(-iwt) },
G(iw) = INT(0,T){ dt g(t)*exp(-iwt) },
g(t) = g(0) + INT(0,t){dt f(t) },

Since f(t) = g'(t), integrating by parts yields

F(iw) = INT(0,T){ dt g'(t)*exp(-iwt) }
= g(T)*exp(-iwT)- g(0) + iw*G(iw)

Since F = iw*G must hold for all w, the following
conditions are required

1. g(0) = 0
2. g(T) = F(0) = 0

For both the sampled CFT and the DFT, the division by
iw at w = 0 is problematical. Therefore,

3. The point w = 0 should be excluded by assuming
G(0) = 0 = F(0) .

> I've been trying this but to no avail
> since I seem to be getting hung up on the symmetry of the
> FFT around the nyquist.
>
> For example:
>
> x = 1:10;
>
> and the integral of x is trivial.

However, the average value of x is not zero.

> But I try to do this via fft, division by iw and then
> inverse fft and end up with junk. I recall seeing this
> done someplace before but I checked the archives and found
> nothing.
>
> Can anybody show how this is done!

N = 10
dt = 1
t = (0:N-1)*dt;
x = t+1; % need to subtract the mean value
m = mean(x) % m = 99/18
x = x - mean(x);
y = ((t+1).^2/2 - (99/18)*t -1/2;

T = N*dt
df = 1/T
f = (0:N-1)*df;
X = fft(x);
Y = fft(y);

P = floor((N+1)/2)
fP = (P-1)*df
fb = f-fP;
Xb = fftshift(X);
Yb = fftshift(Y);

iwb = i*2*pi*fb;

%Since iwb(P) = 0, using Xb./iwb is problematic.
% Instead, investigate

Err = Yb-iwb.*Xb;

Hope this helps.

Greg

Re: integrating with FFT

Greg Heath — Thu, 14 Feb 2008 05:16:14 -0500

On Dec 27 2007, 3:49=A0am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Dec 26, 8:56 pm, "Paul " <p...@ceri.memphis.edu> wrote:
>
> > You can integrate in the frequency domain by dividing the
> > complex FFT by iw.
>
> With caveats.
>
> To be consistent with the periodic DFT, make the following
> assumptions
>
> F(iw) =3D INT(0,T){ dt f(t)*exp(-iwt) },
> G(iw) =3D INT(0,T){ dt g(t)*exp(-iwt) },
> g(t) =3D g(0) + INT(0,t){dt f(t) },
>
> Since f(t) =3D g'(t), integrating by parts yields
>
> F(iw) =3D INT(0,T){ dt g'(t)*exp(-iwt) }
> =A0 =A0 =A0 =A0 =A0 =3D g(T)*exp(-iwT)- g(0) + iw*G(iw)
>
> Since F =3D iw*G must hold for all w, the following
> conditions are required
>
> 1. g(0) =3D 0
> 2. g(T) =3D F(0) =3D 0
>
> For both the sampled CFT and the DFT, the division by
> iw at w =3D 0 is problematical. Therefore,
>
> 3. The point w =3D 0 should be excluded by assuming
> =A0 =A0 G(0) =3D 0 =3D F(0) .
>
> > I've been trying this but to no avail
> > since I seem to be getting hung up on the symmetry of the
> > FFT around the nyquist.
>
> > For example:
>
> > x =3D 1:10;
>
> > and the integral of x is trivial.
>
> However, the average value of x is not zero.
>
> > But I try to do this via fft, division by iw and then
> > inverse fft and end up with junk. =A0 I recall seeing this
> > done someplace before but I checked the archives and found
> > nothing.
>
> > Can anybody show how this is done!
>
> N =3D 10
> dt =3D 1
> t =3D (0:N-1)*dt;
> x =3D t+1; =A0 =A0 =A0 =A0 =A0 =A0% need to subtract the mean value
> m =3D mean(x) =A0 % m =3D 99/18
> x =3D x - mean(x);
> y =3D ((t+1).^2/2 - (99/18)*t -1/2;
>
> T =3D N*dt
> df =3D 1/T
> f =3D (0:N-1)*df;
> X =3D fft(x);
> Y =3D fft(y);
>
> P =3D floor((N+1)/2)
> fP =3D (P-1)*df
> fb =3Df-fP;

Correction

Q =3D ceil((N+1)/2)
fQ =3D (Q-1)*df
fb =3D f-fQ;

> Xb =3D fftshift(X);
> Yb =3D fftshift(Y);
>
> iwb =3D i*2*pi*fb;
>
> %Since iwb(P) =3D 0, using Xb./iwb is problematic.

Since iwb(Q) =3D 0, using Xb./iwb is problematic.

> % Instead, investigate
>
> Err =3D Yb-iwb.*Xb;
>
> Hope this helps.
>
> Greg