power spectral density versus magnitude spectrum

A friend and I have been analyzing a signal x(t) containing both random and
known-period deterministic variations.  So, I wrote a simple Matlab DFT
algorithm for calculating the periodogram of x(t).  My friend asked me to
scale the calculated power and plot my results as an "Amplitude Spectrum",
in order that one could read off the amplitudes of the underlying
sinusoids.  Following Parseval's Theorem I decided to scale the squared
DFT output in order that the average signal power is the integral of my
PSD, (sum(x.^2)/length(x)) = sum(Pxx), where Pxx are the densities.  

At first, it seemed to me that the amplitudes at a given fourier frequency
f could then then be recovered as sqrt(2*Pxx(f)).  Of course, this is bogus
since these are "densities" and energy for f will also be contained in
nearby frequency bins.   For simple simulations where I choose rectangular
windows and generate periodics exactly at fourier frequencies, no problem. 
If real-life were only so simple.  For other windows, energy gets smeared.

So, my questions, I guess, are:

1.) Am I missing something simple here and can the amplitudes actually be
recovered by some simple transformation on the Pxx?

2.) Matlab has a msspectrum method in its signal processing toolbox which
claims to generate something similar to what my friend requested, but
there's not additional information on how it is calculated.  It seems they
have first constructed the autocorrelation for x and from there ...?

Thanks for the help.

Ken

On Aug 13, 5:20 am, "seniork" <seni...@gmail.com> wrote:
> A friend and I have been analyzing a signal x(t) containing both random and
> known-period deterministic variations.  So, I wrote a simple Matlab DFT
> algorithm for calculating the periodogram of x(t).  My friend asked me to
> scale the calculated power and plot my results as an "Amplitude Spectrum",
> in order that one could read off the amplitudes of the underlying
> sinusoids.  Following Parseval's Theorem I decided to scale the squared
> DFT output in order that the average signal power is the integral of my
> PSD, (sum(x.^2)/length(x)) = sum(Pxx), where Pxx are the densities.
>
> At first, it seemed to me that the amplitudes at a given fourier frequency
> f could then then be recovered as sqrt(2*Pxx(f)).  Of course, this is bogus
> since these are "densities" and energy for f will also be contained in
> nearby frequency bins.   For simple simulations where I choose rectangular
> windows and generate periodics exactly at fourier frequencies, no problem.
> If real-life were only so simple.  For other windows, energy gets smeared.
>
> So, my questions, I guess, are:
>
> 1.) Am I missing something simple here and can the amplitudes actually be
> recovered by some simple transformation on the Pxx?
>
> 2.) Matlab has a msspectrum method in its signal processing toolbox which
> claims to generate something similar to what my friend requested, but
> there's not additional information on how it is calculated.  It seems they
> have first constructed the autocorrelation for x and from there ...?
>
> Thanks for the help.
>
> Ken

Ken

Apply your 'at first' approach with Matlab's flattop window. The
flattop windows effectively perform the integration over the 'smeared
region' and provide amplitude accuracy of < 0.01 dB for pure tones,
regardless of bin alignment. There are also other flattop windows with
higher sidelobe rejection than Matlab's.

Whether your Pxx are densities or not depends on what the  transformed
signal is. If the transformed signal consists of samples of a
continuous distribution, the Pxx are densities. If the transformed
signal consists of samples of the sum of bin centered tones, the Pxx
are power spectra. Real signals can contain some of each.

Dale B. Dalrymple
http://dbdimages.com

On Aug 13, 5:20 am, "seniork" <seni...@gmail.com> wrote:
> A friend and I have been analyzing a signal x(t) containing both random and
> known-period deterministic variations.  So, I wrote a simple Matlab DFT
> algorithm for calculating the periodogram of x(t).  My friend asked me to
> scale the calculated power and plot my results as an "Amplitude Spectrum",
> in order that one could read off the amplitudes of the underlying
> sinusoids.  Following Parseval's Theorem I decided to scale the squared
> DFT output in order that the average signal power is the integral of my
> PSD, (sum(x.^2)/length(x)) = sum(Pxx), where Pxx are the densities.
>
> At first, it seemed to me that the amplitudes at a given fourier frequency
> f could then then be recovered as sqrt(2*Pxx(f)).  Of course, this is bogus
> since these are "densities" and energy for f will also be contained in
> nearby frequency bins.   For simple simulations where I choose rectangular
> windows and generate periodics exactly at fourier frequencies, no problem.
> If real-life were only so simple.  For other windows, energy gets smeared.
>
> So, my questions, I guess, are:
>
> 1.) Am I missing something simple here and can the amplitudes actually be
> recovered by some simple transformation on the Pxx?

If you have known frequency sinusoids which are not periodic
within your DFT length, you can calculate a "scalloping
loss" from the bin offset mapped onto the FT of the window
(the relative distance down the side of a sinc lobe for
a rectangular window).  You can then use the reciprocal of
this loss ratio to estimate amplitudes; and double conversely,
use this loss ratio curve to subtract the "leakage" of your
sinusoid from all the adjacent bins, something the use of a
"flat top" window may not do for you, if you need a better
estimate of nearby spectral noise as well.


IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M

Thanks for the reply, Dale.

I tried this and I am a bit more confused, at least on your first point. 
For example, I generated a pure sinusoid with amplitude 7 sampled so that
the frequency f0 of this sinusoid was one of my Fourier frequencies. 
Next, I performed three different windows to obtain:

rectwin: sqrt(2*pxx(f0)) = 7
flattopwin: sqrt(2*pxx(f0)) = 3.606
blackmanharris: sqrt(2*pxx(f0)) = 4.944

Clearly, flattopwin is performing the worst in recovering the amplitude of
7, having the most "smearing" of the three.

I understand your second point.  In fact, we do have random (continuous
spectra) as well as deterministic (discrete) signals present.  But, my
interest is only in exactly quantifying the amplitudes of the
deterministic components.

Ron N. wrote:

"If you have known frequency sinusoids which are not periodic
within your DFT length, you can calculate a "scalloping
loss" from the bin offset mapped onto the FT of the window
(the relative distance down the side of a sinc lobe for
a rectangular window)."

If you know the frequencies, you don't need a DFT.

 "You can then use the reciprocal of
this loss ratio to estimate amplitudes; and double conversely,
use this loss ratio curve to subtract the "leakage" of your
sinusoid from all the adjacent bins, something the use of a
"flat top" window may not do for you, if you need a better
estimate of nearby spectral noise as well.

Pages 28-31 of:
http://www.bksv.com/pdf/Bv0032.pdf
give a similar method of amplitude estimation applied to the Hann
window without knowledge of the frequency.

To perform an accurate subtraction of a tone it would be better to
estimate the tone phase, amplitude and frequency from the complex DFT
coefficients and subtract coherently. For the rectangular window there
is a closed form solution in terms of the complex coefficients of the
bin with peak magnitude and the adjacent bins on each side.

None of these approaches work well without adequate SNR.

Dale B. Dalrymple
http://dbdimages.com

On Aug 13, 12:14 pm, dbd <d...@ieee.org> wrote:
> Ron N. wrote:
>
> "If you have known frequency sinusoids which are not periodic
> within your DFT length, you can calculate a "scalloping
> loss" from the bin offset mapped onto the FT of the window
> (the relative distance down the side of a sinc lobe for
> a rectangular window)."
>
> If you know the frequencies, you don't need a DFT.

If you know all the frequencies.  The OP seems to
imply some mix of known and unknown spectra for which
an estimate of the amplitudes is required.

>  "You can then use the reciprocal of
> this loss ratio to estimate amplitudes; and double conversely,
> use this loss ratio curve to subtract the "leakage" of your
> sinusoid from all the adjacent bins, something the use of a
> "flat top" window may not do for you, if you need a better
> estimate of nearby spectral noise as well.
>
> Pages 28-31 of:http://www.bksv.com/pdf/Bv0032.pdf
> give a similar method of amplitude estimation applied to the Hann
> window without knowledge of the frequency.
>
> To perform an accurate subtraction of a tone it would be better to
> estimate the tone phase, amplitude and frequency from the complex DFT
> coefficients and subtract coherently.

And one can get an even more accurate estimates by using
multiple overlapped fft frames (e.g. one can also estimate
slight amplitude ramps), in the style of a phase vocoder,
before coherent spectral subtraction.  I've been trying
variations on this technique to try and better separate out
and characterize the evolution of the slightly inharmonic
overtones of low octave piano notes. (That's what I get
for trying to tune a cheap spinet by myself and understand
why it still doesn't sound good... :)


IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M
 http://www.nicholson.com/rhn/dsp.html

Thanks guys for all the help.  It seems the more I read and the more I
think I understand, the more confused I become.

I read the article on compensating for the "picket fence" effect using
Hanning and rectangular windows.  But, the numbers just aren't working out
for me in simulation.  E.g., I generated a simple sinusoid having amplitude
7 at a given frequency f0.  I then calculated the PSD with Hanning
windowing, ensuring that the known frequency f0 was one of my Fourier
frequencies.  According to the article, this means that no "amplitude"
correction is necessary and that that the two adjacent bins should be 6 dB
below the peak.  I do not obtain a 6 dB difference between the peak and
adjacent bins (as shown below) and I also don't follow how the article
recommends calculating the amplitude.

pxx(left of f0) = 4.083
pxx(f0) = 16.338
pxx(right of f0) = 4.083

Notice that 16 - 4 = 12 is not 6 dB.  What am I missing? Argh!

Now, since all of the power in my signal is actually 16.338 + 4.083 +
4.083 ~ 24.5 and since 7 = sqrt(2*(24.5)), wouldn't the calculation of the
amplitude have to involve *integrating* the PSD?  Unfortunately, the
article is not clear on this point.  It seems to me the best approach is
to sum up all of the PSD values believed to contribute to a given
frequency (i.e., surrounding bins), multiply by 2 and take the square
root.  I see the utility of a picket fence correction only in the
estimation of the frequency, when it isn't known exactly.

seniork wrote:

   ...

> pxx(left of f0) = 4.083
> pxx(f0) = 16.338
> pxx(right of f0) = 4.083
> 
> Notice that 16 - 4 = 12 is not 6 dB.  What am I missing? Argh!

16.338/4.083 = 4.001

A power ratio of of 4 is 6.02 dB.

   ...

Jerry
-- 
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Ah, yes, delta means ratio in this case.  I can never get used to dB.

Can someone explain the preferred method to convert the
windowed-periodogram estimates Pxx(f) into amplitudes?  Is the preferred
method to apply some window-dependent factor to the peak-value Pxx(f0) so
as to determine power *at* f0 and then from this determine the amplitude? 
Or, is it more appropriate to integrate the Pxx around f0 and convert this
sum into amplitude?  Let's assume we're only talking about the
deterministic case here, no noise.  

I would also like to know what it means when someone says to, "scale the
averaged frequency spectrum in terms of mean squared amplitude."  That is,
given the {Pxx(f)} where I have already made the Parseval scaling
sum(Pxx(f))=sum(x.^2)/length(x), what does the quoted sentence mean to do
to the {Pxx(f)}?

On Aug 15, 11:11 am, "seniork" <seni...@gmail.com> wrote:
> Ah, yes, delta means ratio in this case.  I can never get used to dB.
>
> Can someone explain the preferred method to convert the
> windowed-periodogram estimates Pxx(f) into amplitudes?  Is the preferred
> method to apply some window-dependent factor to the peak-value Pxx(f0) so
> as to determine power *at* f0 and then from this determine the amplitude?
> Or, is it more appropriate to integrate the Pxx around f0 and convert this
> sum into amplitude?  Let's assume we're only talking about the
> deterministic case here, no noise.
>

> I would also like to know what it means when someone says to, "scale the
> averaged frequency spectrum in terms of mean squared amplitude."  That is,
> given the {Pxx(f)} where I have already made the Parseval scaling
> sum(Pxx(f))=sum(x.^2)/length(x), what does the quoted sentence mean to do
> to the {Pxx(f)}?

One way to interpret this phrase is discussed on pp 30-38 of:
http://www.bksv.com/pdf/Bv0031.pdf

That would be 'mean squared amplitude' vs 'power spectral density' vs
'transient energy' interpretations.

Dale B. Dalrymple
http://dbdimages.com