For each of the following exercises create the specified files with your answer(s). Submit your homework using the course submit utility.
This assignment carries out a Fourier transform of real-valued data. There are two applications. The first produces a frequency power spectrum (one-sided power spectral density). The second is used to smooth a signal by simulating the effect of a low pass filter. The data files are in ~p6730/examples/a05/data.
| File | Purpose of procedure |
|---|---|
| spectrum.cc | Driver (main) program. (incomplete) |
| Makefile | Rule for compiling the above. (provided) |
This code calls the C++-compatible Numerical Recipes routine realft.cc for doing the Fourier transform.
All calculations are to be done in double precision.
Takes time series data line-by-line from standard input until an EOF is reached. Each line has two values
t y
where t is the time and y is the amplitude to be
transformed. It computes the correct number of points for the FT,
figures out the time interval dt, pads the data with zeros if
necessary, does the FT and writes out the frequency power spectrum
(one-sided power spectral density) to standard output as a list of the
form
freq power
The frequency power spectrum is essentially the square of the complex
modulus as a function of frequency divided by the number of signal
points (including zero padding). It is normalized so the sum of the
power values equals the sum of the squares of the signal values.
All output (if any) except for the power spectrum is written to standard error. In this way it will be possible to pipe the power spectrum into a plotting utility.
The subprogram that performs the Fourier transform is called realft. This routine in turn calls four1. Please see ~p6730/examples/nmr. The header files are in ~p6730/nr/recipes_cpp/utils.
Here is a description of the calling sequence for realft. We must start with a time series with equal time spacing. Call this dt. Our FFT requires that number of points must be an integer power of 2. That means that the number of points must be 1 or 2 or 4 or 8 or 16... If not, it is necessary to pad the end of the data array with zeros. We then call realft to get the Fourier transform.
C++ linkage
#include <nr.h>
Vec_DP f(n);
int isign;
NR::realft(f, isign);
The arguments have the following meaning:
The output vector from realft contains the Fourier transform, as described above. Special course utilities setft and getft are provided for decoding it. The call
#include <nr.h>
#include <complex.h>
double_complex getft(Vec_I_DP &ft,int k);
Vec_DP f(MAX);
int k;
double_complex ft;
ft = getft(f, k)
returns a complex number
equal to the Fourier transform
at frequency index k.
The call
#include <nr.h>
#include <complex.h>
void setft(Vec_O_DP &ft, int k, double_complex &val);
Vec_DP f(MAX);
int k;
double_complex ft;
setft(f, k, ft);
sets the value of the Fourier transform to ft at frequency
index k.
To generate the time series, use the following command
~p6730/exercises/ode/anhrk4 < input | awk '{print $1,$2}' | tail -8192
where input is either data_chaos or data_regular from
~p6730/examples/a05. You are welcome to experiment with other
input parameters. Here the only difference between these data sets is
the strength g of the harmonic driving term. The awk
and tail commands select the time and angular displacement
values and keep only the last 213 = 8192 data points so any
transients will have died out.
For this exercise, feed the time series for both cases into your spectrum code and plot the resulting spectrum for frequencies on the interval [0,1]. (It might help to make a plot of the logarithm of the spectral density vs frequency.) Put this portion of the spectrum in the files chaos.txt and regular.txt for handing in.
In a file spectrum.txt, also to be submitted, answer the following questions:
| File | Purpose of procedure |
|---|---|
| ftfilter.cc | Driver (main) program. (incomplete) |
| Makefile | Rule for compiling the above. (provided) |
The Fourier transform of the RC filter function is APPROXIMATELY
1/(1-j*2*PI*freq*tau)
where j = sqrt(-1), PI = 3.1415..., and tau is
the time constant for the filter. See the lab exercise for
a more accurate expression. Take your choice for this exercise!
This program first prompts for and reads the name of the input file. Next it prompts for and reads the value of the time constant tau. After that the computation proceeds as in your program spectrum, except that it doesn't compute and report the power spectrum. Instead, it multiplies the FT by the FT of the filter function, as given above, and then retransforms. (Note that at the highest frequency, you will multiply only by the real part of the filter function transform). The filtered time series data should be written to standard output. All other output goes to standard error.
Additional comments:
Please note that you should be careful to pad the input signal with zeros over a time interval sufficiently long that you minimize "wrap-around pollution". For the RC filter let's choose a padding time interval equal to 4*tau. (If that padding doesn't make the number of data points a power of 2, you will have to pad further to make it a power of 2.)
Run your program with the test files test1 and test2 and with a variety of time constants, tau = 0.1, 0.5, ..., as you wish. One data file is a square wave pulse. The other is a "delta" function spike. Notice that the time interval drifts a bit in this data, but ignore the drift. Plot the filtered signal to see the effect of the filter. In a file filter.txt, to be submitted, explain your observations.
Run your program with the file pulse1 and a time constant tau = 0.1.. This data represents the photocathode current in a photomultiplier tube, while it was recording a light pulse. It was recorded by the Utah cosmic ray group. The first value on each line is the time and the second is the signal amplitude. Notice that the time intervals between adjacent measurements are not precisely the same, but drift a little. So we take the average time interval.
The pulse timing, pulse height, and pulse area are important for detecting cosmic rays. The detailed wiggles result from amplifier noise and are not interesting. Smoothing the data is useful as a preliminary step in determining the useful parameters. The filter performs the smoothing. Plot your result to see the effect of the filter. Put your output in a file pulse1_filter.txt for submitting.
Using gnuplot, create a plot with both the original and filtered data displayed. Create a gnuplot command file pulse.gpl for this purpose and hand it in. Construct your gnuplot command file so it builds the plot from the original data file and the result file pulse1_filter.txt so we will be able to reconstruct your plot.