辅导data留学生编程、讲解Python编程、Python程序调试
- 首页 >> Matlab编程 General tips for computing and plotting discrete fourier transforms (DFT)
Let g(t) be some time signal that’s sampled at to get a discrete array/list g = [ … ]
You don’t need to code your own DFT, use numpy: A = numpy.fft.fft(g) * dt
You should also use numpy.fft.fftshift( A ) to shift the fft output such that the 0-frequency component is centered (see why
here https://docs.scipy.org/doc/numpy/reference/routines.fft.html#background-information) which you probably want when plotting
You can use f_axis = numpy.fft.fftshift( numpy.fft.fftfreq(len(g), dt) ) to create the frequency axis for plotting
the shifted spectrum
Keep in mind the fft output is in general complex, so to compare two fourier transforms (e.g. DFT vs CFT) you should compare either
the real and imaginary parts ( ), or the phase and amplitude ( ). Amplitude plots are most useful for this lab, show
all 4 aspects though if you want.
Note numpy fft assumes the time signal starts at , if yours doesn’t you should center it at zero. If you don’t then the complex
components (x and y, or phase) will be off, but the amplitude should not change (why? analytically, recall that axis shifts in either
domain are equivalent to complex exponential scaling, which has amplitude 1, in the other domain).
If your time signal g is centered at zero, a hack to ‘rotate it’ to start at zero (and then take the fft and fftshift that) is to do: fftshift(
fft( ifftshift(g) ) ) (you may see mention of this online)
1 Fourier transform of Gaussian Functions (6 pts)
A common function used for the convolution of time series data is the Gaussian function
where is the half duration.
1. Plot for and sec on the same graph with domain and .
2. The analytical formula for the Fourier transform of is
Compute the discrete Fourier transform (DFT) for both sampled time series, and compare them to the analytical for both
's on the same graph.
Hints:
As numpy fft assumes signal starts from time 0, you can use the shift property of Fourier transform to first shift the to start
from zero, and after fftshift(fft()) operations, multiply the spectrum by complex exponential sinusoid function.
You need to sample the theoretical curve with w_axis = 2*pi*f_axis , or else rewrite it as if you'd rather
sample it with f_axis
As a guide (so you can be confident of your fft utilization for the remainder of the lab), we expect that the amplitudes (use
numpy.abs(…) ) of the discrete FT and the continuous FT essentially match. The phase won’t necessarily match.
3. Comment on the effect of filtering a general input time function by (i.e. convolution of with ), and explain the
difference in filtered output after applying Gaussian functions with or secs.
4. Comment on how this is related to the time-frequency uncertainty principle (a signal cannot be infinitesimally sharp both in time and
frequency).
2 Fourier transform of Window Functions (6 pts)
A continuous time signal can be truncated into a signal of finite length by window functions :
Typical window functions include:
Boxcar function
Hann window
Now let seconds, and sample both window functions by seconds:
1. Plot both window functions on the same graph.
2. Calculate the Fourier transform of both functions by numpy fft() . Pay extra attention to how you interpret the corresponding
frequencies of output results from python. (Hint: fftshift() may be useful. Also pay attention to the length of the input signal (> 10
sec), as it dictates the frequency resolution for the spectrum.)
3. Plot the Fourier transform of both functions in the appropriate frequency range on the same graph.
4. Based on the FTs, comment on the effect of truncating a continuous time series by either window on its frequency spectrum
compared to the original spectrum .
5. Speculate on the advantages and disadvantages of boxcar and Hann window functions for truncation.
3 Radial Distribution Function (12 pts)
Background
Liquids have no fixed internal structure. Yet they do have some short range order in the sense that they have preferred intermolecular
spacings, which are determined by the locations of minima in the intermolecular potentials. The microscopic structure of liquids is often
characterized by a quantity known as the Radial Distribution Function , which is essentially the probability (Relative to the average
probability, which means that tends to 1 at large , where the neighbour is too far away to feel any interaction.) that a molecule has a
neighbouring molecule at distance . Typically shows a value that approaches zero at small since molecules cannot occupy the
same space; it also shows a peak at the preferred distance of nearest neighbours, and secondary peaks at preferred distances of more
distant neighbours. If a suitable collimated beam of particles (e.g. X-rays or neutrons) is sent through a sample of the liquid, some of the
particles are scattered. The number of particles scattered through a given angle is related to the Fourier Transform of evaluated at the
wavenumber corresponding to the transverse momentum transfer associated with that scattering angle. Kittel derives this relationship in
Chapter 17 of Introduction to Solid State Physics.
If this all sounds complicated, all you need to know here is that something called the Structure Factor is effectively measured by
looking at the scattered intensity as a function of scattering transverse wavenumber (proportional to scattering angle), and that the Radial
Distribution Function is related to it by
where is liquid number density (number of atoms per unit volume, computable from the three constants mentioned in the introduction),
is wavenumber, and is radius.
1. You may have noticed some resemblance between expression (6) and the Fourier transform. First show that the integration part
can be rewritten as
Hint: The structure factor S(k) is even, since there should be no reason why scattering intensity would be different for one direction (+k)
compared to its opposite (−k). Using the fact that S(k) is even may be useful.
2. Now we can make some connections between the Radial Transfer Function and the Fourier Transform, if we substitute and
. What is the Fourier transform of ? Is a real, imaginary or general complex function? Is it even or odd? How will
these affect ? Is that what you expect? Plot as a function of ranging from Å to Å based on argon.py (i.e.
import and use the variables defined there).
Hint: In constructing from argon.py , you should make an "even" array twice the length (minus 1) of YanData. YanData
represents the structure factor (i.e. ) for argon sampled at the defined in the argon.py file. It's specifically sampled
from k = 0 to k = len(YanData) * dk , so create an even function out to the same length in the negative direction (i.e. the "kaxis"
it's sampled on would be len(YanData)*dk, ... 0, ..., + len(YanData)*dk ).
3. Write a Python function [gn, rn] = RDFcalc(S, dk, rho) to calculate Radial Distribution Function from Structure Factor
data, sampled at , and density , and output the RDF vector and its corresponding radial distance vector .
Hint: for Python fft() and ifft() functions, realize that the values of the Fourier Transform corresponding to negative
frequencies are stored in the second half of the arrays given to ( ifft ) or obtained from it ( fft ). You also have to study the
difference between the DFT and FT to multiply the right factors.
4. With the data provided in argon.py , compute the corresponding Radial Distribution Function . Plot your results for from to
Å. Over what range of radius can you trust your result?
Hint: To check if your results make sense, recall that is related to the probability that a molecule has a neighbouring molecule at
distance , therefore, should be close to when , i.e. two molecules can not occupy the same space, and you can set
. Recall . Also note the unit used in formula (6).
5. From the you computed, estimate the average molecular radius of liquid argon. Give your reasoning and state what accuracy
you can justify for your estimation.
6. Now we explore the effect of sampling range. Yan sampled in wavenumber out to Å , and he could have saved
himself work by not collecting as much data, i.e., reducing . But how much could he have reduced the sampling length ,
while still seeing distinct peaks in the Radial Distribution Function? Also explain theoretically what you observe.
Hint: Plot on top of the obtained in part 4, the 's you compute for a series of values. You can try half each time to
look for changes. For the theoretical explanations for part 6 and 7, realize the interchangability of ( ).
7. To explore the effect of data sampling, let's assume Yan decided to save his work by sampling less often (i.e. increasing ). How
large a can he use to be able to still recover the first two peaks clearly? State your answers and a theoretical justification for what
you expect to see if you increase too much.
Hint: Plot on top of the obtained from argon.py data, the 's you obtain when you subsample the same dataset. Try
doubling each time to observe the effect of coarser sampling.
Let g(t) be some time signal that’s sampled at to get a discrete array/list g = [ … ]
You don’t need to code your own DFT, use numpy: A = numpy.fft.fft(g) * dt
You should also use numpy.fft.fftshift( A ) to shift the fft output such that the 0-frequency component is centered (see why
here https://docs.scipy.org/doc/numpy/reference/routines.fft.html#background-information) which you probably want when plotting
You can use f_axis = numpy.fft.fftshift( numpy.fft.fftfreq(len(g), dt) ) to create the frequency axis for plotting
the shifted spectrum
Keep in mind the fft output is in general complex, so to compare two fourier transforms (e.g. DFT vs CFT) you should compare either
the real and imaginary parts ( ), or the phase and amplitude ( ). Amplitude plots are most useful for this lab, show
all 4 aspects though if you want.
Note numpy fft assumes the time signal starts at , if yours doesn’t you should center it at zero. If you don’t then the complex
components (x and y, or phase) will be off, but the amplitude should not change (why? analytically, recall that axis shifts in either
domain are equivalent to complex exponential scaling, which has amplitude 1, in the other domain).
If your time signal g is centered at zero, a hack to ‘rotate it’ to start at zero (and then take the fft and fftshift that) is to do: fftshift(
fft( ifftshift(g) ) ) (you may see mention of this online)
1 Fourier transform of Gaussian Functions (6 pts)
A common function used for the convolution of time series data is the Gaussian function
where is the half duration.
1. Plot for and sec on the same graph with domain and .
2. The analytical formula for the Fourier transform of is
Compute the discrete Fourier transform (DFT) for both sampled time series, and compare them to the analytical for both
's on the same graph.
Hints:
As numpy fft assumes signal starts from time 0, you can use the shift property of Fourier transform to first shift the to start
from zero, and after fftshift(fft()) operations, multiply the spectrum by complex exponential sinusoid function.
You need to sample the theoretical curve with w_axis = 2*pi*f_axis , or else rewrite it as if you'd rather
sample it with f_axis
As a guide (so you can be confident of your fft utilization for the remainder of the lab), we expect that the amplitudes (use
numpy.abs(…) ) of the discrete FT and the continuous FT essentially match. The phase won’t necessarily match.
3. Comment on the effect of filtering a general input time function by (i.e. convolution of with ), and explain the
difference in filtered output after applying Gaussian functions with or secs.
4. Comment on how this is related to the time-frequency uncertainty principle (a signal cannot be infinitesimally sharp both in time and
frequency).
2 Fourier transform of Window Functions (6 pts)
A continuous time signal can be truncated into a signal of finite length by window functions :
Typical window functions include:
Boxcar function
Hann window
Now let seconds, and sample both window functions by seconds:
1. Plot both window functions on the same graph.
2. Calculate the Fourier transform of both functions by numpy fft() . Pay extra attention to how you interpret the corresponding
frequencies of output results from python. (Hint: fftshift() may be useful. Also pay attention to the length of the input signal (> 10
sec), as it dictates the frequency resolution for the spectrum.)
3. Plot the Fourier transform of both functions in the appropriate frequency range on the same graph.
4. Based on the FTs, comment on the effect of truncating a continuous time series by either window on its frequency spectrum
compared to the original spectrum .
5. Speculate on the advantages and disadvantages of boxcar and Hann window functions for truncation.
3 Radial Distribution Function (12 pts)
Background
Liquids have no fixed internal structure. Yet they do have some short range order in the sense that they have preferred intermolecular
spacings, which are determined by the locations of minima in the intermolecular potentials. The microscopic structure of liquids is often
characterized by a quantity known as the Radial Distribution Function , which is essentially the probability (Relative to the average
probability, which means that tends to 1 at large , where the neighbour is too far away to feel any interaction.) that a molecule has a
neighbouring molecule at distance . Typically shows a value that approaches zero at small since molecules cannot occupy the
same space; it also shows a peak at the preferred distance of nearest neighbours, and secondary peaks at preferred distances of more
distant neighbours. If a suitable collimated beam of particles (e.g. X-rays or neutrons) is sent through a sample of the liquid, some of the
particles are scattered. The number of particles scattered through a given angle is related to the Fourier Transform of evaluated at the
wavenumber corresponding to the transverse momentum transfer associated with that scattering angle. Kittel derives this relationship in
Chapter 17 of Introduction to Solid State Physics.
If this all sounds complicated, all you need to know here is that something called the Structure Factor is effectively measured by
looking at the scattered intensity as a function of scattering transverse wavenumber (proportional to scattering angle), and that the Radial
Distribution Function is related to it by
where is liquid number density (number of atoms per unit volume, computable from the three constants mentioned in the introduction),
is wavenumber, and is radius.
1. You may have noticed some resemblance between expression (6) and the Fourier transform. First show that the integration part
can be rewritten as
Hint: The structure factor S(k) is even, since there should be no reason why scattering intensity would be different for one direction (+k)
compared to its opposite (−k). Using the fact that S(k) is even may be useful.
2. Now we can make some connections between the Radial Transfer Function and the Fourier Transform, if we substitute and
. What is the Fourier transform of ? Is a real, imaginary or general complex function? Is it even or odd? How will
these affect ? Is that what you expect? Plot as a function of ranging from Å to Å based on argon.py (i.e.
import and use the variables defined there).
Hint: In constructing from argon.py , you should make an "even" array twice the length (minus 1) of YanData. YanData
represents the structure factor (i.e. ) for argon sampled at the defined in the argon.py file. It's specifically sampled
from k = 0 to k = len(YanData) * dk , so create an even function out to the same length in the negative direction (i.e. the "kaxis"
it's sampled on would be len(YanData)*dk, ... 0, ..., + len(YanData)*dk ).
3. Write a Python function [gn, rn] = RDFcalc(S, dk, rho) to calculate Radial Distribution Function from Structure Factor
data, sampled at , and density , and output the RDF vector and its corresponding radial distance vector .
Hint: for Python fft() and ifft() functions, realize that the values of the Fourier Transform corresponding to negative
frequencies are stored in the second half of the arrays given to ( ifft ) or obtained from it ( fft ). You also have to study the
difference between the DFT and FT to multiply the right factors.
4. With the data provided in argon.py , compute the corresponding Radial Distribution Function . Plot your results for from to
Å. Over what range of radius can you trust your result?
Hint: To check if your results make sense, recall that is related to the probability that a molecule has a neighbouring molecule at
distance , therefore, should be close to when , i.e. two molecules can not occupy the same space, and you can set
. Recall . Also note the unit used in formula (6).
5. From the you computed, estimate the average molecular radius of liquid argon. Give your reasoning and state what accuracy
you can justify for your estimation.
6. Now we explore the effect of sampling range. Yan sampled in wavenumber out to Å , and he could have saved
himself work by not collecting as much data, i.e., reducing . But how much could he have reduced the sampling length ,
while still seeing distinct peaks in the Radial Distribution Function? Also explain theoretically what you observe.
Hint: Plot on top of the obtained in part 4, the 's you compute for a series of values. You can try half each time to
look for changes. For the theoretical explanations for part 6 and 7, realize the interchangability of ( ).
7. To explore the effect of data sampling, let's assume Yan decided to save his work by sampling less often (i.e. increasing ). How
large a can he use to be able to still recover the first two peaks clearly? State your answers and a theoretical justification for what
you expect to see if you increase too much.
Hint: Plot on top of the obtained from argon.py data, the 's you obtain when you subsample the same dataset. Try
doubling each time to observe the effect of coarser sampling.