CTK Exchange
Front Page
Movie shortcuts
Personal info
Awards
Reciprocal links
Terms of use
Privacy Policy

Interactive Activities

Cut The Knot!
MSET99 Talk
Games & Puzzles
Arithmetic/Algebra
Geometry
Probability
Eye Opener
Analog Gadgets
Inventor's Paradox
Did you know?...
Proofs
Math as Language
Things Impossible
My Logo
Math Poll
Other Math sit's
Guest book
News sit's

Recommend this site

Manifesto: what CTK is about Search CTK Buying a book is a commitment to learning Table of content Things you can find on CTK Chronology of updates Email to Cut The Knot Recommend this page

CTK Exchange

Subject: "spectral analysis question"     Previous Topic | Next Topic
Printer-friendly copy     Email this topic to a friend    
Conferences The CTK Exchange College math Topic #562
Reading Topic #562
alan
guest
Mar-10-06, 06:04 PM (EST)
 
"spectral analysis question"
 
   I'm trying to teach myself a little bit about fourier series and spectral analysis. I'm using matlab to create a plot of the power spectrum of a simple periodic signal, but I don't quite understand why the plot looks like it does. The signal is y=sin(50*2*pi*t)+sin(75*2*pi*t)+sin(125*2*pi*t). The power spectrum has spikes at 50, 75, 125 as expected, but there are also spikes at 375, 425, 450. I'm using the matlab function FFT(x,N) which computes a discrete fourier transform of the signal x using N points. I'm using N=500, so obviously the three extra spikes are related to that number, but I don't understand exactly how. My guess is it has something to do with the corresponding negative frequencies, but I don't understand how this implementation of the DFT is treating those frequencies.

The code I'm using:

function spec()
dt=.002;
N=500;
f1=50;
f2=75;
f3=125;

t = 0:dt:0.6;
y = sin(f1*2*pi*t)+sin(f2*2*pi*t)+sin(f3*2*pi*t);
Y = fft(y,N);
Pyy = Y.* conj(Y)/N;
f = (0:(N-1))/(N*dt);
plot(f,Pyy)

The code is pretty simple, it'should be easy to understand for anyone familiar with programming (just in case, conj(x) is the complex conjugate of x, and the notation 0:(N-1) produces an array of N points starting at 0). The resulting plot is here:

Any help would be appreciated.


  Alert | IP Printer-friendly page | Reply | Reply With Quote | Top

  Subject     Author     Message Date     ID  
spectral analysis question alan Mar-10-06 TOP
  RE: spectral analysis question alexb Mar-11-06 1
     RE: spectral analysis question alan Mar-11-06 2
         RE: spectral analysis question alexb Mar-11-06 3
         RE: spectral analysis question mr_homm Mar-11-06 4
             RE: spectral analysis question alexb Mar-12-06 5
                 RE: spectral analysis question alexb Mar-12-06 6
                     RE: spectral analysis question mr_homm Mar-13-06 8
             RE: spectral analysis question alan Mar-12-06 7

Conferences | Forums | Topics | Previous Topic | Next Topic
alexb
Charter Member
1800 posts
Mar-11-06, 10:21 AM (EST)
Click to EMail alexb Click to send private message to alexb Click to view user profileClick to add this user to your buddy list  
1. "RE: spectral analysis question"
In response to message #0
 
   Have you tried to plot the function itself?

The problem my be in what Matlab sees when you tell it the function is a sine. There are at least three ways to extend the sine defined in half the usual period for positive t to the negative domain: making it either even or odd or zero. You'll get different spikes in the Fourier transform for each of the three.

I have a hunch that you problem has something to do with setting

dt = 0.002.

Have you tried dt = 0.001?

Also, why did not you take N = 512?



  Alert | IP Printer-friendly page | Reply | Reply With Quote | Top
alan
guest
Mar-11-06, 12:44 PM (EST)
 
2. "RE: spectral analysis question"
In response to message #1
 
   The function doesn't seem to care whether N is a power of two or not, and when I change it to 512 there isn't any discernible difference. Changing dt to .001 just changes the scale of the frequency axis, and the three spikes go to 1000-50, 1000-75, 1000-125. The signal function itself looks as I expect it to look.
For a signal sin(f*t), this algorithm should produce a spike at +f and -f, correct? I think those three spikes must be the negative frequencies, plotted in the wrong place somehow. Either that or a sampling artifact, but they remain no matter what the sample rate is.


  Alert | IP Printer-friendly page | Reply | Reply With Quote | Top
alexb
Charter Member
1800 posts
Mar-11-06, 03:58 PM (EST)
Click to EMail alexb Click to send private message to alexb Click to view user profileClick to add this user to your buddy list  
3. "RE: spectral analysis question"
In response to message #2
 
   There is a certain idiosynchrasy in the Matlab implementation which you should learn about. I just looked at their manual at

https://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.html

It'shows an example with two frequencies and random noise, but otherwise, is very much like yours:

==== the last part of the example ====

The power spectrum, a measurement of the power at various frequencies, is

Pyy = Y.* conj(Y) / 512;

Graph the first 257 points (the other 255 points are redundant) on a meaningful frequency axis:

f = 1000*(0:256)/512;
plot(f,Pyy(1:257))
title('Frequency content of y')
xlabel('frequency (Hz)')

=====

Just pay attention to how they execute the plotting. And why half the points are redundant, or what they mean by that


  Alert | IP Printer-friendly page | Reply | Reply With Quote | Top
mr_homm
Member since May-22-05
Mar-11-06, 11:17 PM (EST)
Click to EMail mr_homm Click to send private message to mr_homm Click to view user profileClick to add this user to your buddy list  
4. "RE: spectral analysis question"
In response to message #2
 
   What you are seeing is not an artifact of the sampling size, nor is it a peculiarity of MatLab's implementation of the discrete Fourier transform. Rather, it is an inherent feature of the DFT itself, regardless of the particular algorithmic implementation used.

Here are the relevant properties of Fourier (and Fourier-like) transforms in general:

Transforms of periodic functions result in discrete functions (this is fairly obvious, since if a periodic function can be decomposed into a sum of sinusoids at all, it must be decomposed into a sum containing only frequencies that are multiples of the original function's frequency).

Transforms of aperiodic functions result in continuous functions (or actually functions that are nonzero on a dense set, of which continuous functions are a subclass).

Inverse transforms exist, and they have the same mathematical form as the forward transforms, except for some sign or scaling constant differences.

Based on these facts, for any pair of functions that are transforms of each other, if one is periodic, the other is discrete, and if one is aperiodic, the other is continuous, and this works both ways. Therefore, there are 4 distinct versions of the Fourier transform, based on whether the functions are discrete (D), continuous (C), periodic (P) or aperiodic (A):


CA ---Fourier Transform------------------>CA
CP ---Fourier Series--------------------> DA
DA ---Discrete Time Fourier Transform---> CP
DP ---Discrete Fourier Transform--------> DP

Note that the DFT is inherently periodic, so the complete information content of the DFT is displayed within any given single period. It is your choice which period to display. However, since the DFT is a complex Fourier transform, trig functions such as sine and cosine will of course have amplitude peaks at both positive and negative frequencies in their DFT, and because they are real functions, these complex amplitudes will be conjugates of each other.

Therefore, the most natural interval on which to display the DFT is one that is symmetrical around zero. For an N point DFT, this would range from -N/2 up to N/2 - 1, assuming N is even. However, because of the way MatLab is set up, it likes to use vectors that are indexed by nonnegative integers, so the period that is chosen is from 0 to N-1. Hence the information that is located at negative frequencies will be displayed in the upper half of the spectrum. Just imagine that you cut off the upper half of the spectrum, move it over, and reattach it at the left side of the lower half of the spectrum. I think that answers your question about the graph.

There are a couple of other points as well. You probably noticed that the spikes are not perfectly sharp, but frequencies near 50, 75, and 125 also have nonzero amplitudes. This is NOT an inherent property of the DFT. As in a true Fourier transform, these peaks should be delta functions (all value concentrated at a single point), but the DFT in MatLab shows some widening of the peaks. This is an artifact of the calculation, as you can verify by computing the DFT of a pure sine wave by hand (an easy computation if you visualize the terms in the sum as vectors in the complex plane).

Also, as to the selection of the number of sample points, there is no need to make N a power of 2, because there are versions of the Fast Fourier Transform algorithm that can handle other values for N. Basically, the FFT algorithm is so fast because it divides the points into subsets and then recursively calls itself on these smaller sets. It is easy to write an FFT algorithm that divides the points into two subsets, but it is perfectly possible (and has been done) for arbitrary partitions into equal size subsets. There is an nice implementation of this version of the FFT in C language source code, free for noncommercial use. I forget its name, but if you do a Google search for the phrase "fastest FFT in the west" you will find it, as that's the slogan on the website.

It is still true that the fastest values of N are the powers of 2, but basically any N that factors into several small factors will be nearly as fast. What you want to avoid at all costs is a large prime number value for N. For example, 499 is prime, and will take several times as long to compute as 500, which factors as 2*2*5*5*5.

(My job is to coach students through this kind of stuff at the university, so if you have any further questions, please ask, it's no trouble.)

Hope this helps,

--Stuart Anderson


  Alert | IP Printer-friendly page | Reply | Reply With Quote | Top
alexb
Charter Member
1800 posts
Mar-12-06, 05:53 PM (EST)
Click to EMail alexb Click to send private message to alexb Click to view user profileClick to add this user to your buddy list  
5. "RE: spectral analysis question"
In response to message #4
 
   1. Why the manual example does not show the phantom frequencies?

2. You wrote "However, because of the way MatLab is set up, it likes to use vectors that are indexed by nonnegative integers, so the period that is chosen is from 0 to N-1." If it's not an idiosyncrasy then what is?


  Alert | IP Printer-friendly page | Reply | Reply With Quote | Top
alexb
Charter Member
1800 posts
Mar-12-06, 05:57 PM (EST)
Click to EMail alexb Click to send private message to alexb Click to view user profileClick to add this user to your buddy list  
6. "RE: spectral analysis question"
In response to message #5
 
   >1. Why the manual example does not show the phantom
>frequencies?

Ah, I understand. They dropped the right part of the positive interval because it would show the phantom frequencies inherited from the negative domain.

>2. You wrote "However, because of the way MatLab is set up,
>it likes to use vectors that are indexed by nonnegative
>integers, so the period that is chosen is from 0 to N-1." If
>it's not an idiosyncrasy then what is?


  Alert | IP Printer-friendly page | Reply | Reply With Quote | Top
mr_homm
Member since May-22-05
Mar-13-06, 08:03 AM (EST)
Click to EMail mr_homm Click to send private message to mr_homm Click to view user profileClick to add this user to your buddy list  
8. "RE: spectral analysis question"
In response to message #6
 
   >>1. Why the manual example does not show the phantom
>>frequencies?
>
>Ah, I understand. They dropped the right part of the
>positive interval because it would show the phantom
>frequencies inherited from the negative domain.
>
Yes, although these frequencies are not always redundant. If the original function is complex rather than purely real, then half the information content of the transform is in the negative frequencies. Of course, in most practical cases, the function is some physical signal that is a function of time, and so it is real. In that case, the negative amplitudes are the conjugates of the positive ones, and so they can be constructed (if wanted) from the positive ones. The manual seems to assume that this will always be the case, and that therefore the negative frequencies will be redundant. The manual should not make these assumptions, of course.

>>2. You wrote "However, because of the way MatLab is set up,
>>it likes to use vectors that are indexed by nonnegative
>>integers, so the period that is chosen is from 0 to N-1." If
>>it's not an idiosyncrasy then what is?

Well, it is and it isn't. Since the transform is periodic, one must choose some particular period to display, and the choice is essentially arbitrary. I would prefer to see an interval that is symmetric around the origin, as that makes the interpretation easier, but that's just a matter of preference. It is I think equally reasonable to take a period that starts at zero. Where the manual fails is in not documenting that there is a choice being made at all.

By the way, I did not mean to say that there were no idiosyncracies in MatLab; it is full of them, which steepens its learning curve unnecessarily. What I meant was that there was nothing unusual about its implementation of the FFT algorithm itself. How to display the results of the calculation is a separate issue, and here I agree with you -- I would have preferred a symmetric interval.

--Stuart Anderson


  Alert | IP Printer-friendly page | Reply | Reply With Quote | Top
alan
guest
Mar-12-06, 09:57 PM (EST)
 
7. "RE: spectral analysis question"
In response to message #4
 
   That was very helpful, thank you. I knew that the transform of a periodic signal would produce two delta functions at the positive and negative frequencies, but I didn't make the connection between that, and Matlab not using negative indices. I did see the statement in the manual about the redundancy, but I wanted to see for myself what the full plot of the output of the function.


  Alert | IP Printer-friendly page | Reply | Reply With Quote | Top

Conferences | Forums | Topics | Previous Topic | Next Topic

You may be curious to have a look at the old CTK Exchange archive.
Please do not post there.

|Front page| |Contents|

Copyright © 1996-2018 Alexander Bogomolny

Search:
Keywords:

Google
Web CTK