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');
> 

>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.

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.
�����������������������������������������������������������������������

>
>"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.

"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

"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

>
>"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.

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.
�����������������������������������������������������������������������

>
>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!

"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