Sign in

Not a member? | Forgot your Password?

Search compdsp

Search tips

Free PDF Downloads

A Quadrature Signals Tutorial: Complex, But Not Complicated

Understanding the 'Phasing Method' of Single Sideband Demodulation

Complex Digital Signal Processing in Telecommunications

Introduction to Sound Processing

C++ Tutorial

Introduction of C Programming for DSP Applications

Fixed-Point Arithmetic: An Introduction

Cascaded Integrator-Comb (CIC) Filter Introduction

Discussion Groups

FFT Spectral Analysis Software

Free Online Books

See Also

Embedded SystemsFPGA

Discussion Groups | Comp.DSP | MATLAB code for Power Spectral Density

There are 12 messages in this thread.

You are currently looking at messages 1 to .


Is this discussion worth a thumbs up?

0

MATLAB code for Power Spectral Density - ngeva0 - 2005-12-09 20:23:00



Re: MATLAB code for Power Spectral Density - john - 2005-12-10 05:32:00



Re: MATLAB code for Power Spectral Density - Mike Yarwood - 2005-12-10 10:42:00



Re: MATLAB code for Power Spectral Density - ngeva0 - 2005-12-10 13:35:00



Re: MATLAB code for Power Spectral Density - Jerry Avins - 2005-12-10 13:48:00



Re: MATLAB code for Power Spectral Density - ngeva0 - 2005-12-10 14:11:00



Re: MATLAB code for Power Spectral Density - Glennaebad - 2005-12-10 17:45:00



Re: MATLAB code for Power Spectral Density - Mike Yarwood - 2005-12-10 20:27:00

"ngeva0" <n...@hotmail.com> wrote in message 
news:E...@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


Re: MATLAB code for Power Spectral Density - ngeva0 - 2005-12-10 23:49:00

>
>"ngeva0" <n...@hotmail.com> wrote in message 
>news:E...@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.


Re: MATLAB code for Power Spectral Density - Jerry Avins - 2005-12-11 00:02:00



| 1 | |