Plotting Fourier spectrum versus frequency of a signal

Question

I have looked around here, and i am sure this has been answered, but i don't understand it. The thing is, I have taken a introductory signal processing course, and we had to use Mathematica, and i had no previous experience in Mathematica, so this made a lot of things hard for me.

Now i need to understand how to visually present the spectral content of a sampled signal. So i try something simple,

sampls = Table[Sin[n*(2 \[Pi])/1000*4], {n, 0, 2000}];

where we can assume that the sampling period is 1/1000s and we have a sinus with a frequency of 4 Hz. This signal is sampled for 2 seconds, meaning we get 8 periods of the signal.

Now my naive attempt is to do this (as i saw something like this in some example somewhere):

 ListLinePlot[Abs[FourierDCT[sampls]], PlotRange -> {{0, 100}, {0, 30}}]

Now i end up with this:

which i can't make any sense out of...

What does it mean?

My dream is to be able to plot frequency on the x-axis and amplitude on the y-axis, so that i could present the teacher with a nice narrow peak with height 1 right above 4 on such a plot.

My mathematical knowledge is rather weak, i'm sorry to say, so be gentle with me and explain as concrete as possible. This is not primarily a mathematics course.

Answer 1

This is how I would do this. Define frequencies and sampling rate precisely. Then use Periodogram because it takes SampleRate as an option and rescales frequency axis automatically. Read up Docs on Periodogram - see examples there.

data = Table[{t, Sin[2 Pi 697 t] + Sin[2 Pi 1209 t]}, {t, 0., 0.1, 
    1/8000.}];

ListLinePlot[data, AspectRatio -> 1/4, Mesh -> All, 
 MeshStyle -> Directive[PointSize[Small], Red], Frame -> True]

Periodogram[data[[All, 2]], SampleRate -> 8000, Frame -> True, 
 GridLines -> {{697, 1209}, None}, 
 FrameTicks -> {{All, All}, {{697, 1209}, All}}, 
 GridLinesStyle -> Directive[Red, Dashed], AspectRatio -> 1/4]

Your case specifically is also very simple - sampling rate 1 - or Automatic - and frequency max at 4/1000:

sampls = Table[Sin[n*(2 \[Pi])/1000*4], {n, 0, 2000}];

Periodogram[sampls, Frame -> True, GridLines -> {{4/1000.}, None}, 
 GridLinesStyle -> Directive[Red, Dashed], 
 PlotRange -> {{0, .05}, All}]

Answer 2

I like @Vitaliy's answer, but here's another approach using Fourierinstead of Periodogram.

time = 2;
tinc = 0.001;
sampls = Table[Sin[n*(2 Pi) 4], {n, 0, 2, tinc}];
nyquist = 1/(2 tinc)
len = Length@sampls;
ListLinePlot[Sqrt[4/len] Abs@Fourier[sampls], 
 PlotRange -> {{0, 10}, All}, DataRange -> {0, (len - 1)/time}]

Briefly, I construct a sample set by assuming the units of the Table iterator are seconds and the sampling rate is defined by the increment step (1/tinc). We should always keep in mind the Nyquist frequency so we know the upper limit of frequencies to be determined by the transform. I find DataRange helpful as it avoids having to create multidimensional lists (which isn't a burden in this case) and shows how the x-axis of the transform is related to the original data. The amplitude is determined by the sample size and can be obtained by multiplying the transformed data by Sqrt[4/] as noted in the documentation for Fourier.

In addition to the documentation, I find this resource helpful when implementing FT in Mathematica.

Answer 3

Jenson is correct in his comment above ("Use FourierDST"). The reason is that your waveform is a sine wave - that is, it is zero at the ends of your time interval. All cosine harmonics are 1 at n=0, and for the waveform shown 1 at the end of the interval. Therefore, approximating a sine wave with a sum of cosines is not a trivial exercise, and you get a mess.

Technically speaking, the problem is that your "basis" functions are poorly chosen if you try to use cosines (even harmonics with respect to reflection about 0) to approximate sines (odd harmonics with respect to 0)

Your transform result has two large-ish cosine waves with a difference frequency of 2, i.e. approximately cos(F+delta)t + cos(F-delta)t. Their sum gives cos(F t)cos(delta t)-sin(F t)sin(delta t) + cos(F t)cos(delta t)+sin(F t)sin(delta t), or 2*cos(F t)cos(delta t), although this is only approximate since the two amplitudes are not quite equal. Off the top of my head, I think delta (= 1) gives one period of amplitude modulation over the time window, so this pair of cosine waveforms sums to a waveform that is -1 at the ends of the time interval, and has a +1 amplitude hump in the middle - that is, a modulated cosine wave.

The next two farther apart peaks give (approximately) a modulated cosine that has 2 humps, and so forth. Because they are cosines, none of them are zero at the ends of the interval, although summing enough of them probably gets close to zero. At any rate, summing up enough of them eventually can give an approximation to the sine wave.

Answer 4

I'm not big on Mathematica, but... your plot looks familiar.

The double-spike in the middle happens when your sinewave frequency isn't some nice multiple of the sample rate. As a result, the "energy" gets "smeared" across the frequency spectrum.

When doing a Fourier Transform, this also happens when sampling signals that suddenly "start" and "stop". The classic way to treat this is to apply a "window function". In signal processing, these have names like Rectangle, Triangle, Hamming / Hanning, and Blackman-Harris. They basically attenuate, or taper, the signals at the start and finish of your sample period, while leaving the ones in the middle alone.

The shape of the window can linear (Triangle), square root, power, or cosine. Check the Wikipedia page for some example pictures: http://en.wikipedia.org/wiki/Window_function