Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
Drexel University or The Math Forum.


yubo
Posts:
1
Registered:
6/17/13


Re: FFT and differentiation
Posted:
Jun 17, 2013 4:34 PM


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/21)/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 npoint fft > > t = delt*(0:n1) > f = delf*(0:n1) % array when ffting > f = delf*(n/2:n/21) % fftshifted array, n even > f = delf*((n1)/2:(n1)/2) % fftshifted array, n odd > > with delf*delt = 1/n. (1) > > If you set > > n = 128; > delt = .01; > t = delt*(0:n1); > 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/21); % 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.3785e012 > > (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(((tt0).^2)/sigma^2); > exact_ans = (2*(tt0)/sigma^2).*a; > > and redoing the same process gives > > maxerror = 6.7076e010 > > 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 iffting 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 > 9094th Ave. N, Apt. D > Seattle, WA 98109 USA



