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

# MATLAB code for Power Spectral Density

Started by ●December 10, 2005

Posted by ●December 10, 2005

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

Posted by ●December 10, 2005

"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

Posted by ●December 10, 2005

> >ngeva0 wrote: >> Can anyone please check if my code is correct? SINE512 is just a filename.>> 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!

Posted by ●December 10, 2005

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

Posted by ●December 10, 2005

> >"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/2watt>and your PSD display may either register a 1 volt (peak) input sinewaveas 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 sinewaveinto>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.

Posted by ●December 10, 2005

"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/2watt> and your PSD display may either register a 1 volt (peak) input sinewave as0> 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

Posted by ●December 11, 2005

"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/bandwidthin 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

Posted by ●December 11, 2005

> >"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 gettingclose>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 inthat>fft output bin is spread evenly over the bin , all the >fft.*conj(fft)*appropriate scaling dependent on the number of samples inthe>fft does is give you power in each bin ; it doesn't say anything aboutthe>distribution of power within a bin, for your 50000 samples per secondreal>input , 512 point fft each bin is 50000/512 Hz wide ~100 Hz so the dB/HzPSD>value is about 20 dB less than the dB/Bin value _assuming_ that the poweris>spread evenly through every Hz in that bin - effectively you areconverting>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 toget>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 yourFFT>so you can expect to see some smearing in your FFT output which is >presumably reduced by your hanning window. I am surprised that thespectrum>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 thescope?> >>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 ispower/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 Hzbut>I dont know where you get this from ; at something under 100 Hz per binyou>should get about 0.005. > >> I use my MATLAB code, I also got a spike at 12KHz, but the magnitudeis>> 21. If done correctly, these 2 PSD values should match. But theydon't.> >> To test further, I used different number of samples. Somehow I couldonly>> have a max of 1000 samples (maybe due to the buffer size of theequipment>> 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 itdoesn'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 abit>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.

Posted by ●December 11, 2005

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