How to calculate PSD from FFT?

[gaurav Nanda]

I have calculated FFT of a signal but I have no idea how to calculate PSD from FFT. I have been following the posts and I have come across two ways to find PSD from FFT. The outcome is different in both the cases. Please let me know which one is the right way to do it.

(1) Y=FFT
PSD= (FFT).^2

(2) Y=FFT

PSD= (T/length(y))*FFT

where :
T=sampling time
T=1/Fs,
Fs= sampling frequency
length(y)=signal length


Thanks,
Gaurav

[Nasser M. Abbasi]

On 4/20/2012 6:01 AM, gaurav Nanda wrote:
> I have calculated FFT of a signal but I have no idea how to calculate PSD from FFT.
>I have been following the posts and I have come across two ways to find PSD from FFT.
>The outcome is different in both the cases. Please let me know which one is the
>right way to do it.
>
> (1) Y=FFT
> PSD= (FFT).^2
>

fft is complex?

so, take the absolute value of fft, then square that. This will
give you the PSD.


> (2) Y=FFT
>
> PSD= (T/length(y))*FFT
>

I do not know what that is.

> where :
> T=sampling time
> T=1/Fs,
> Fs= sampling frequency
> length(y)=signal length
>
>
> Thanks,
> Gaurav

--Nasser

[gaurav Nanda]

This is how i am doing it:

x=[0,1,2,3,4,5,6,7,8,9]; % time
y=[1,1.1,2.2,1.3,1.4,1.3,2.9,1.4,1].*10^6; % curremt
plot(x,y);
Fs=50; %sampling frequency
T=1/Fs; % sampling time
L=length(y); %length of the signal
L1=x/T;
t=(0:L-1)*T; % time vector

NFFT=2^nextpow2(L);

y1=fft(y-mean(y))/L;

y2=2*abs(y1(1:NFFT/2+1)); % absolute value of fft

Nc= (y2).^2; % power spectral density

fc=Fs/2*linspace(0,1,NFFT/2+1); %frequency

plot(fc,Nc);

--------------------
is it the correct way to find the PSD?

Also, I need to do the curve fitting and find the coefficients. Could you please tell me how can I do that?

Nc= A+B/f^p <<<<<< how can i find this out?

Thanks,
Gaurav


"Nasser M. Abbasi" <n...@12000.org> wrote in message <jmrgo5$l77$1...@speranza.aioe.org>...

[dbd]

On Friday, April 20, 2012 4:01:12 AM UTC-7, gaurav Nanda wrote:
> I have calculated FFT of a signal but I have no idea how to calculate PSD from FFT.

As I have posted here before:

-begin quote-
The PSD, ESD and power spectrum can be calculated via fft based
methods in Matlab.
Manufacturers of dynamic signal analyzers have provided these
functions, properly scaled, for years. Some have been nice enough to
accurately document their functions and make and keep the
documentation available.

Take a look at "Choose your Units!" from B&K:

http://www.bksv.com/doc/bo0438.pdf

and for more detail, "Signals and Units" on page 29 of:

http://www.bksv.com/doc/bv0031.pdf

For an discussion of the signal processing,
consider pages 5-21 of:

http://www.rssd.esa.int/SP/LISAPATHFINDER/docs/Data_Analysis/GH_FFT.pdf

Particularly "3 Introduction" on page 5 and "9 Scaling the results" on
page 15.
I consider "13 Testing the Algorithm" on page 21 and on as a more
practically oriented discussion of pwelch than the Matlab docs.
-end quote-

Dale B. Dalrymple

[Greg Heath]

N = length(x)
absX = abs(fft(x));
Pxx = absX.^2/N; % Mean-Squared Power Spectrum
PSD = Pxx/Fs; % Power Spectral Density


% Parceval (df = Fs/N):
% Pav = average power


Pav = mean(x.^2) = df*sum(PSD)
Pav = mean(x.^2) = mean(Pxx)
Pav = mean(x.^2) = var(x)+mean(x)^2

Hope this helps.

Greg

[gaurav Nanda]

Thanks a lot, I also need to do the curve fitting. could you please tell me how to FIT THE SPECTRUM IN LOGARITHMIC SCALE by a 1/f^(alpha) function?


Greg Heath <g.h...@verizon.net> wrote in message <b47cb16e-9b87-46ef...@d4g2000vbn.googlegroups.com>...

[dbd]

On Monday, April 23, 2012 4:33:05 AM UTC-7, gaurav Nanda wrote:
> Thanks a lot, I also need to do the curve fitting. could you please tell me how to FIT THE SPECTRUM IN LOGARITHMIC SCALE by a 1/f^(alpha) function?

If you wish to resolve a PSD that falls off, that is, alpha < 0, you must use a window function that has sidelobes falling off faster. The default rectangular window has sidelobes that fall off at 1/f^-1 and so cannot resolve a PSD that falls off at faster than 1/f^-1.

A common choice is the "maximum sidelobe rolloff" window family. See:

http://www.compdsp.com/presentations/Dalrymple/dbd.pdf

pages 6,8,9.

Dale B. Dalrymple

[dbd]

Well, except for my sign error on alpha.

How about:
If you wish to resolve a PSD that falls off, that is, alpha > 0, you must use a window function that has sidelobes falling off faster. The default rectangular window has sidelobes that fall off at 1/f^1 and so cannot resolve a PSD that falls off at faster than 1/f^1.

Dale B. Dalrymple