Help with signals stuff (FFT, PSD)

[KillerCharlie]

I basically need to analyze the frequency content of some functions - these are simple functions, not noisy signals or anything. I just want to see the frequency makeup of the functions - basically make a PSD plot. I'm not an EE so I'm having some trouble doing this.

Could anyone help me? I'm using Matlab. I've played around a bit using the fft function. It turns out okay on some functions but not others. Right now the function I want to look at is the normal distribution function. Later I might look at some other equations. Help would be much appreciated.

 

[jmcoreymv]

Generate a vector of whatever function you want to find the frequency components of. Then do fft(vector), then do plot(vector).

 

[PurdueRy]

Should be working fine, if its a nice looking function(not incredibly hard) you can't use the fourier transform equations and do it by hand

 

[KillerCharlie]

Originally posted by: jmcoreymv
Generate a vector of whatever function you want to find the frequency components of. Then do fft(vector), then do plot(vector).

You mean plot(fft(vector))? Say the function is X=sin(10*pi*t). There should be a spike at 10pi... but if I just plot the fft it doesn't make any sense.

 

[Stojakapimp]

i believe in the newer matlabs you use the pwelch command to plot PSD's. You'll want to store that in some variable like P=pwelch(data) or something and then plot that with the correct frequency vector. So I think you'll have to know what sampling rate was used for the signals.

 

[KillerCharlie]

Originally posted by: Stojakapimp
i believe in the newer matlabs you use the pwelch command to plot PSD's. You'll want to store that in some variable like P=pwelch(data) or something and then plot that with the correct frequency vector. So I think you'll have to know what sampling rate was used for the signals.

Sampling rate isn't a problem... I'm just using equations instead of taking actual data. If it makes any difference, here's the equation:

sigma=0.2;
X=(1/(sigma*sqrt(2*pi)))*exp(-(t).^2/(2*sigma^2));

 

[Stojakapimp]

Originally posted by: KillerCharlie

Originally posted by: jmcoreymv
Generate a vector of whatever function you want to find the frequency components of. Then do fft(vector), then do plot(vector).

You mean plot(fft(vector))? Say the function is X=sin(10*pi*t). There should be a spike at 10pi... but if I just plot the fft it doesn't make any sense.

Because if you do plot(fft(vector)) then you aren't using any sort of reference for the time axis. You'll see a spike somewhere, but your x-axis is simply incrementing up to the length of your vector. Therefore, you need to do plot(frequency_vector,fft(vector)), although I believe in actuality you want to use abs(fft(vector)) and then take only half of that vector, as I believe the second half is simply repetitive. I could be wrong though, so you're probably better off using the pwelch function

 

[WhoBeDaPlaya]

You want to be careful that you sample exactly one period of the source signal.
Taking anything more or less will give you an erroneously high noise floor and low SNR / SFDR.
Yeah, windows (Blackman, etc.) are supposed to alleviate the need for exact sampling, but they still don't provide stellar results.

 

[KillerCharlie]

Originally posted by: Stojakapimp

Originally posted by: KillerCharlie

Originally posted by: jmcoreymv
Generate a vector of whatever function you want to find the frequency components of. Then do fft(vector), then do plot(vector).

You mean plot(fft(vector))? Say the function is X=sin(10*pi*t). There should be a spike at 10pi... but if I just plot the fft it doesn't make any sense.

Because if you do plot(fft(vector)) then you aren't using any sort of reference for the time axis. You'll see a spike somewhere, but your x-axis is simply incrementing up to the length of your vector. Therefore, you need to do plot(frequency_vector,fft(vector)), although I believe in actuality you want to use abs(fft(vector)) and then take only half of that vector, as I believe the second half is simply repetitive. I could be wrong though, so you're probably better off using the pwelch function

I can get good results if I plot, say, sin(20*pi*t). However, if I plot something with a much lower frequency (like sin(pi*t) or sin(.5*pi*t)), it gets cut off in the lower end when I try to do the PSD plot - the spike is cut off on the left.

 

[Stojakapimp]

ok....let's do an example.

Say your signal is x=sin(10*pi*t). We know that this signal has a frequency of 5 Hz, since in general for sin(2*pi*t*w), w is your frequency in Hz. So now create some time variable t. So we have:

t=-10:.05:10;
x=sin(10*pi*t);

you can plot this signal if you want by doing plot(t,x). Anyways, to find the frequency of this signal (which we know is 5 Hz), you can use pwelch. Since in my example i am increasing my time interval by .05 seconds, that gives me a sampling frequency of 1/.05=20. So 20 is my sampling frequency. To use pwelch, I type:

pwelch(x,[],[],[],20)

This should give you a plot of the psd with a spike at 5 Hz.

 

[hypn0tik]

Originally posted by: KillerCharlie

Originally posted by: Stojakapimp
i believe in the newer matlabs you use the pwelch command to plot PSD's. You'll want to store that in some variable like P=pwelch(data) or something and then plot that with the correct frequency vector. So I think you'll have to know what sampling rate was used for the signals.

Sampling rate isn't a problem... I'm just using equations instead of taking actual data. If it makes any difference, here's the equation:

sigma=0.2;
X=(1/(sigma*sqrt(2*pi)))*exp(-(t).^2/(2*sigma^2));

That's a Gaussian distribution.

Mathworld

 

[KillerCharlie]

Originally posted by: hypn0tik

Originally posted by: KillerCharlie

Originally posted by: Stojakapimp
i believe in the newer matlabs you use the pwelch command to plot PSD's. You'll want to store that in some variable like P=pwelch(data) or something and then plot that with the correct frequency vector. So I think you'll have to know what sampling rate was used for the signals.

Sampling rate isn't a problem... I'm just using equations instead of taking actual data. If it makes any difference, here's the equation:

sigma=0.2;
X=(1/(sigma*sqrt(2*pi)))*exp(-(t).^2/(2*sigma^2));

That's a Gaussian distribution.

Mathworld

Yes... yes it is....

 

[Stojakapimp]

Originally posted by: KillerCharlie

Originally posted by: Stojakapimp

Originally posted by: KillerCharlie

Originally posted by: jmcoreymv
Generate a vector of whatever function you want to find the frequency components of. Then do fft(vector), then do plot(vector).

You mean plot(fft(vector))? Say the function is X=sin(10*pi*t). There should be a spike at 10pi... but if I just plot the fft it doesn't make any sense.

Because if you do plot(fft(vector)) then you aren't using any sort of reference for the time axis. You'll see a spike somewhere, but your x-axis is simply incrementing up to the length of your vector. Therefore, you need to do plot(frequency_vector,fft(vector)), although I believe in actuality you want to use abs(fft(vector)) and then take only half of that vector, as I believe the second half is simply repetitive. I could be wrong though, so you're probably better off using the pwelch function

I can get good results if I plot, say, sin(20*pi*t). However, if I plot something with a much lower frequency (like sin(pi*t) or sin(.5*pi*t)), it gets cut off in the lower end when I try to do the PSD plot - the spike is cut off on the left.

If you have sin(.5*pi*t), then the frequency of that signal is .25 Hz. I used pwelch on that and got a spike at the beginning at approximately .25 Hz. Sure the spike is somewhat cut off on the left side, but that's because that happens at 0 Hz. You can't go below that

 

[Stojakapimp]

Ok, by using pwelch on your signal, It looks like a lowpass filter that flattens out at about 5 Hz. Is this what you're getting?

 

[KillerCharlie]

Originally posted by: Stojakapimp

If you have sin(.5*pi*t), then the frequency of that signal is .25 Hz. I used pwelch on that and got a spike at the beginning at approximately .25 Hz. Sure the spike is somewhat cut off on the left side, but that's because that happens at 0 Hz. You can't go below that

The problem is that the frequencies I'll be looking at in one of the equations is 1/3 Hz, 1/6 Hz, and 1/12 Hz so I'm having trouble getting the desired resolution to see them all properly. It looks like playing around with the interval and step size helps though.

 

[Stojakapimp]

Originally posted by: KillerCharlie

Originally posted by: Stojakapimp

If you have sin(.5*pi*t), then the frequency of that signal is .25 Hz. I used pwelch on that and got a spike at the beginning at approximately .25 Hz. Sure the spike is somewhat cut off on the left side, but that's because that happens at 0 Hz. You can't go below that

The problem is that the frequencies I'll be looking at in one of the equations is 1/3 Hz, 1/6 Hz, and 1/12 Hz so I'm having trouble getting the desired resolution to see them all properly. It looks like playing around with the interval and step size helps though.

Most definately....when you use the pwelch command, the graph that comes up will go up to a frequency that is half of the sampling frequency. This is called the Nyquist frequency. So in my examply, since my sampling frequency was 20, the pwelch graph will go up to 10 Hz. Since you are looking at frequencies that are so low, try using a sampling frequency of 2, so that the Pwelch graph goes to 1 Hz. Might be worth a try.

 

[KillerCharlie]

I got the artificial broadband wave consisting of those really low frequencies to come out beautifully, now I just need to plot it against the correct frequency vector.

 

[Stojakapimp]

Originally posted by: KillerCharlie
I got the artificial broadband wave consisting of those really low frequencies to come out beautifully, now I just need to plot it against the correct frequency vector.

If you know what your sampling frequency is, try doing something like this:

freq_vector=(samplingFreq/2)*(0:length(vector)-1/length(vector));