Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

Plot of samples

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:

output from above statement

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.

share|improve this question
1  
For startes you should use FourierDST - the discrete sine transform. But the peak's position still depends on the sample size. –  jenson Apr 7 at 9:17
1  
Why should i use the sine tranform rather than the cosine? Could you explain a bit further? –  mickey Apr 7 at 10:52
    
Yuo need to add the corresponding frequency range. Have a look at my code in mathematica.stackexchange.com/questions/18082/… (Oh, and I use Fourier) –  Peltio Apr 7 at 16:52
    
Also see an detailed explanation of the frequency arrangement in the output of Fourier here. –  xslittlegrass Apr 7 at 17:20
    
If you want to call somebody in the comment, you need to write something like @jenson . BTW, I'm also interested in an explanation for the bad result of FourierDCT. (I roughly know that it's because Sin[x] is a odd function, but don't know how to explain it precisely. ) –  xzczd Apr 8 at 3:13

4 Answers 4

up vote 20 down vote accepted

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]

enter image description here

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}]
share|improve this answer
    
+1 In my Mma version an Integer is required in the SampleRate option –  belisarius Apr 7 at 11:37
    
@belisarius thanks, very sharp eye, - corrected to integer ;-) –  Vitaliy Kaurov Apr 7 at 11:40

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}]

Mathematica graphics

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.

share|improve this answer

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.

share|improve this answer
    
but in the general case, how would i know which version maps best to any given signal? Or is it as others have proposed always like this since we use windows and the windowed signal are always 0 at it's ends? –  mickey Apr 8 at 9:37

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

share|improve this answer
1  
Er……can you show an example fixing OP's code with the method you suggested here? –  xzczd Apr 8 at 3:03
1  
Sorry - like I said, I'm not familiar with Mathematica. However... I'd suggest changing the sample size to 2048: Fast Fourier Transforms in particular prefer multiples of 2 as sample size. Then I'd try a simple [triangle] window: OUT = Data * X / 1024 for X = points 0 to 1023, OUT = Data * (1-X) for points X = 1024 to 2047. –  Alan Campbell Apr 13 at 3:54
    
Your suggestion for changing sample size to 2048 and using a triangle does work (if I've understood it correct). Maybe you can consider improving your answer when you become more familiar with Mathematica, most people here are not familiar with signal processing after all :D –  xzczd Apr 14 at 2:39

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.