Hi,
If you're looking to learn how to compute accurate spectral estimates, a
particular favorite text around here is Stoica & Moses (Intro to Spectral
Analysis).
If you want to use preprogrammed software for this, I encourage you to use
two methods in the Signal Processing Toolbox (if you have it). Type "help
spectrum/psd" and "help spectrum/msspectrum" to learn how to compute PSD and
mean-square spectrum estimates. These are two well-defined estimates that
are useful in practice. Quick examples taken from the methods cited above
to get you started:
PSD EXAMPLE (if you're looking for a spectral density on, say, a noisy comms
signal or a general continuous spectral distribution):
% Spectral analysis of a signal that contains a 200Hz cosine plus
noise.
Fs = 1000; t = 0:1/Fs:.296;
x = cos(2*pi*t*200)+randn(size(t));
h = spectrum.welch; % Instantiate a welch object.
psd(h,x,'Fs',Fs); % Plot the one-sided PSD.
MSS EXAMPLE (if you have distinct sines, e.g. discrete spectra, and want
amplitude levels - sounds like your example):
% In this example we will measure the power of a deterministic
% power signal which has a frequency component at 200Hz. We'll use a
% signal with a peak amplitude of 3 volts therefore, the theoretical
% power at 200Hz should be 3^2/2 volts^2 (watts) or 6.53dB.
Fs = 1000; t = 0:1/Fs:.2;
x = 3*cos(2*pi*t*200);
h = spectrum.welch; % Instantiate a welch object.
% Plot the one-sided Mean-square spectrum.
h.FFTLength = 'UserDefined';
msspectrum(h,x,'Fs',Fs,'NFFT',2^14);
Regards,
--Don
"ngeva0" <ngeva0@hotmail.com> wrote in message
news:ne-dnR4SaOiNswfenZ2dnUVZ_tqdnZ2d@giganews.com...
> Can anyone please check if my code is correct? SINE512 is just a file
> name.
> This file has 512 data points and it's a sine wave. Thank you!!
>
> %Sampling frequency
> Fs=50000;
>
> %# of samples in the data
> datasize=size(SINE512);
> numsample=datasize(1);
>
> %Windowing
> H=hann(numsample);
> W=H.*(SINE512(:,2));
>
> %Fourier Transform
> FFTX=fft(W,numsample);
>
> %Power: magnitude^2
> X=FFTX(1:floor(numsample/2)).*conj(FFTX(1:floor(numsample/2)));
>
> %Bandwidth
> BW=1.5*Fs*numsample;
>
> %PSD=magnitude^2/bandwidth
> PSD=X/BW;
>
> %Computing the corresponding frequency values
> Omega=Fs*(0:numsample-1)/numsample;
> Omega=Omega(1:floor(numsample/2));
>
> %Plot PSD
> H=plot(Omega,PSD);
> set(H,'Color','BLACK');
>
Reply by ngeva0●December 11, 20052005-12-11
>ngeva0 wrote:
>
>> ... Can you explain what does the "2 volts
>> peak to peak" mean? I don't get that. Thank you.
>
>Consider the function y = A*sin(wt). the positive peak value is +A and
>is easy to see on an oscilloscope. The value for y = zero is hard to
>estimate accurately, but the negative peak, -A, is easy to see. We
>measure the waveform as having a peak-to-peak value of 2*A and often
>describe it that way. For any function, peak to peak is the entire range
>that it needs.
>
>Jerry
>--
>Engineering is the art of making what you want from things you can get.
>�����������������������������������������������������������������������
>
Thank you.
Reply by Jerry Avins●December 11, 20052005-12-11
ngeva0 wrote:
> ... Can you explain what does the "2 volts
> peak to peak" mean? I don't get that. Thank you.
Consider the function y = A*sin(wt). the positive peak value is +A and
is easy to see on an oscilloscope. The value for y = zero is hard to
estimate accurately, but the negative peak, -A, is easy to see. We
measure the waveform as having a peak-to-peak value of 2*A and often
describe it that way. For any function, peak to peak is the entire range
that it needs.
Jerry
--
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������
Reply by ngeva0●December 11, 20052005-12-11
>
>"ngeva0" <ngeva0@hotmail.com> wrote in message
>news:E5ednYmLdYH7tQbeRVn-gg@giganews.com...
>> >
><snip>
>>>>> %Bandwidth
>>>>> BW=1.5*Fs*numsample;
>
> where did you derive this from?
>Fs/numsample is the width of each of the frequency bins in the FFT output
>where does the 1.5 come from and why are you dividing (later) by
>FS*numsample instead of by fft/numsample - I can only see this getting
close
>if you need to scale your fft output power by dividing by numsample^2 (
>which is what I have to do with scilab ) - I can't see any reason for the
>factor of 2/3.
>
>
>>>>>
>>>>> %PSD=magnitude^2/bandwidth
>in dBX/Hz if bandwidth is in Hz - this assumes that the total power in
that
>fft output bin is spread evenly over the bin , all the
>fft.*conj(fft)*appropriate scaling dependent on the number of samples in
the
>fft does is give you power in each bin ; it doesn't say anything about
the
>distribution of power within a bin, for your 50000 samples per second
real
>input , 512 point fft each bin is 50000/512 Hz wide ~100 Hz so the dB/Hz
PSD
>value is about 20 dB less than the dB/Bin value _assuming_ that the power
is
>spread evenly through every Hz in that bin - effectively you are
converting
>power in a 100 Hz wide range to average power in each 1 Hz of that range.
>Your test sinusoid's PSD is not flat in that range so applying this
>bandwidth correction factor will mess up your displayed values for
>sinusoidal inputs, if you are putting in white noise which is ( averaged
>over many independent FFTs) flat across the bin, then you can expect to
get
>a good estimate of PSD by applying this ( -10*log10(Fs/numsample) )
>bandwidth correction factor.
>
><more snip>
>>
>> I did a test with a sine wave. I used a 12KHz sine wave with 1 volt
>> amplitude.
>On the scope did this appear as 2 Volts peak to peak ?
>
>On the spectrum analyzer, it shows a spike at 12KHz with a peak
>> value of -3dBV(this peak value changes if I use a different sine wave
>> frequency).
>At 12 kHz and 50000 samples /sec 512 samples is ~ 122.9 cycles in your
FFT
>so you can expect to see some smearing in your FFT output which is
>presumably reduced by your hanning window. I am surprised that the
spectrum
>analyser display changes significantly as you change your input frequency
>though ; what differences do you see if you put in 3125 Hz sinewave
>compared with 6250 Hz ? Do they both show the same magnitude on the
scope?
>
>>I asked my professor and he said that I could convert that to
>> PSD by using V^2 = 10^(dBV/10) and V^2/BW = PSD. V^2 is
power/magnitude
>> square. If I use this equation to convert -3dBV to PSD, I got .0021.
>
>!0^(-0.3) ~ .5 so you expect your bandwidth per bin to be around 240 Hz
but
>I dont know where you get this from ; at something under 100 Hz per bin
you
>should get about 0.005.
>
>> I use my MATLAB code, I also got a spike at 12KHz, but the magnitude
is
>> 21. If done correctly, these 2 PSD values should match. But they
don't.
>
>> To test further, I used different number of samples. Somehow I could
only
>> have a max of 1000 samples (maybe due to the buffer size of the
equipment
>> I use). So I sampled 1000 data points. But to do the FFT, you need the
>> number of samples to be the power of 2 I think.
>not in matlab it just does a dft , you can input 23 points if you want,
>provided that you don't declare what size fft you want because it will
>truncate or pad with zeros to the declared length otherwise. What it
doesn't
>seem to do is scale the outputs so that parsevals theorem is met ( i.e.
>total energy in the fft o/p is > total energy in ).
>> If anything you think is confusing, please let me know.
>I've done the best I can to let you know what is confusing me by
>interleaving comments above ( but please cross-check my arithmetic - I'm
>rotten at it.) - mainly confused that your bandwidth scaling factor is a
bit
>odd and seems to incorporate the FFT scaling too. You might be clearer in
>your own mind if you separate the two.
>
>Best of Luck - Mike
>
>
>
>
>
Thanks for your very detail explaination. I need to look into it before
providing any further information. Can you explain what does the "2 volts
peak to peak" mean? I don't get that. Thank you.
Reply by Mike Yarwood●December 10, 20052005-12-10
"ngeva0" <ngeva0@hotmail.com> wrote in message
news:E5ednYmLdYH7tQbeRVn-gg@giganews.com...
> >
<snip>
>>>> %Bandwidth
>>>> BW=1.5*Fs*numsample;
where did you derive this from?
Fs/numsample is the width of each of the frequency bins in the FFT output
where does the 1.5 come from and why are you dividing (later) by
FS*numsample instead of by fft/numsample - I can only see this getting close
if you need to scale your fft output power by dividing by numsample^2 (
which is what I have to do with scilab ) - I can't see any reason for the
factor of 2/3.
>>>>
>>>> %PSD=magnitude^2/bandwidth
in dBX/Hz if bandwidth is in Hz - this assumes that the total power in that
fft output bin is spread evenly over the bin , all the
fft.*conj(fft)*appropriate scaling dependent on the number of samples in the
fft does is give you power in each bin ; it doesn't say anything about the
distribution of power within a bin, for your 50000 samples per second real
input , 512 point fft each bin is 50000/512 Hz wide ~100 Hz so the dB/Hz PSD
value is about 20 dB less than the dB/Bin value _assuming_ that the power is
spread evenly through every Hz in that bin - effectively you are converting
power in a 100 Hz wide range to average power in each 1 Hz of that range.
Your test sinusoid's PSD is not flat in that range so applying this
bandwidth correction factor will mess up your displayed values for
sinusoidal inputs, if you are putting in white noise which is ( averaged
over many independent FFTs) flat across the bin, then you can expect to get
a good estimate of PSD by applying this ( -10*log10(Fs/numsample) )
bandwidth correction factor.
<more snip>
>
> I did a test with a sine wave. I used a 12KHz sine wave with 1 volt
> amplitude.
On the scope did this appear as 2 Volts peak to peak ?
On the spectrum analyzer, it shows a spike at 12KHz with a peak
> value of -3dBV(this peak value changes if I use a different sine wave
> frequency).
At 12 kHz and 50000 samples /sec 512 samples is ~ 122.9 cycles in your FFT
so you can expect to see some smearing in your FFT output which is
presumably reduced by your hanning window. I am surprised that the spectrum
analyser display changes significantly as you change your input frequency
though ; what differences do you see if you put in 3125 Hz sinewave
compared with 6250 Hz ? Do they both show the same magnitude on the scope?
>I asked my professor and he said that I could convert that to
> PSD by using V^2 = 10^(dBV/10) and V^2/BW = PSD. V^2 is power/magnitude
> square. If I use this equation to convert -3dBV to PSD, I got .0021.
!0^(-0.3) ~ .5 so you expect your bandwidth per bin to be around 240 Hz but
I dont know where you get this from ; at something under 100 Hz per bin you
should get about 0.005.
> I use my MATLAB code, I also got a spike at 12KHz, but the magnitude is
> 21. If done correctly, these 2 PSD values should match. But they don't.
> To test further, I used different number of samples. Somehow I could only
> have a max of 1000 samples (maybe due to the buffer size of the equipment
> I use). So I sampled 1000 data points. But to do the FFT, you need the
> number of samples to be the power of 2 I think.
not in matlab it just does a dft , you can input 23 points if you want,
provided that you don't declare what size fft you want because it will
truncate or pad with zeros to the declared length otherwise. What it doesn't
seem to do is scale the outputs so that parsevals theorem is met ( i.e.
total energy in the fft o/p is > total energy in ).
> If anything you think is confusing, please let me know.
I've done the best I can to let you know what is confusing me by
interleaving comments above ( but please cross-check my arithmetic - I'm
rotten at it.) - mainly confused that your bandwidth scaling factor is a bit
odd and seems to incorporate the FFT scaling too. You might be clearer in
your own mind if you separate the two.
Best of Luck - Mike
Reply by Glennaebad●December 10, 20052005-12-10
"Mike Yarwood" <mpyarwood@btopenworld.com> wrote in message
news:dnet0d$d1k$1@nwrdmz03.dmz.ncs.ea.ibs-infra.bt.com...
>
> "john" <johns@xetron.com> wrote in message
> news:1134210735.813800.60750@g47g2000cwa.googlegroups.com...
> >
> > ngeva0 wrote:
> >> Can anyone please check if my code is correct? SINE512 is just a file
> >> name.
> >> This file has 512 data points and it's a sine wave. Thank you!!
> >>
> >> %Sampling frequency
> >> Fs=50000;
> >>
> >> %# of samples in the data
> >> datasize=size(SINE512);
> >> numsample=datasize(1);
> >>
> >> %Windowing
> >> H=hann(numsample);
> >> W=H.*(SINE512(:,2));
> >>
> >> %Fourier Transform
> >> FFTX=fft(W,numsample);
> >>
> >> %Power: magnitude^2
> >> X=FFTX(1:floor(numsample/2)).*conj(FFTX(1:floor(numsample/2)));
> >>
> >> %Bandwidth
> >> BW=1.5*Fs*numsample;
> >>
> >> %PSD=magnitude^2/bandwidth
> >> PSD=X/BW;
> >>
> >> %Computing the corresponding frequency values
> >> Omega=Fs*(0:numsample-1)/numsample;
> >> Omega=Omega(1:floor(numsample/2));
> >>
> >> %Plot PSD
> >> H=plot(Omega,PSD);
> >> set(H,'Color','BLACK');
> >
> > I am happy to report that your code behaves exactly as programmed.
> > Congratulations.
> >
> But you may want to look at :
> http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.html#998360
>
> where you can see that they normalise power in the bin by dividing by the
> number of points in the DFT and get powers of about 2*pi for their unit
> peak amplitude sinewaves.
>
> A 1Volt peak amplitude sinewave across a 1 ohm resistor dissipates 1/2
watt
> and your PSD display may either register a 1 volt (peak) input sinewave as
0
> dBV (w.r.t. 1V pk), -3.01 dBV (w.r.t. 1 V r.m.s) or -6.02 dBV (w.r.t. 1
> Vpk-pk but this is very unlikely). So put a known amplitude sinewave into
> your instrument (with as little extra noise as possible) and see what it
> displays.
>
> None of the above? has it scaled by a further 10*log10(512/50000) ~ -19.9
> dB ?
>
> Now you know what the instrument is doing you can decide what your matlab
> routine should do.
>
>
> Best of Luck - Mike
>
>
Remember its the AREA under the PSD that gives power - not the peak values.
For a sine wave on a bin this means you need to multiply by the width of the
bin.
Glen
Reply by ngeva0●December 10, 20052005-12-10
>
>"john" <johns@xetron.com> wrote in message
>news:1134210735.813800.60750@g47g2000cwa.googlegroups.com...
>>
>> ngeva0 wrote:
>>> Can anyone please check if my code is correct? SINE512 is just a file
>>> name.
>>> This file has 512 data points and it's a sine wave. Thank you!!
>>>
>>> %Sampling frequency
>>> Fs=50000;
>>>
>>> %# of samples in the data
>>> datasize=size(SINE512);
>>> numsample=datasize(1);
>>>
>>> %Windowing
>>> H=hann(numsample);
>>> W=H.*(SINE512(:,2));
>>>
>>> %Fourier Transform
>>> FFTX=fft(W,numsample);
>>>
>>> %Power: magnitude^2
>>> X=FFTX(1:floor(numsample/2)).*conj(FFTX(1:floor(numsample/2)));
>>>
>>> %Bandwidth
>>> BW=1.5*Fs*numsample;
>>>
>>> %PSD=magnitude^2/bandwidth
>>> PSD=X/BW;
>>>
>>> %Computing the corresponding frequency values
>>> Omega=Fs*(0:numsample-1)/numsample;
>>> Omega=Omega(1:floor(numsample/2));
>>>
>>> %Plot PSD
>>> H=plot(Omega,PSD);
>>> set(H,'Color','BLACK');
>>
>> I am happy to report that your code behaves exactly as programmed.
>> Congratulations.
>>
>But you may want to look at :
>http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.html#998360
>
>where you can see that they normalise power in the bin by dividing by the
>number of points in the DFT and get powers of about 2*pi for their unit
>peak amplitude sinewaves.
>
>A 1Volt peak amplitude sinewave across a 1 ohm resistor dissipates 1/2
watt
>and your PSD display may either register a 1 volt (peak) input sinewave
as 0
>dBV (w.r.t. 1V pk), -3.01 dBV (w.r.t. 1 V r.m.s) or -6.02 dBV (w.r.t. 1
>Vpk-pk but this is very unlikely). So put a known amplitude sinewave
into
>your instrument (with as little extra noise as possible) and see what it
>displays.
>
>None of the above? has it scaled by a further 10*log10(512/50000) ~
-19.9
>dB ?
>
>Now you know what the instrument is doing you can decide what your matlab
>routine should do.
>
>
>Best of Luck - Mike
>
>
>
Hi Mike:
I did a test with a sine wave. I used a 12KHz sine wave with 1 volt
amplitude. On the spectrum analyzer, it shows a spike at 12KHz with a peak
value of -3dBV(this peak value changes if I use a different sine wave
frequency). I asked my professor and he said that I could convert that to
PSD by using V^2 = 10^(dBV/10) and V^2/BW = PSD. V^2 is power/magnitude
square. If I use this equation to convert -3dBV to PSD, I got .0021. When
I use my MATLAB code, I also got a spike at 12KHz, but the magnitude is
.21. If done correctly, these 2 PSD values should match. But they don't.
To test further, I used different number of samples. Somehow I could only
have a max of 1000 samples (maybe due to the buffer size of the equipment
I use). So I sampled 1000 data points. But to do the FFT, you need the
number of samples to be the power of 2 I think. So I only used the first
512 points and deleted the rest data. I am not sure if that makes any
difference if I sampled 512 data points directly instead of doing this
way. The PSD peak value changes when I used different number of samples,
512 and 1000, when I used my MATLAB code to plot it. I think this makes
sense because the frequency range decreases with increases number of
samples, so the PSD peak value has to increases, and vice versa ( because
the range of frequeny times the PSD magnitude = power?)
If anything you think is confusing, please let me know. Thanks for you
advice.
Reply by Jerry Avins●December 10, 20052005-12-10
ngeva0 wrote:
...
>>I am happy to report that your code behaves exactly as programmed.
>>Congratulations.
>>
>>John
>>
>>
>
>
> I was not sure if the BW should be 1.5*Fs*numsample or just 1.5*Fs. So I
> guess what I have is correct. Thanks!
Don't get too happy. John says that it does what it's programmed to do,
not that what it does is right.
Jerry
--
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������
Reply by ngeva0●December 10, 20052005-12-10
>
>ngeva0 wrote:
>> Can anyone please check if my code is correct? SINE512 is just a file
name.
>> This file has 512 data points and it's a sine wave. Thank you!!
>>
>> %Sampling frequency
>> Fs=50000;
>>
>> %# of samples in the data
>> datasize=size(SINE512);
>> numsample=datasize(1);
>>
>> %Windowing
>> H=hann(numsample);
>> W=H.*(SINE512(:,2));
>>
>> %Fourier Transform
>> FFTX=fft(W,numsample);
>>
>> %Power: magnitude^2
>> X=FFTX(1:floor(numsample/2)).*conj(FFTX(1:floor(numsample/2)));
>>
>> %Bandwidth
>> BW=1.5*Fs*numsample;
>>
>> %PSD=magnitude^2/bandwidth
>> PSD=X/BW;
>>
>> %Computing the corresponding frequency values
>> Omega=Fs*(0:numsample-1)/numsample;
>> Omega=Omega(1:floor(numsample/2));
>>
>> %Plot PSD
>> H=plot(Omega,PSD);
>> set(H,'Color','BLACK');
>
>I am happy to report that your code behaves exactly as programmed.
>Congratulations.
>
>John
>
>
I was not sure if the BW should be 1.5*Fs*numsample or just 1.5*Fs. So I
guess what I have is correct. Thanks!
Reply by Mike Yarwood●December 10, 20052005-12-10
"john" <johns@xetron.com> wrote in message
news:1134210735.813800.60750@g47g2000cwa.googlegroups.com...
>
> ngeva0 wrote:
>> Can anyone please check if my code is correct? SINE512 is just a file
>> name.
>> This file has 512 data points and it's a sine wave. Thank you!!
>>
>> %Sampling frequency
>> Fs=50000;
>>
>> %# of samples in the data
>> datasize=size(SINE512);
>> numsample=datasize(1);
>>
>> %Windowing
>> H=hann(numsample);
>> W=H.*(SINE512(:,2));
>>
>> %Fourier Transform
>> FFTX=fft(W,numsample);
>>
>> %Power: magnitude^2
>> X=FFTX(1:floor(numsample/2)).*conj(FFTX(1:floor(numsample/2)));
>>
>> %Bandwidth
>> BW=1.5*Fs*numsample;
>>
>> %PSD=magnitude^2/bandwidth
>> PSD=X/BW;
>>
>> %Computing the corresponding frequency values
>> Omega=Fs*(0:numsample-1)/numsample;
>> Omega=Omega(1:floor(numsample/2));
>>
>> %Plot PSD
>> H=plot(Omega,PSD);
>> set(H,'Color','BLACK');
>
> I am happy to report that your code behaves exactly as programmed.
> Congratulations.
>
But you may want to look at :
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.html#998360
where you can see that they normalise power in the bin by dividing by the
number of points in the DFT and get powers of about 2*pi for their unit
peak amplitude sinewaves.
A 1Volt peak amplitude sinewave across a 1 ohm resistor dissipates 1/2 watt
and your PSD display may either register a 1 volt (peak) input sinewave as 0
dBV (w.r.t. 1V pk), -3.01 dBV (w.r.t. 1 V r.m.s) or -6.02 dBV (w.r.t. 1
Vpk-pk but this is very unlikely). So put a known amplitude sinewave into
your instrument (with as little extra noise as possible) and see what it
displays.
None of the above? has it scaled by a further 10*log10(512/50000) ~ -19.9
dB ?
Now you know what the instrument is doing you can decide what your matlab
routine should do.
Best of Luck - Mike