PSD Calculation

HI All,
I am an engg student doing a project on FFT. Now I am in the process
of learning DSP.I am stuck with  a problem related to the calculation
of PSD. Could anybody please help me out ??
I have got a 16 point FFT Processor.Now after FFT Processiing I do
have 16 complex o/ps which are in the frequency domain.I have to
calculate the PSD over these points. How do i do that ? I should not
use MATLAB for this. I have to write the code in VHDL.
I downloaded some material which says that After FFT processing I have
to calculate the Power which is the square of the samples, then divide
it by the sampling freuency.
Some other documents says that If I take the squre of the o/p samples,
I will get the PSD. So I am confused over this.

Thanks in advance.

Saumya

saumyajit_tech@yahoo.co.in (Saumyajit) wrote in message news:<9d5c33e7.0312080518.6d08d5e9@posting.google.com>...
> HI All,
> I am an engg student doing a project on FFT. Now I am in the process
> of learning DSP.I am stuck with  a problem related to the calculation
> of PSD. Could anybody please help me out ??
> I have got a 16 point FFT Processor.Now after FFT Processiing I do
> have 16 complex o/ps which are in the frequency domain.I have to
> calculate the PSD over these points. How do i do that ? I should not
> use MATLAB for this. I have to write the code in VHDL.
> I downloaded some material which says that After FFT processing I have
> to calculate the Power which is the square of the samples, then divide
> it by the sampling freuency.
> Some other documents says that If I take the squre of the o/p samples,
> I will get the PSD. So I am confused over this.

First, you have to find the documentation of your FFT routine. Formally, 
the Discrete Fourier Transform (DFT) is defined as something like

              1    N-1             2pi nk
    X(k) = ------- sum x(n) exp(-j-------- )
           sqrt(N) n=0                N

and the Inverse DFT as 

              1    N-1             2pi nk
    x(n) = ------- sum X(k) exp( j-------- ).
           sqrt(N) k=0                N

Note the 1/sqrt(N) factor that pops out in both formulas.

The FFT/IFFT routines usually "simplify" the transforms by dividing 
by N in either the forward or inverse transform. Matlab does, for 
instance, not scale the sum at all in the forward transform but 
divides by N in the inverse transform.

You have to make sure that you skale the DFT with 1/sqrt(N). 

As for scaling with the frequency, we usually want the PSD estimate 
in terms of physichal, not normalized, frequency. So when you have 
an N-point DFT, it spans a physical bandwidth of Fs Hz (Fs is sampling 
frequency) and its frequency bins are dF= Fs/N wide. 

You want a power spectral *density* estimate that relates to the 
physical world, so you want to know how much power each frequency 
bin represents. Remember, power spectral density relates to how 
much power is found inside a given bandwidt, so you need to multiply 
the coefficient in each bin with the bandwidth of the bin, dF.

These scaling factors are often skipped (or treated very sloppy) 
since people often are content with finding an estimate for the 
shape of the spectrum, and do not necessarily care too much about 
the numbers. 

Rune

Saumyajit wrote:

> HI All,
> I am an engg student doing a project on FFT. Now I am in the process
> of learning DSP.I am stuck with  a problem related to the calculation
> of PSD. Could anybody please help me out ??
> I have got a 16 point FFT Processor.Now after FFT Processiing I do
> have 16 complex o/ps which are in the frequency domain.I have to
> calculate the PSD over these points. How do i do that ? I should not
> use MATLAB for this. I have to write the code in VHDL.
> I downloaded some material which says that After FFT processing I have
> to calculate the Power which is the square of the samples, then divide
> it by the sampling freuency.
> Some other documents says that If I take the squre of the o/p samples,
> I will get the PSD. So I am confused over this.
>
> Thanks in advance.
>
> Saumya

You need the nearest approximation which is called the periodogram - found
as follows.

Suppose X(k) are the FFT frequency samples found in the usual way (with no
division by N or sqrt(N) by the way).
Then the periodogram is mag(X(k))^2/N
where mag is the sqrt(real^2+imag^2) ie the usual modulus of a complex
number. Here we square it and divide by N, the number of FFT points. There
is another way you can do it if you know X(k). You can do this

P(i)=P(i-1)*beta+(1-beta)*X(k)*X'(k)

where X'(k) is the complex conjugate of X(k) at FFT frame i and beta is a
forgetting factor 0<beta<1 the choice of which needs a little
experimenting with. P(i) is the periodogram at FFT frame i. The forgetting
fact beta should be close to unity for good smoothing and perhaps nearer
to 0.6? if you need good tracking ie the data is highly non-stationary
like speech. If it is staionary data then a beta close to unity will make
the periodogram 'converge' to the PSD.

Tom