Power Spectral Density of a Signal in decibels

Hello,

I am calculating the Power Spectral Density of a signal using fft as recommended in the matlab demo section. I would like to display the fft analysis in the same manner as you would see on a third octave band analyser sound level meter.

I have the following code:

M1Afft = fft(M1A,length(M1A));

PowM1Afft = M1Afft.*conj(M1Afft)/length(M1A);

plot(f,PowM1Afft(1:length(M1A)))

axis([0,2000,0,1.5])

title('Measurement 1')

xlabel('Frequency (Hz)')

I am unsure of the Y axis units, but i would like them in decibels based on sound power of the frequency. I am also having a quite a bit of problems getting a histogram to produce frequency bins in the manner of a third octave band analyser.

Any help would be greatly appreciated.

Many Thanks, Andy

Answer 1

[Wayne King]

Do you have the Signal Processing Toolbox?

Fs = 1e3;
t = 0:.001:1-0.001;
x = cos(2*pi*100*t)+randn(size(t));
plot(psd(spectrum.periodogram,x,'Fs',1e3,'NFFT',length(x)));

[Andy]

Thanks for your help Wayne,

I do have the signal processing toolbox, however, i think i may have explained the question wrong, or may have been using the wrong code to start with.

I am playing some white noise through a speaker array, and have four microphones at various points around the room.

I am looking to display the measured data as Lp, or sound pressure in decibels on the Y axis, and frequency on the X axis in the manner seen on most third octave band analysers.

I think i may have been wrong in using the spectral density.

Regards,

Andy

[Wayne King]

Hi Andy, in order to do that you will need more information like the microphone sensitvity. See

http://www.mathworks.com/matlabcentral/answers/10962-how-to-convert-pascal-to-db-knowing-the-microphone-sensitivity

[Andy]

Thanks

Answer 2

[Wayne King]

If you want to do the above with just fft()

    xdft = fft(x);
    % even length x
    psdest = abs(xdft(1:length(xdft)/2+1)).^2;
    psdest = (1/(Fs*length(x))).*psdest;
    psdest(2:end-1) = 2*psdest(2:end-1);
    xdft = (1/(Fs*length(x))).*xdft;
    f = 0:Fs/length(x):Fs/2;
    plot(f,10*log10(psdest)); grid on;

[Andy]

I am confused as to why there are negative values in the x axis. If a signal is generated, then all values would be positive. Or am i missing something?

[Andy]

Y axis,sorry.

[Wayne King]

It just means you have less than 1 (but greater than 0)