Get documents about "Whats the deal with FFTs
Document Sample
What’s the deal with FFTs?
Occasionally in physics it is convenient to decompose functions into a sum of sin and cos
components. This is particularly useful when there is a periodic solution, like a wave on
a string where the solution may be sin(n*2x/L) and n is the harmonic of the wave.
Wave solutions are common in quantum mechanics because of the wave functions, as
well as in general harmonic oscillator solutions. In fact, wave packets often have terms
e2ft in them which is the main ingredient of the Fourier Transform. The definition of a
fourier transform is:
H ( f ) h(t )e 2ift dt OR H ( f ) h( x)e 2ix / dx ,
depending on whether you want to transform a function of time, h(t) or position, h(x).
The inverse transform is defined to be:
h(t ) H ( f )e 2ift
dt OR h( x) H ( f )e 2ix / dx .
So, we can go back and forth between the frequency domain and the time or position
domain as needed.
But, why bother with the frequency domain at all when the computer is there to do the
work for us? The basic reason is that there are several integrals that can be evaluated far
more efficiently in the frequency domain than in the time or position domain. Today
we’ll explore the convolution integral.
Convolution Integral
The convolution of g(x) and h(x) is denoted g*h and is defined by
g h g ( )h( x )d
where is the distance that is integrated over to leave a function of x. For example, if
g(x)=sin(x) and h(x)=cos(x) then the convolution is g h sin( ) cos(x )d .
500
400
300
200
100
g*h
0
-100
-200
-300
-400
-500
0 2 4 6 8 10 12 14
x
The cool thing is though, in Fourier space the integral is a simple multiplication.
g h G( f ) H ( f ) “Convolution Theorem”
There are many situations where a convolution that takes more than half an hour using
loops can be executed in less than a minute in Fourier space.
Detector Response:
Whenever we measure a signal what we see is the signal smeared out by the detector
response. Let’s look at a detector with a simple RC response =5 seconds (shown in the
left figure below). If we measure a square signal (blue curve on the right) with this
response we will find that the signal is blurred a bit by the detector response (green curve
on the right). This isn’t too much of a problem if the detector response is a lot shorter
than the signal, but we often find ourselves probing the limits of what our detector can
measure. In these cases its useful to use computer simulations of the response convolved
with the signal to understand how the blurring affects the science measurement.
detector response true and measured signal
0.8
1
0.6
response
signal
0.4
0.5
0.2
0 0
0 50 100 0 50 100
time (s) time (s)
What really happens when we convolve a signal with a detector response? You can think
of it as taking a single bin, signal(i) and putting parts of it into other bins based on the
shape of the response. In the following figure, each point in the “square signal” is
distributed into additional bins based on the shape of the response. These blurred
contributions are shown in the dotted curves. The sum of the dotted lines is then
displayed in green showing our simulated measurement again.
1
signal decomposition
0.8
0.6
0.4
0.2
0
0 10 20 30 40 50 60 70 80 90 100
time (s)
So, how do we do this using matlab? Lets look at a signal with a Gaussian shape.
% gaussian signal
t = 0:0.01:100;
tz = 50; sigma = 2;
signal = 1/sqrt(2*pi)/sigma*exp(-(t-tz).^2/sigma^2);
% detector response
tau = 5;
response = exp(-t/tau);
response = response/sum(response); % normalize to 1.
% convolution
measSig = convn(signal,response);
plot(t,signal, t,measSig(1:length(t)));
The blue curve is the original signal and the green curve is the signal convolved with the
detector response. Notice that the convolution is performed using Matlab’s convn
function. This function does the following:
% convolution using fft directly
kSig = fft([signal zeros(size(response))]); % signal padded with 0
kRes = fft([response zeros(size(signal))]); % response padded with 0
kMeas = kSig.*kRes; % multiply in frequency space
measSig2 = ifft(kMeas); % inverse FFT to get back to time space
plot(t,signal, t,measSig2(1:length(t)))
This version of convolution is performed using the fft of the signal and resolution,
multiplying the coefficients in frequency space and then transforming the result back into
the time domain with an inverse fft called ifft. When it is plotted you’ll see they give
exactly the same result.
M14_FFT1 Assignment
1. Convolve a signal that looks like the error function with the detector response above.
erfSig = 0.5*(erf(t-tnot) + 1);
How would you determine if the measured signal was the expected “erf” function or
some other shape? If you were using this detector to measure the time when the signal
transitioned from 0 to 1, how would you do this?
2. Cameras and telescopes measure images. Image resolution comes from a variety of
sources that may include: atmosphere distortion, diffraction from the edge of the optical
aperture, or CCD electronic resolution. Lets model all these effects as a Gaussian:
[xg,yg] = meshgrid(-20:20, -20:20);
sigma = 2;
response = exp(-(xg.^2 + yg.^2)/sigma);
response = response/sum(sum(response));
surf(xg,yg,response);
Pick your favorite image of the Mandelbrot set and convolve it with Gaussian resolution.
Include a surface plot of the Gaussian surface in your write-up. Include the original and
convolved image. Describe the science impact of loosing detail due to detector
resolution. Try it again with sigma =0.5 and sigma = 10.
3. (optional) Download a jpg image from the internet and convolve it with a 2D
resolution. You’ll want to read the image into matlab pixels using imread. Most images
that you download will have an array for each color, [red,green, blue]. You’ll have to
smear each color independently. One of my favorite images is the Hubble Deep Ultra-
deep field, hubble ultra-deep field. Its about 7 MB and takes awhile to convolve. You
may find it amusing to convolve your image with something more complicated than a
single symmetric Gaussian. You’ll notice all sorts of problems with the colors when you
really get to changing them. These are problems that Kodac generally corrects when it
prints photos. You can color correct Matlab pictures using the histeq function. This
makes pictures look good but they will no longer be radiometrically accurate.
