Steve on Image Processing

Previously in my Fourier transforms series I've talked about the continuous-time Fourier transform and the discrete-time Fourier transform. Today it's time to start talking about the relationship between these two.

Let's start with the idea of sampling a continuous-time signal, as shown in this graph:

Mathematically, the relationship between the discrete-time signal and the continuous-time signal is given by:

(When I write equations involving both continuous-time and discrete-time quantities, I will sometimes use a subscript "c" to distinguish them.)

The sampling frequency is (in Hz) or (in radians per second).

The discrete-time Fourier transform of is related to the continuous-time Fourier transform of as follows:

But what does that mean? There are two key pieces to this equation. The first is a scaling relationship between and : . This means that the sampling frequency in the continuous-time Fourier transform, , becomes the frequency in the discrete-time Fourier transform. The discrete-time frequency corresponds to half the sampling frequency, or .

The second key piece of the equation is that there are an infinite number of copies of spaced by .

Let's look at a graphical example. Suppose looks like this:

Note that equals zero for all frequencies . This is what we mean when we say a continuous-time signal is band-limited. The frequency is called the bandwidth of the signal.

The discrete-time Fourier transform of looks like this:

where . As I mentioned before, normally only one period of is shown:

For this example, then, between and looks just like a scaled version of .

Next time we'll consider what happens when doesn't look like . In other words, we're about to tackle aliasing.

Get the MATLAB code

Published with MATLAB® 7.9

Last fall Loren wrote a blog post about a new syntax in R2009b for ignoring function inputs or function outputs. For example, suppose you call sort and only want the second output. In R2009b you can do that using this syntax:

 [~, idx] = sort(A);

The tilde character here is used as a placeholder to indicate that you don't need the first output argument of sort.

For another example, suppose you write a function that has to take three input arguments (because another function is always going to pass it three arguments), but your function doesn't need the second argument. Then you can write the first line of your function this way:

 function out = myfun(A, ~, C)

This new syntax has drawn a startling amount of discussion, both in the comments on Loren's post as well as in the comp.soft-sys.matlab newsgroup. (See the thread "Getting indexes of rows of matrix with more than n repetitions," for example.)

Responses to this new syntax have fallen into roughly four categories:

The passionate arguments in comp.soft-sys.matlab caught my eye and have prompted me to add my two cents here. Although I have no particular expectation of changing anyone's mind, I thought it might be interesting to address one particular question that was expressed well by Matt Fig: "If something already works, why complicate things?"

Matt means that we already have at least a couple of ways to ignore an output variable. The first is to use a variable name that (hopefully) makes the programmer's intent clear:

  [unused, idx] = sort(A);

The second technique takes advantage of the way MATLAB assigns function outputs from left to right:

  [idx, idx] = sort(A);

So why go to the trouble of introducing a new syntax?

Well, you'll probably get a somewhat different answer from every MATLAB developer you ask. Some might even agree. (I'll note, however, that this was one of the least-controversial MATLAB syntax proposals ever considered by the MATLAB language team.)

The proposed syntax was originally considered years ago, sometime around 2001. The proposal was approved internally at that time but then didn't get implemented right away because of competing priorities.

To understand why we eventually did decide to go ahead and implement it, you have to understand how this fairly minor syntactic issue is connected to our long-term efforts to do something much more important: provide automatic (and helpful!) advice to MATLAB users about their code. If you've used the MATLAB Editor at all over the last few years, you've probably observed how it marks code sections and makes suggestions to you. These suggestions fall into several categories, such as errors, possible errors, performance, new features, etc.

One of the things we try to flag is potential programming errors. A coding pattern that very often indicates a programming error is when you save a computed value by assigning it to a variable, but then you never use that saved value. At the very least, that pattern may indicate "dead code," or code that was doing something useful at one time but is now just cruft.

Both of the conventions mentioned above ([unused,idx] = sort(A) and [idx,idx] = sort(A)) exhibit this pattern of computing and saving a value and then never using it.

That might not be obvious for the [idx,idx] = sort(A) case, but it's true. The first output of sort is assigned to the variable idx. Then the second output is assigned to idx, causing the first output value to be discarded without ever being used.

These cases aren't programming errors, though, because we've given you no other way to ignore outputs. We can't automatically and reliably distinguish between the intentional [idx,idx] = sort(A) and the similar-looking [y,y] = foobar(x) that is the result of a typo.

I think it's important for the programmer to communicate his or her intent very clearly (which is why I tend to prefer the [unused,idx] = sort(A) convention). It is useful to have a way to communicate intent syntactically instead of by convention. And the new syntax helps us, in a small way, progress toward our long-term goal of helping MATLAB users write better MATLAB code.

OK, fire away!

Get the MATLAB code

Published with MATLAB® 7.9

In the last two posts in my Fourier transform series I discussed the continuous-time Fourier transform. Today I want to start getting "discrete" by introducing the discrete-time Fourier transform (DTFT).

The DTFT is defined by this pair of transform equations:

Here x[n] is a discrete sequence defined for all n:

I am following the notational convention (see Oppenheim and Schafer, Discrete-Time Signal Processing) of using brackets to distinguish between a discrete sequence and a continuous-time function. n is unitless. The frequency-domain variable, , is continuous with units of radians.

Note that is periodic with period .

Here are a few common transform pairs:

Unit Impulse

DTFT of Unit Impulse

Rectangular Pulse

DTFT of Rectangular Pulse

Note that the DTFT of a rectangular pulse is similar to but not exactly a sinc function. It resembles the sinc function between and , but recall that is periodic, unlike the sinc function.

Cosine

DTFT of Cosine

The DTFT of a discrete cosine function is a periodic train of impulses:

I updated the above plot on 6-Jan-2010 to show the location of the impulses. -SE

Because of the periodicity of it is very common when plotting the DTFT to plot it over the range of just one period: . For example, the DTFT of the rectangular pulse will most often be shown like this:

Next time I'll discuss the relationship between the continuous-time and the discrete-time Fourier transforms. Until then, Happy New Year everyone!

Get the MATLAB code

Published with MATLAB® 7.9

Last week I showed a couple of continuous-time Fourier transform pairs (for a cosine and a rectangular pulse). Today I want to follow up by discussing one of the ways in which reality confounds our expectations and causes confusion. Specifically, when we're talking about real signals and systems, we never truly have an infinitely long signal.

To refresh your memory, here are the ideal cosine signal and its continuous-time Fourier transform plots again:

The dots at the left and right of the cosine plot are meant to remind you that the cosine signal is defined for all t.

But what if we modify the cosine signal so that we have only a finite number of periods? Let's start with a single period, for example.

% N = number of desired cosine periods
N = 1;

t1 = -N*pi;
t2 = N*pi;
t = linspace(t1, t2, 1000);
x = cos(t);
plot(t,x)
ylim([-1.25 1.25])
xlabel('t')
title('x(t) = cos(t)')

In signal processing jargon you'll see this kind of signal referred to as a windowed cosine with a rectangular window. That is, the signal is an infinite-extent cosine multiplied by a finite-extent rectangular pulse.

If I've done my math correctly, the continuous-time Fourier transform of a cosine signal truncated to a finite number of periods is:

where N is the number of periods.

(And if I haven't done my math correctly, I'm sure a sharp-eyed reader will let me know. I feel like I'm back at Georgia Tech working on homework sets. Except that I'm doing it during the day instead of late at night.)

omega = linspace(-5.5, 5.5, 500);
X = sin(pi*N*(omega-1))./(omega-1) + sin(pi*N*(omega+1))./(omega+1);
plot(omega, X)
title('X(\Omega)')
xlabel('\Omega (rad/s)')

Well, that doesn't look very much like the continuous-time Fourier transform for the ideal cosine! If you run this code yourself, try zooming in to look closely at those peaks. You'll see they aren't even located precisely at 1.0 rad/s.

Let's try 10 periods.

N = 10;
t1 = -N*pi;
t2 = N*pi;
t = linspace(t1, t2, 1000);
x = cos(t);
plot(t,x)
ylim([-1.25 1.25])
xlabel('t')
title('x(t) = cos(t)')

omega = linspace(-5.5, 5.5, 500);
X = sin(pi*N*(omega-1))./(omega-1) + sin(pi*N*(omega+1))./(omega+1);
plot(omega, X)
title('X(\Omega)')
xlabel('\Omega (rad/s)')

That's closer, but there's still a lot of funny stuff going on. The more periods we use, the closer the Fourier transform gets to the ideal case, but for any finite number of periods we're going to get this decaying ripple effect.

Later, when we start constructing synthetic spatial sinusoid images (like in this post), and what we see in the frequency domain looks "messy," this is one of the big reasons. All our images have finite extent! We'll be coming back to this idea in future posts.

Have you heard the term bandlimited? It refers to a signal that has no frequency energy above a certain cut-off frequency. Ponder this notion: No finite-extent signal is bandlimited. That is, any signal that is nonzero over only a finite portion of the time domain has nonzero frequency energy across the entire frequency domain.

I think that next time I'll be ready to start talking about the discrete-time Fourier transform, or DTFT. Conceptually we are traveling methodically toward the discrete Fourier transform, or DFT, which is what the MATLAB function fft computes. Along the way we'll figure out how all three forms (continuous-time Fourier transform, discrete-time Fourier transform, and discrete Fourier transform) relate to each other.

Get the MATLAB code

Published with MATLAB® 7.9

In my planned (well, partially planned) discussion on Fourier transforms, I'll focus on three of the four types I listed in my November 23 post:

The existence of multiple transform flavors, as well as the details of their relationships, is at the heart of much of the confusion on this topic.

Let's start with the continuous-time Fourier transform. (When the context makes it clear whether I'm talking about the continuous-time or the discrete-time flavor, I'll often just use the term Fourier transform.)

The continuous-time Fourier transform is defined by this pair of equations:

There are various issues of convention and notation in these equations:

All these variations exist because we don't already have enough to worry about in the rest of our lives.

With the equations I use, the frequency domain unit is angular frequency (radians/second).

There are so many useful Fourier transform properties and transform pairs that it's hard for to pick the bare minimum necessary for the ideas I want to convey. Today I'll just show you two of the most essential Fourier transform pairs in signal processing applications.

Here's a cosine signal:

Plots corrected December 14 thanks to help from Mark Andrews.

And here is its Fourier transform:

This is what most people who have some knowledge of the Fourier transform expect to see. A signal containing a single frequency (here the frequency is 1 rad/s) has all its frequency domain energy concentrated at that single frequency.

The second pair is a rectangular pulse in the time domain and a sinc function in the frequency domain.

I'd like to call your attention especially to the dots ("...") at the left and right on the cosine and sinc function plots. The dots are there to remind you that these functions have infinite extent. That is, they are nonzero over the entire domain.

That's an important thing to keep in mind. I'll come back to that point next time.

One final note: I've started using the category "Fourier transforms" for posts on this topic. You can see all the related posts by clicking on the category link on the right side of the page.

Get the MATLAB code

Published with MATLAB® 7.9