Matlab activity 7: Sliding FFTs

Our previous MATLAB activities have dealt primarily with the study of periodic signals. The theory of Fourier series tells us that periodic signals can be represented as sums of the simplest kinds of periodic signals — sines and cosines.

This makes less sense for aperiodic signals, the signals we’ve began studying last week. To see this, let’s look at an example:

t=1/4096:1/4096:1;
f1 = 300;
x=cos(2*pi*f1*t.^2);
plot(t,x)

(Note the t.^2 inside the cosine.) A cosine wave of (constant) frequency f is given described by the equation x=\cos(2\pi\cdot f\cdot t). The function x=\cos(2\pi\cdot t\cdot t) is a cosine function of “variable” frequency. Near t, it’s frequency is roughly t. (Such a sound is called a “chirp”. There is actually a MATLAB command chirp that makes these kinds of signals.) Take a look at some “zoomed in” plots:

x1=x(2000:2199);
plot(x1);
hold on;
x2=x(3000:3199);
plot(x2,’red’);
x3=x(3897:4096);
plot(x3,’green’);

(Maybe try looking at them on separate plots too.) These look more like standard sinusoids. Approximately what are their frequencies?

If these zoomed in pieces of the signal x do indeed behave like sinusoids, then their FFTs should gives sharp spikes close to their frequencies. Do they?

X1=fft(x1);
X2=fft(x2);
X3=fft(x3);
hold off;
hold on;
plot(abs(X1))
plot(abs(X2),’red’);
plot(abs(X3),’green’);

Compare these plots with that of the FFT of the entire signal x:

hold off
X4=fft(x);
plot(abs(X4));

What is this plot trying to tell us? A lot of frequencies seem to be present. The spectrum is spread out between frequencies of 0 and like 900 cycle/second, where there is a steep dropoff. (This is sort of strange because our “variable frequency cosine wave” has frequency varying between 0 and 300 cycles/second.) From the plot in the time domain (the plot of x vs t), we saw that x oscillates very slowly near t=0 and that the rate of oscillation increases as t increases. This information cannot be read off this last plot. (At least I don’t see how to do it…)

The FFT is a great tool for reading off information about the frequency domain composition of signals so long as this composition isn’t varying too quickly. This means we should only perform FFTs on relatively short time intervals. Since our “variable frequency cosine wave” has its frequency varying between 0 and 30, and we know we should always sample at more than two times as great as the highest frequency component, let’s do FFTs of length 256=2^8. The simplest thing to do is to break our signal up into 16=length(x)/256 blocks of 256 samples, take the FFT of each, and “assemble the results”.

A = zeros(256,16);

for k = 1:16
xx = x((k-1)*256+1:k*256);
A(:,k)=fft(xx);
end

image(abs(A’))

The vertical axis along the bottom tells us which block’s FFT we’re looking at. The horizontal axis is frequency (scaled in a strange way, let’s not bother with that for the moment.) The color intensity indicates the strength (or amplitude or volume) of a given frequency within a given block. You should see a bright line sloping downwards from the top left and another one sloping downwards from the top right. The brightness indicates that our signal contains a dominant frequency component that is increasing as we move from block to block, i.e., as time progresses. This is the behavior we expect from our variable (increasing) frequency cosine wave. (The bright line sloping downwards from the top left comes from the symmetry of the Fourier transform. It gives no new information.)

This picture is still very rough, and it’s tough to read quantitative information off of it. First, the resolution is a little grainy because we only have 16 blocks along the horizontal axis. We can remedy this to an extent by allowing out blocks to overlap: Instead of our first FFT being samples 1 to 256, our second being 257 to 512, etc., we instead slide along our FFT window of length 256 by a smaller increment each time: Our first FFT would still be calculated from samples 1 to 256, but our second may be say from 128 to 384, and the third from 256 to 512, … I won’t make you code this process by hand. MATLAB has a built in function for exactly this: spectrogram. (The picture it gives you is also called a spectrogram.) Here’s the syntax:

spectrogram(x,256,128,256,4096)

(Note that the spectrogram has the axes interchanged: It has frequency on the horizontal axis and time on the vertical.) The first 256 is the “window” which for us at the moment is the same as the other 256, just the length of the FFT (the number of samples used to compute it). 128 is the overlap between each successive FFT. Thus, as stated there is 50% overlap between each FFT. 4096 is the sampling rate — just the length of x since our t ran from 0 to 1. Play around with the parameters and see what happens.

Try a different x:

x=cos(2*pi*f1*t.^3);
spectrogram(x,256,128,256,4096);

In the next weeks, we’ll touch on a few ways to optimize these spectrograms — make them smoother, have better resolution, less noise, etc. But for now I just wanted to introduce the notion. The spectrogram is one of the most useful and pervasive methods for visualizing frequency composition of aperiodic signals. The key is that it shows three parameters simultaneously — time, frequency (frequencies present in the signal at a given time) and intensity/amplitude/power of the signal.

The takeaway sentence from this exercise is: Since frequency composition of aperiodic signals is roughly constant over short enough time intervals, we can use our tools (FFT) for analyzing periodic signals — those whose frequency composition is constant over long intervals — in blocks for analyzing aperiodic signals.

This entry was posted in Uncategorized. Bookmark the permalink.

One Response to Matlab activity 7: Sliding FFTs

  1. Doug says:

    This article is very interesting, but is
    hard to find in google. I found it on 14 spot. You can reach google top-10 easily using one
    handy wordpress plugin and increase targeted traffic
    many times. Just search in google for:
    Aemikimi’s Rank Plugin

Leave a Reply

Your email address will not be published. Required fields are marked *