Rusty
Posts:
108
Registered:
6/16/05


Re: derivative of discrete fourier transform interpolation
Posted:
Oct 1, 2005 5:35 AM


<stevenj@alum.mit.edu> wrote in message news:1128133526.658039.28980@g44g2000cwa.googlegroups.com... > Peter Spellucci wrote: >> the FFT gives a real sin/cos series only for special N and in this case >> your sum >> is running from N/2 to N/2 . the real sin/cos (Fourier ) sum can then be >> differentiated in the usual manner with no trouble. > > Special N?? The DFT of real inputs x_n, for any N, can be interpreted > as the amplitudes of a real sin/cos series (with real derivatives). > > To see this, you pair the output X_k with X_{nk}=X_k^*, by the > Hermitian symmetry of the output. Equivalently, you are using the > aliasing property to think of the X_{nk} output, for k > n/2, as the > X_{k} amplitude (i.e. a negative frequency k), and so you get > complexconjugate pairs of sinusoids. (The k=n/2 Nyquist element, for > even n, must be treated specially. Since it is purely real, it can be > thought of as 1/2 k=n/2 and 1/2 k=+n/2.) > > In order to take the derivative, you need to realize that the > interpolation implied by the DFT, and hence the slope, is not unique > because of aliasing. Normally, however, you want the interpolation > corresponding to exactly the aliasing described above: the k > n/2 > outputs are *negative* frequency amplitudes. > > This choice means that your frequencies run from N/2+1 to N/2 (for > even N). Not only does it guarantee real derivatives from real > inputs, but it also corresponds to the interpolation with the *minimal* > meansquare slope.
I was trying to remember the trick with the middle coefficient. Here is the MatLab code which puts half its amplitude in the upper band and half in the lower, giving a real valued interpolation. The ripples are a bit excessive unless the function starts off well band limited, which could be done by FFT.
PS: why is Matlab alone in not using array indices starting at zero ? It makes a pigs ear of DFT's Rusty
N=16
x=(1:N)'
x=x
X=fft(x)
s=0.25
ix=0;
for k=1:N/2
ix=ix+X(k)*exp(2*pi*1i*(k1)*s/N);
end
for k=1:N/21
ix=ix+X(N+1k)*exp(2*pi*1i*k*s/N);
end
ix=ix+X(N/2+1)*cos(2*pi*s/2);
ix=ix/N

