Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » Software » comp.soft-sys.matlab

Topic: FFT and differentiation
Replies: 2   Last Post: Jun 17, 2013 4:34 PM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View  
yubo

Posts: 1
Registered: 6/17/13
Re: FFT and differentiation
Posted: Jun 17, 2013 4:34 PM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

If I change the function to exp(-2t), then it doesn't work.
The following is my code. X2 should be equal to X3. X1 should be equal to X4. But they are not equal. Could you help me explain why?
N = 512;
Fs = 100;
ts = 0;
dt = 0.01;
td = 2;
t = (ts:dt:td);
x1 = exp(-2*t);
x2 = -2*exp(-2*t);

figure
plot(t,x1,'r',t,x2,'b')
X1 =fftshift(fft(x1,N));
X2 = fftshift(fft(x2,N));
f_fft = (-N/2:N/2-1)/N*Fs;
omega = 2*pi*f_fft;
X3 = (X1.*omega)*1i;
X4 = X2./(omega*1i);

figure
plot(f_fft,abs(X1),'r',f_fft,abs(X2),'b')
legend('X1','X2')

figure
plot(f_fft,abs((X1)),'r',f_fft,abs((X2)),'b',f_fft,abs((X3)),'g',f_fft,abs((X4)),'m')
legend('X1','X2','X3 = iwX1','X4 = X2/iw')



Dave Goodmanson <dgoodmanson@worldnet.att.net> wrote in message <6inmnf$q22@bgtnsc03.worldnet.att.net>...
> PAR wrote:
> > > >: Hi,
> > > >:
> > > >: taking an FFT of a signal and then multiplying by i*w and then doing the
> > > >: inverse FFT is supposed to be equivalent to the derivative of the signal.
> > > >:
> > > >: Help !!!
> > > >:
> > > >: I can't get this to work on even the simplist of signal.

> <snip>
>
> PAR--
>
> For many functions it's certainly possible to differentiate by means of
> fft. I don't agree that the function has to be something obviously
> periodic such as a sine wave. The fft thinks *every* function is
> periodic, and many "nonperiodic" functions fit perfectly well into this
> scheme. (See the gaussian example below).
>
> To begin with, I believe neither Matthias's nor Hans's frequency array
> is correct. If delt is the time step, and delf is the frequency step,
> then for an n-point fft
>
> t = delt*(0:n-1)
> f = delf*(0:n-1) % array when ffting
> f = delf*(-n/2:n/2-1) % fftshifted array, n even
> f = delf*(-(n-1)/2:(n-1)/2) % fftshifted array, n odd
>
> with delf*delt = 1/n. (1)
>
> If you set
>
> n = 128;
> delt = .01;
> t = delt*(0:n-1);
> T = delt*n; % length of the time record
> a = sin(2*pi*10*t/T); % periodic in T
>
> and proceed basically as Matthias did (changing the notation slightly),
> things turn out OK:
>
> b = fftshift(fft(a));
> delf = 1/T; % same as (1)
> f = delf*(-n/2:n/2-1); % n is even
> w = 2*pi*f;
> c = b.*(i*w);
> da_dt = ifft(ifftshift(c));
> plot(t,da_dt);
> % find the error compared with the actual answer:
> exact_ans = (2*pi*10/T)*cos(2*pi*10*t/T);
> maxerror = max(abs(da_dt - exact_ans))
>
> maxerror = 1.3785e-012
>
> (If your version of Matlab doesn't have ifftshift, you can use fftshift
> instead if n is even, or use Hans's suggested mfile, or make a copy of
> fftshift and change all occurrences of "ceil" to "floor".)
>
> A gaussian is an example of a "nonperiodic" function. Replacing a and
> exact_ans with
>
> t0 = .5; sigma = .1;
> a = exp(-((t-t0).^2)/sigma^2);
> exact_ans = -(2*(t-t0)/sigma^2).*a;
>
> and redoing the same process gives
>
> maxerror = 6.7076e-010
>
> so the method definitely can work. The key is that both the gaussian
> and its transform fall off quickly before reaching the edges of their
> respective arrays (in the case of frequency, the shifted array with f=0
> in the center). The fft considers array end points b(1) and b(n) to be
> next door neighbors, but the function is so small and flat out there
> that neither the function nor its derivative has a step discontinuity.
> Even after multiplying by iw, you can say the same thing about the
> resulting function c. So ifft-ing back into the time domain gives the
> right answer.
>
> The method won't work for a function that falls off too slowly either in
> time or frequency (you can't differentiate t^2 this way). There is the
> additional requirement that c(f) = i*w*b(f) also falls off quickly
> enough. A few functions b(f) are inherently bad and fall off like 1/|w|
> or some such, but many functions do fall off quickly enough for this
> method to work.
>
> Note that you are free to scale your time array so as to vary the
> percentage of the array taken up by a(t). The wider you make a(t), the
> narrower b(f) becomes. This is good if you intend to multiply by iw.
> But if the function a(t) gets too wide, there will be a discontinuity at
> the two neighboring points a(1) and a(n) of the time array. Also, if
> the frequency function is too narrow, the frequency array effectively
> becomes too coarse to properly describe b(f). It's a tradeoff to get
> the right time and frequency array scale.
>
> --Dave
>
> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
>
> David Goodmanson dgoodmanson@worldnet.att.net
> 909-4th Ave. N, Apt. D
> Seattle, WA 98109 USA




Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.