[PDF] [PDF] Notes on FFT-based differentiation - MIT Mathematics

[PDF] La transformation de Fourier rapide sous Matlab : fft - ENSTA Paris

sous Matlab : fft, ifft, fftshift, ifftshift Karsten Plamann, février 2018 1 Définition et syntaxe de fft et ifft1 Y = fft(X) et X = ifft(Y) calculent la transformation de Fourier 

[PDF] FFT Tutorial

Matlab's FFT function is an effective tool for computing the discrete Fourier transform of a signal The following code examples will help you to understand the 

[PDF] Introduction A Simple Example

Introduction In this tutorial, we will discuss how to use the fft (Fast Fourier Transform) command plot the frequency spectrum using the MATLAB fft command

[PDF] CS425 Lab: Frequency Domain Processing 1 Discrete Fourier

the fast Fourier transform (FFT) is a fast algorithm for computing the discrete MATLAB has three related functions that compute the inverse DFT: 0 ifft 1 ifft2

[PDF] Analyse de Fourier en pratique

Spectre d'un signal sinusoïdal sous Matlab avec l'application SignalAnalyzer Les signaux sont échantillonnés Les spectres sont calculés par FFT : calcul de 

[PDF] Fourier series in MATLAB

Find the Fourier coefficients using your MATLAB function: plot the Fourier coefficients vs frequency 5 The FFT Despite the fact that we presented the discrete 

[PDF] Notes on FFT-based differentiation - MIT Mathematics

In Matlab, the 1/N normalization is moved from the DFT (Matlab's fft function) to the IDFT (Matlab's ifft function), which doesn't affect any of the procedures in this 

[PDF] IMAGE PROCESSING IN FREQUENCY DOMAIN USING MATLAB

15 sept 2008 · FFT: Fast Fourier Transform IFFT: Inverse Fast Fourier Transform MSQE: Mean Square Quantization Error MATLAB: Matrix Laboratory ROI:

[PDF] RECOVERING SIGNAL FREQUENCIES WITH THE - Kenneth L Ho

The DFT is computed in Matlab by the command fft, for the fast Fourier X = FFTFREQ(X,T,DT) returns the discrete Fourier transform of X sampled over the  

[PDF] The Fundamentals of FFT-Based Signal Analysis and Measurement

The Fast Fourier Transform (FFT) and the power spectrum are powerful tools for analyzing and measuring signals from plug-in data acquisition (DAQ) devices

[PDF] fiche de lecture a cp

[PDF] fiche de lecture compréhension cp a imprimer

[PDF] fiche de lecture cp a imprimer pdf

[PDF] fiche de lecture cp gratuite a imprimer

[PDF] fiche de lecture cp pdf

[PDF] fiche de lecture cp son a

[PDF] fiche de lecture pour cp a imprimer

[PDF] fiche de poste comptable

[PDF] fiche de poste en anglais

[PDF] fiche de poste exemple word

[PDF] fiche de poste modèle

[PDF] fiche de vocabulaire français par thème pdf

[PDF] fiche lecture ce1 pdf

[PDF] fiche lecture ce2 pdf

[PDF] fiche pedagogique de l'enseignant primaire au maroc

Notes on FFT-based differentiation

Steven G. Johnson, MIT Applied Mathematics

Created April, 2011, updated May 4, 2011.

Abstract

A common numerical technique is to differentiate some sampled functiony(x)via fast Fourier trans-

forms (FFTs). Equivalently, one differentiates an approximate Fourier series. Equivalently, one differen-

tiates a trigonometric interpolation. These are also known asspectraldifferentiation methods. However,

the implementation of such methods is prey to common confusions due to thealiasingphenomenon inherent in sampling, and the treatment of the maximum-frequency (Nyquist) component is especially

tricky. In this note, we review the basic ideas of spectral differentiation (for equally spaced samples) and

a correct resolution of the aliasing problems for implementing operations dydx ,d2ydx 2,ddx c(x)dydx , andr2.

One surprising consequence of aliasing for the Nyquist component is that the sampled second-derivative

operation isnotequivalent to two sampled first-derivative operations! The operationddx c(x)dydx is par-

ticularly subtle because of the interplay of aliasing and desired algebraic properties of the differential

operator. (Readers who simply want to know what algorithms to follow, rather than the details of the derivations, can read the first two introductory paragraphs and then skip to the labelledAlgorithm boxes.)

1 Background

In spectral methods for differential equations, considering one dimension here for simplicity, one has a

periodic functiony(x)with periodLthat one conceptually expands as a Fourier series: y(x) =1X k=1Y ke2iL for Fourier coefficientsYk=1L 0e2iL kxy(x)dx. One then wishes to apply differential operators likeddx 2dx

2, and more generallyddx

c(x)ddx [for some periodic functionc(x)]. Differentiation is performed term-by- term

1in Fourier domain, and multiplication by functionsc(x)is done by transforming back to space domain

for the multiplication (equivalent to convolution in Fourier domain). For example, ddx y(x) =y0(x) =1X k=1 2iL kYk 2iL kx;1

You may have been disturbed by rumors from real analysts that differentiating a Fourier series term-by-term is not always

valid. Don"t worry. For one thing, even if we required pointwise convergence of the Fourier series, differentiating term-by-term

is valid for they(x)that typically appear in practice, whereyis usually continuous andy0is at least piecewise differentiable.

More generally, in practical applications of Fourier analysis, such as for PDEs, we are ordinarily not interested in pointwise

convergence-we only care about "weak" convergence (equality when both sides are integrated against smooth, localized test

functions), in which case differentiation term-by-term is always valid for any generalized functiony.

which is just a pointwise multiplication of eachYkby a term proportional tok. To implement this on a computer, one approximates the Fourier series by adiscrete Fourier transform (DFT). That is, we replace the functiony(x)byNdiscretesamplesyn=y(nL=N)forn= 0;1;:::;N1, andYkis approximated by: k=1N N1X n=0y ne2iN nk; a DFT (although sign and normalization conventions for the DFT vary

2). The inverse transform (IDFT) to

computeynfromYkis very similar: n=N1X k=0Y ke+2iN nk:

On a computer, the nice thing about these expressions is that all theYkcan be computed from theyn, or

vice versa, in(NlogN)operations by afast Fourier transform(FFT) algorithm. This makes working with

Fourier series practical: we can quickly transform back and forth between space domain [where multiplying

by functions likec(x)is easy] and Fourier domain [where operations like derivatives are easy]. Almost

invariably, FFT implementations compute DFTs and IDFTs in forms similar to these equations, with the kcoefficients arranged "in order" fromk= 0toN1, and this ordering turns out to make the correct implementation of FFT-based differentiation more obscure. In order to compute derivatives likey0(x), we need to do more than expressyn. We need to use the

IDFT expression to define a continuous interpolation between the samplesyn-this is calledtrigonometric

interpolation-and then differentiate this interpolation. At first glance, interpolating seems very straightfor-

ward: one simply evaluates the IDFT expression at non-integern2R. This indeed definesaninterpolation,

but it is not theonlyinterpolation, nor is it thebestinterpolation for this purpose. The reason that there

is more than one interpolation is due toaliasing: any terme+2iN nkYkin the IDFT can be replaced by +2iN n(k+mN)Ykfor any integermand still give thesamesamplesyn, sincee2iN nmN=e2inm= 1for any integersmandn. Essentially, adding themNterm tokmeans that the interpolated functiony(x)just oscillatesmextra times in between the sample points, i.e. betweennandn+ 1, which has no effect on

nbut has a huge effect on derivatives. To resolve this ambiguity, one imposes additional criteria-e.g. a

bandlimited spectrum and/or minimizing some derivative of the interpolatedy(x)-and we will then explore

the consequences of the disambiguated interpolation for various derivative operations.

2 Trigonometric interpolation

What is the right way to do trigonometric interpolation? We can define a more arbitrary interpolated functiony(x)(but still limited to onlyNfrequency components) by substitutingn=Nx=Linto the IDFT2

In our FFTW implementation of FFTs, this DFT is termed anFFTW_FORWARDtransform, but the1=Nnormalization is

omitted and must be supplied by the user. (In computing derivatives by the procedures described below, the1=Nfactor can

be combined at no cost with the multiplications by2=L.) The IDFT is termed anFFTW_BACKWARDtransform. In Matlab, the

1=Nnormalization is moved from the DFT (Matlab"sfftfunction) to the IDFT (Matlab"sifftfunction), which doesn"t affect

any of the procedures in this note. Beware that Matlab uses indices that start with 1 rather than 0, however: ifz = fft(y)in

Matlab, thenz(k+1)corresponds toNYk.

and allowing any arbitrary aliasing integermkfor eachYk:3 y(x) =N1X k=0Y ke2iL (k+mkN)x: Regardless of the values ofmk, this gives the same sample pointsyn=y(nL=N), but changingmkgreatly

modifiesy(x)in between the samples. In order to uniquely determine themk, a useful criterion is that we

wish tooscillate as little as possiblebetween the sample pointsyn. One way to express this idea is to assume

thaty(x)isbandlimitedto frequenciesjk+mkNj N=2. Another approach, that gives the same result

given this form ofy(x)(i.e., what are the "best"Nfrequency components?), is tominimize the mean-square

slope jy0(x)j2dx=1L N1X k=02iL (k+mkN)Yke2iL (k+mkN)x2 N1X k=0N1X 0=01L 02iL (k+mkN)Yke2iL (k+mkN)x2iL (k0+mk0N)Yk0e2iL (k0+mk0N)xdx 2N1X k=0jYkj2(k+mkN)2; where in the last step we have used the orthogonality of the Fourier basis functions: 2iL (k+mkN)xe2iL (k0+mk0N)xdx=

0ifk+mkN6=k0+mk0N, which for0k;k0< Nboils down to being zero fork6=k0(and givingL

fork=k0). From this expression, for a given set of coefficientsYk, the mean-square slope is minimized by

choosingmkthat minimizes(k+mkN)2for eachk.3

Unfortunately, this is actually not the most general possible aliasing, because we are assigning onlyonemkto eachk. More

generally, we could break eachYkamong arbitrarily many aliasesk+mNwith amplitudesum;kYk, such thatP mum;k= 1. That is, we could write the most general trigonometric interpolation as: y(x) =N1X k=0Y m=1u m;ke2iL (k+mN)x:

The mean-square slope is then

2PN1 k=0jYkj2P1 m=1jum;kj2(k+mN)2. We can eliminate theP mum;k= 1constraint by solving foru0;k= 1P m6=0um;k, yielding a mean-square slope 2N1X k=0jYkj22 m6=0u m;k 2+X m6=0jum;kj2(k+mN)23

Minimizing this unconstrained expression is a straightforward calculus exercise, and yieldsum;k=ak=(k+mN)2, where

k= 1=P

m1(k+mN)2: an infinite number of frequency components! (In fact, this interpolation is especially problematic,

because itssecondderivativey00actually diverges! By allowing infinitely many aliases, we are allowinganyperiodic function,

and the minimial-slope periodic interpolant is simplypiecewise-linearinterpolation, with infinitey00at the "kinks" at each

sample.) To restrict ourselves to a single nonzeroum;kfor eachk, one possiblity is to restrict ourselves to one frequency

jk+mNjperYkso that we are asking for the "best" choice ofNfrequency components in the minimal-slope sense. Of course,

if we further assume that the spectrum is bandlimited tojk+mNj N=2, that completely determines the aliasing question

by itself (except for thek=N=2component, which can be determined by e.g. requiring that the interpolant be real given real

n). Alternatively, rather than anya prioribandlimiting assumption, we can obtain the desired result below if we minimize

the mean-squarep-th derivativejy(p)(x)j2dx, and then take thep! 1limit. If0k < N=2, then(k+mkN)2is minimized formk= 0. IfN=2< k < N, then(k+mkN)2is minimized formk=1. Ifk=N=2(for evenN), however, there is an ambiguity: eithermk= 0or1 gives the same value(k+mkN)2= (N=2)2. For thisYN=2term (the "Nyquist" term), we need to consider a more general class of aliasing: we can arbitrarily split up theYN=2term betweenm= 0[2iL x= +iL Nx, positive frequency] andm=1[2iL (N2

N)x=iL

Nx, negative frequency]:

quotesdbs_dbs2.pdfusesText_4