







CTK Exchange
alan
guest

Mar1006, 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:(N1))/(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:(N1) produces an array of N points starting at 0). The resulting plot is here: Any help would be appreciated. 

Alert  IP 
Printerfriendly page 
Reply 
Reply With Quote  Top 


alexb
Charter Member
1800 posts 
Mar1106, 10:21 AM (EST) 

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 
Printerfriendly page 
Reply 
Reply With Quote  Top 



alan
guest

Mar1106, 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 100050, 100075, 1000125. 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 
Printerfriendly page 
Reply 
Reply With Quote  Top 




alexb
Charter Member
1800 posts 
Mar1106, 03:58 PM (EST) 

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 
Printerfriendly page 
Reply 
Reply With Quote  Top 




mr_homm
Member since May2205

Mar1106, 11:17 PM (EST) 

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 Fourierlike) 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 N1. 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 
Printerfriendly page 
Reply 
Reply With Quote  Top 






mr_homm
Member since May2205

Mar1306, 08:03 AM (EST) 

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 N1." 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 
Printerfriendly page 
Reply 
Reply With Quote  Top 




alan
guest

Mar1206, 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 
Printerfriendly page 
Reply 
Reply With Quote  Top 



You may be curious to have a look at the old CTK Exchange archive. Please do not post there.
Front page
Contents
Copyright © 19962018 Alexander Bogomolny

