Thread Subject: integrating with FFT

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!

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

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

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