Differences between power spectrum and power spectral density
Dan
I was told to measure the power spectral density of a signal, i would like to know the advantages of using power spectral density over the others.
Moreover, i would like to know the differences between the power spectrum and power spectral density?
Thanks
Rune Allnor
> I am currently working on some kind of spectral analysis
> I was told to measure the power spectral density of a signal, i would like to know the advantages of using power spectral density over the others.
Which ones?
> Moreover, i would like to know the differences between the power spectrum and power spectral density?
It's only terminology. "Power Spectral Density" is the academically
correct term, whereas "Power Spectrum" is what one uses in day-to-day
language. Like somebody is named 'Daniel' in official documents but
go by 'Dan' among the friends.
Rune
Dan
By the way, can you suggest me some ways to get the average value before getting the power spectral density?
Some materials just mention the averaging can help a lot on getting a much accurate result of it
Thanks,
Dan
Godzilla
I would suggest you read about spectral analysis before you start using 'black boxes' that Matlab provides.
Greg Heath
> I am currently working on some kind of spectral analysis
> I was told to measure the power spectral density of a signal,
> i would like to know the advantages of using power spectral density over the others.
What others?
> Moreover, i would like to know the differences between the
> power spectrum and power spectral density?
The spectrum is the plot of the PSD vs frequency.
See my post on Power
http://groups.google.com/group/comp.soft-sys.matlab/msg/2222327db2ea7f51
Hope this helps.
Greg
dbd
Appropriate units to apply to your calculations may vary with the
nature of your signals.
Take a look at "Choose your Units!" from B&K:
http://www.bksv.com/doc/bo0438.pdf
Dale B. Dalrymple
Val Schmidt
> It's only terminology. "Power Spectral Density" is the academically
> correct term, whereas "Power Spectrum" is what one uses in day-to-day
> language. Like somebody is named 'Daniel' in official documents but
> go by 'Dan' among the friends.
I might be wrong, but I am pretty sure this is completely not true. Power Spectral Density has units of Power/Hz. But Power Spectrum has units of Power. This is how I think it should look in MATLAB
x = signal;
fs = sample rate;
N=length(x);
T = 1/fs * N; % recording duration
X = fft(x)/N;
psd = abs(X).^2;
fo = 1/T; % frequency bin size of the fft
ps =psd * 1/T; % power spectrum
You can check this by sticking a sign wave in for x with a period, T, different than 1 second. Then compare the peak value of the Power Spectrum (ps) with the mean square power value measured on the time series like this:
mean(x.^2)
Val Schmidt
> http://groups.google.com/group/comp.soft-sys.matlab/msg/2222327db2ea7f51
I am curious, could there be an error on this post in calculating the Power Spectral Density? From the post:
PSD = abs(X).^2/N; % Power Spectral Density
But I think this should be:
PSD = abs(X).^2/N^2; % Power Spectral Density
I might be wrong...
-Val
Martin Trauth
I found the Sign Proc Tbx Manual quite confusing since they mix different definitions and normalizations of the power spectral density while describing different methods for spectral analysis, in particular the periodogram method. Moreover, the corresponding MATLAB functions are not easy to read since they interfinger each other in a very complex manner, some of them are from the 80's, some of them have recently experienced the transition to object-oriented functions. In the early 90's as I was a PhD student I used SPECTRUM, later PSD and CPSD, then PERIODOGRAM became very popular that it often used in my field, now we have the SPECTRUM.ESTMETHOD set of functions. Thus, I can understand why people are confused a bit while trying do understand what a powerspectral density is and how it is computed with MATLAB.
As the author of a popular textbook on the use of MATLAB in Earth Sciences, I am currently rewriting the chapter on time series analysis. Reading through that chapter of the 2nd Edition you will see that even I was confused by the MATLAB set of functions for computing the PSD. The third edition will see a major rewrite of this chapter. In many other examples, I rather use simple MATLAB code to teach students how complex methods such as the Lomb-Scargle Algorithm works - here I translated the original FORTRAN code of Scargle into MATLAB rather than using the read-to-use functions available on the web.
Well, here is what I did yesterday that compares a simple MATLAB code to compute the PSD from the normalized Fouriertransform of the autocorrelation sequence and compare this with the PSD that comes out of the MATLAB function PERIODOGRAM. Any thoughts on this from the friends of the PSD around here?
With kind regards, Martin
-----------------------------------------------------------------------------------------------------
% Choose sampling frequency
Fs = 8;
% Define time vector and generate synthetic signal
t = 1/Fs :1/Fs : 1000/Fs;
x = 2*sin(2*pi*t/50) + sin(2*pi*t/15) + 0.5*sin(2*pi*t/5);
% Define nfft that is the next power of two close to the length of the
% signal
nfft = 2^nextpow2(length(t));
% Compute autocorrelation sequence of x
Corrxx = corrmtx(x,0);
% Compute Fourier transform of x using nfft
Xxx = fft(Corrxx,nfft);
% Frequency vector
f = 0 : Fs/(nfft-1) : Fs;
% Powerspectral density scaled by sampling frequency Fs and multiplied by
% two since we compute the twosided spectrum using the FFT but the function
% periodogram uses the oneside spectrum in the Nyquist frequency range
Pxx = 2 * abs(Xxx).^2 /Fs;
% Plot the Blackman-Tukey powerspectral density estimate
subplot(1,2,1)
plot(f,Pxx), grid
axis([0 0.5 0 max(Pxx)])
% Compute the onesided periodogram using nfft and Fs; you can choose the
% option 'twosided' to see the difference between the onesided and twosided
% powerspectral density estimate
[Pxx,f] = periodogram(x,[],nfft,Fs);
% Plot the Periodogram powerspectral density estimate
subplot(1,2,2)
plot(f,Pxx), grid
axis([0 0.5 0 max(Pxx)])
xlabel('Frequency')
ylabel('Power')
title('Autospectrum')
Martin Trauth
> See my post on Power
> http://groups.google.com/group/comp.soft-sys.matlab/msg/2222327db2ea7f51
> Unfortunately, the help files of pwelch and periodogram
> do not indicate that the corresponding PSDs are scaled
> by a factor of Fs, i.e.,
> PSDw = (1/Fs)*pwelch(x,ones(N,1),0,N,Fs,'twosided')
> PSDp = (1/Fs)*periodogram(x,ones(N,1),N,Fs,'twosided')
Page 2-10 of the Sign Proc Tbx User's Guide says that periodogram is scaled by Fs, and since Welch's method is based on periodogram it is also scaled by Fs. Am I wrong with this? See my other post to this thread.
Cheers ... Martin
Greg Heath
TideMan
All this worry about scaling is unnecessary if you approach the
problem from the other end.
Parseval's Theorem says that the area under the spectrum must equal
the variance.
For a spectrum S(f) at intervals of freq df, derived in whatever way
you like, the area is trapz(S)*df
Therefore, if the original signal was y(t), you must scale S(f) by var
(y)/(trapz(S)*df).
Doing it this way, you don't have to worry about ensemble or frequency
averaging or dividing by N, or sqrt(N), or whatever.
dbd
> On Mar 24, 1:22 pm, "Martin Trauth" <tra...@geo.uni-potsdam.de> wrote:
>
> > Greg Heath <he...@alumni.brown.edu> wrote in message
> > > See my post on Power
> > >http://groups.google.com/group/comp.soft-sys.matlab/msg/2222327db2ea7f51
> > > Unfortunately, the help files of pwelch and periodogram
> > > do not indicate that the corresponding PSDs are scaled
> > > by a factor of Fs, i.e.,
> > > PSDw = (1/Fs)*pwelch(x,ones(N,1),0,N,Fs,'twosided')
> > > PSDp = (1/Fs)*periodogram(x,ones(N,1),N,Fs,'twosided')
>
> > Page 2-10 of the Sign Proc Tbx User's Guide says that periodogram is scaled by Fs, and since Welch's method is based on periodogram it is also scaled by Fs. Am I wrong with this? See my other post to this thread.
>
> > Cheers ... Martin
>
.>
.> http://en.wikipedia.org/wiki/Spectral_density
>
> Hope this helps.
>
> Greg
Thanks for the warning Greg. It is a pretty shoddy presentation. It
illustrates the continuous formulas but fails to correctly relate the
discrete applications (any use of DFT) to the continuous forms. The
units for power spectral density and energy spectral density are given
but power is absent. The dependencies of the choice of unit on the
type of the signal providing the energy: stationary, stochastic for
psd, transient for esd and tonal for power, are absent.
Your warnings help
Dale B. Dalrymple
Martin Trauth
> averaging or dividing by N, or sqrt(N), or whatever.
This is very true, and probably important for those of us who are warned by Greg everyday to read the textbooks about spectral analysis. The posts, however, were more about the definition of the PSD in the MATLAB routines. And these are not easy to find in the manuals and in the MATLAB codes.
Just type "edit periodogram" and count the minutes how long it takes you to find the scaling of the PSD, as Greg mentioned. I am not a MATLAB fanatics, not an applied mathematician, I am an earth scientist who is much more interested in detecting major climate shifts during the past five million years that might have influenced human evolution than in exploring MATLAB code.
I like MATLAB as a tool, as it saves time compared to my early days as a researcher in 1992 when I used the brilliant book "Numerical Recipes in FORTRAN" and tried to compile the FORTRAN code for a FFT. But whoever R. Losada is who wrote the code and manuals for the spectrum tools, he/she is not doing a good job here. Just see line 91 in periodogram and you know what I mean. Is this how the visible code of a commercial tool should look like?
Kind regards, Martin
PS I am not using Windows :o)
TideMan
wrote:
> I am an earth scientist who is much more interested in detecting major climate shifts during the past five >million years that might have influenced human evolution than in exploring MATLAB code.
And I am an engineer whose interest is similar, but using sea-level
records (on much shorter timescales!).
But, given your interest, shouldn't you be using wavelet analysis, not
Fourier analysis? Using Fourier analysis removes all variability with
time. Indeed, a basic assumption for Fourier analysis is non-
stationarity, yet you are looking for changes with time. With wavelet
analysis, you can retain some variability with time, by doing a
partial transformation into frequency space. This allows you to
determine how a portion of the signal (over a range of timescales)
varies with time.
Martin Trauth
> records (on much shorter timescales!).
sounds interesting!
> But, given your interest, shouldn't you be using wavelet analysis, not
> Fourier analysis? Using Fourier analysis removes all variability with
> time.
Thanks for the suggestion! I recently worked on marine dust records to disprove the hypothesis of a drier climate that pushed human evolution. Since this work was about recalculating another person's analysis, we used Blackman-Tukey, Wavelets and Lomb-Scargle (the data are originally unevenly spaced) spectral methods besides lots of other methods to detect breakpoints in regression lines, transitions in time and frequency domain. As you know from analyzing sea-level records, we need to switch methods quite often depending on the types of data. Therefore, I am always happy to receive suggestions in this newsgroup. One of my former PhD student does recurrrence plots, another PhD student tries complex bayesian networks, which also seems to be great fun!
Well, back to our origins now ...
Greg Heath
> On Mar 25, 12:26 pm, Greg Heath <he...@alumni.brown.edu> wrote:
>
>
>
>
>
> > On Mar 24, 1:22 pm, "Martin Trauth" <tra...@geo.uni-potsdam.de> wrote:
>
> > > Greg Heath <he...@alumni.brown.edu> wrote in message
> > > > See my post on Power
> > > >http://groups.google.com/group/comp.soft-sys.matlab/msg/2222327db2ea7f51
> > > > Unfortunately, the help files of pwelch and periodogram
> > > > do not indicate that the corresponding PSDs are scaled
> > > > by a factor of Fs, i.e.,
> > > > PSDw = (1/Fs)*pwelch(x,ones(N,1),0,N,Fs,'twosided')
> > > > PSDp = (1/Fs)*periodogram(x,ones(N,1),N,Fs,'twosided')
>
> > > Page 2-10 of the Sign Proc Tbx User's Guide says that periodogram is scaled by Fs, and since Welch's method is based on periodogram it is also scaled by Fs. Am I wrong with this? See my other post to this thread.
>
> > > Cheers ... Martin
>
> > Also of interst:
>
> >http://en.wikipedia.org/wiki/Spectral_density
>
> > Hope this helps.
>
> > Greg
>
> All this worry about scaling is unnecessary if you approach the
> problem from the other end.
> Parseval's Theorem says that the area under the spectrum must equal
> the variance.
You are assuming zero mean?
Hope this helps.
Greg
dbd
> On Mar 25, 10:02 pm, "Martin Trauth" <tra...@geo.uni-potsdam.de>
> wrote:
>
> > I am an earth scientist who is much more interested in detecting major climate shifts during the past five >million years that might have influenced human evolution than in exploring MATLAB code.
>
> And I am an engineer whose interest is similar, but using sea-level
> records (on much shorter timescales!).
>
> But, given your interest, shouldn't you be using wavelet analysis, not
> Fourier analysis? Using Fourier analysis removes all variability with
.> stationarity, yet you are looking for changes with time. With
wavelet
.> analysis, you can retain some variability with time, by doing a
> partial transformation into frequency space. This allows you to
> determine how a portion of the signal (over a range of timescales)
> varies with time.
Tideman
Thank you for demonstrating the danger of mismatching the continuous/
infinite versions of Fourier analysis with discrete applications. No
one with real data has ever calculated the integral from -inf to +inf
in time because they are still waiting for their data collection to
complete. Calculations like the periodogram are applied to finite sets
of sampled data. Most interpretation of the periodogram relies on
signal changes being small during the data set of each DFT
calculation, but over longer times, signals may change greatly between
short finite sampled data sets. The application of 'Fourier
techniques' to finite sampled data sets with DFTs often appears in the
literature as the 'Short Time Fourier Transform (STFT)' which produces
a true time-frequency distribution.People who blithely mismingle
infinite/continuous/integral characteristics with finite/sampled/
summation applications or give references to that kind of shoddy
misrepresentation don't get anywhere near the abuse they deserve.
The wavelet analysis literature is filled with justifications based on
this kind of misrepresentation instead of comparisons to STFT
technique which would be useful, meaningful and appropriate.
Or maybe it's close to April Fools Day and I missed an attempt at
humor. Everyone can make their own call on that.
Dale B. Dalrymple
Steven Lord
"Martin Trauth" <tra...@geo.uni-potsdam.de> wrote in message
news:gqcrua$t8h$1...@fred.mathworks.com...
If you have a concern about whether the code is correct or a suggestion for
a way to improve the product, please let our Technical Support staff know so
they can pass your feedback on to Development. [They also track how
frequently users have contacted us about a particular problem/enhancement
and uses that information to help Development prioritize features and bug
fixes.] If you have a concern/suggestion about the documentation, you can
also give that feedback to Technical Support or you can go through the
feedback links in the upper- and lower-right corners of the documentation
pages (both on the website and in the Help Browser in sufficiently recent
releases of MATLAB.)
--
Steve Lord
sl...@mathworks.com
Martin Trauth
> If you have a concern ...
The code is probably correct but not well written, as the manual about it. I have contacted my German MathWorks representative in the morning already who forwarded my message directly to Ric Losoda who has written most of the code for spectral analysis. Also, I have contacted the MathWorks book program for assistance in rewriting the corresponding chapter in my textbook.
This discussion came up after I have checked the newsgroup for other people having problems with the definition of the power spectral density with MATLAB. I assumed that it was my problem but now it turns out to be other's problem as well. That's why I then emailed my German MathWorks contact and book program.
"Hope this helps", using Greg's words.
Regards, Martin
TideMan
dbd:
If people applied Fourier analysis only in the way you describe using
STFT, I would have no problem. Indeed, I use precisely this method
myself when calculating real-time ocean-wave spectra for my website
(http://www.mulgor.co.nz/MarsPt/WaveSpectra.html).
But there are continual examples given in this forum of people
applying the FFT to non-stationary time series like the monthly
Southern Oscillation Index (SOI). These are the ones that would be
better analysed using wavelets.
dbd
Tideman
This is a discussion of a STFT technique: the periodogram. I think
that comparing an infinite/continuous Fourier technique with infinite
time scale (which is not what the periodogram does) to a finite
duration tool is misleading in this context.
If people are applying a measurement tool with the wrong time scale
to their data, we should recommend a better time scale.
Whether the proper time scale for sub-band decomposition is a
constant or a function of frequency for an analyst's data might make
an interesting discussion topic for another thread comparing STFT and
wavelet implementations.
The discussion in this thread seems to be about determining whether
the periodogram function performs the proper calculation for someone's
data. That depends on whether the function can be used to correctly
calculate the PSD (I think it can) and whether the PSD is really the
proper tool.
Dale B. Dalrymple
ch li
> On Mar 1, 5:13=A0am, Dan <chunhoish...@yahoo.com.hk> wrote:
> > I am currently working on some kind of spectral analysis
> > I was told to measure the power spectral density of a signal,
> the others.
>
> What others?
>
> > Moreover, i would like to know the differences between the
> > power spectrum and power spectral density?
>
> The spectrum is the plot of the PSD vs frequency.
>
> See my post on Power
>
> http://groups.google.com/group/comp.soft-sys.matlab/msg/2222327db2ea7f51
>
> Hope this helps.
>
> Greg
Dear Greg,
Excuse me, in your post.
"PSD = abs(X).^2/N; % Power Spectral Density"
I don't realize this equation.
I think the abs(X) is the amplitude of bin in frequency domain.
why the power of each bin isn't abs(X)^2?
And thank the post is very useful for me.
lich
I have gathered some expert's perspective and derived the following
According to Greg Heath's post
http://groups.google.com/group/comp.soft-sys.matlab/msg/2222327db2ea7f51
and dbd's post
Take a look at "Choose your Units!" from B&K:
http://www.bksv.com/doc/bo0438.pdf
and Val Schmidt's post
Define N : number of FFT point; T : signal time domain period
then get dt = T/N, df = 1/T, Fs = 1/dt
We can derive some power related equations
For easy to read, I substitue x for magnitude of time domain signal,
and X for magnitude of frequency domain signal.
Meanwhile, the unit has marked behind the equation
Time Domain:
total signal energy Pxt = sum(dt*x^2) ........................................... (J) [1]
average signal power Pxtav = Pxt/T = sum (x^2*dt/T)
= sum (x^2*(T/N)/T) = sum (x^2/N) ..... (W) [2]
if we make Fourier Transform scaling (1/N) in IFFT synthesis eq.(as adopted in MATLAB)
then the Parseval's equality is derived as sum[x^2] = 1/N*sum[X^2]
and then take a look at average power definition
http://www.stanford.edu/class/ee179/summaries/lecture5.pdf
According to the average power is identical and Parseval equation above
we can get that
Frequency Domain:
average signal power Pxfav = sum[(X^2)/(N^2)] ............................ (W) [3]
then the total signal energy as
Pxf = Pxfav*T = sum[(X^2)/(N^2)*T]
= sum[(X^2) / (N^2) * (T^2) / T]
= sum[(X^2) / (N^2) * (T^2) * df] .......................... (J) [4]
From the eq. [4]
the (X^2) / (N^2) * (T^2) * df ...................................................... (J) [5]
can be thought as the energy of each FFT frequency bin.
then (X^2) / (N^2) * (T^2) ........................................................... (J/Hz) [6]
can be thought as the energy at the freq. sample over 1 Hz, (ESD)
further, we transfer the eq. [6] straight to physical meaning
according to Fs = 1/dt = 1/(T/N) = N/T
the ESD = (X^2)/(Fs^2) ................................................................ (J/Hz) [7]
So the PSD = ESD/T ...................................................................... (W/Hz) [8]
Assume the load 1 ohm is applied,
therefore, sqrt(ESD) = X/Fs is the frequency amplitude response over a Hz
and then X is the frequency amplitude response over Fs (the freq span we observed)
If we look back the IFFT we used (with 1/N scaling), and according to the X meaning.
We can get that X/N means the freq. amp. response over an FFT freq. bin.
Or we derive it from X/Fs (over a Hz)
X/Fs = X/N*T .............................................................................. (V/Hz) [9]
X/N = X/N*T/T = X/Fs/T = X/Fs*df (Amp. per Hz * An FFT bin freq. span) .. (V) [10]
Again, X/N means the freq. amp. response over an FFT freq. bin.
Therefore, in IFFT synthesis eq., the value of the freq. samples has to be normalized
to represent its own bin's power.
If we define X/N as the freq. amp. spectrum.
the energy spectrum is (X/N/df)^2*df = (X^2)/(N^2)/df .................... (J/bin) [11]
the power spectrum is (X^2)/(N^2)/df/T = (X^2)/(N^2) ................... (W/bin) [12]
(Note: the summation of energy spectrum is total energy of signal
the summation of power spectrum is average power of signal
refer to the above derived eq.)
About power spectrum explanation, please see the following
http://mathworld.wolfram.com/PowerSpectrum.html
Finally, if we try the others FT scaling method, we can get "X" in different explanation.
This explanation of scaling, power spectrum and power spectral density is
in reference to many experts' perspectives.
But I'm not sure the statement is right.
Please give some comments on unreasonable point.
Thanks in advance!
Lich
Martin Trauth
I have noticed a mistake in the code I have posted some time ago. I have mixed to recipes, one compared the Blackman-Tukey Autospectrum with Periodogram, the other one compared a step-by-step recipe for the periodogram with the MATLAB function periodogram. The BT spectrum is the FFT of the autocorrelation sequence. Here are two ways to compute the autocorrelation sequence with CORRMTX or XCORR.
clear, Fs = 1;
t = 1/Fs :1/Fs : 1000/Fs; t = t';
x = 2*sin(2*pi*t/50) + sin(2*pi*t/15) + 0.5*sin(2*pi*t/5);
[X,Corrxx1] = corrmtx(x,length(x)-1); Corrxx1 = Corrxx1(:,1);
Corrxx2 = xcorr(x,'biased'); Corrxx2 = Corrxx2((end+1)/2:end,1);
subplot(2,1,1), plot(Corrxx1)
subplot(2,1,2), plot(Corrxx2)
Below is the comparison of a the step-by-step periodogram with the function PERIODOGRAM. Comments are welcome.
Kind regards, Martin
-----
% Comparison of Periodogram step by step with with Periodogram using
% MATLAB function PERIODOGRAM. First with flexible sampling frequency
clear
% Choose sampling frequency
Fs = 1;
% Define time vector and generate synthetic signal
t = 1/Fs :1/Fs : 1000/Fs; t = t';
x = 2*sin(2*pi*t/50) + sin(2*pi*t/15) + 0.5*sin(2*pi*t/5);
% Define nfft that is the next power of two close to the length of the
% signal
nfft = 2^nextpow2(length(t));
% Compute Fourier transform of x using nfft. This is the two-sided Fourier
% spectrum of Corrxx
Xxx = fft(x,nfft);
% Two-sided power spectral density scaled by sampling frequency Fs and length of data series
Pxx2 = abs(Xxx).^2 /Fs /length(x);
% Use one-sided power spectral density only
Pxx = [Pxx2(1); 2*Pxx2(2:512)];
% Frequency vector
f = 0 : Fs/(nfft-1) : Fs/2; f = f';
% Plot the Blackman-Tukey powerspectral density estimate
subplot(1,2,1)
plot(f,Pxx), grid
axis([0 0.5 0 max(Pxx)])
% Compute the onesided periodogram using nfft and Fs; you can choose the
% option 'twosided' to see the difference between the onesided and twosided
% powerspectral density estimate
[Pxx,f] = periodogram(x,[],nfft,Fs,'onesided');
Harish
Hi...
The PSD using x and Corrxx1 will be same? Because when i used it, i
found the similar plot but difference in Values.
Kind regards
Harish
Rune Allnor
> The PSD using x and Corrxx1 will be same? Because when i used it, i
> found the similar plot but difference in Values.
Lots of people ignore some scaling factors involved
in spectrum estimates. For single-frame spectra,
try ans scale with N or sqrt(N) where N is the number
of data points (which might be different from the FFT
length.) You might also have to scale by another factor 2.
These factors don't matter unless you work with calibrated
analysis. Very few people do that.
Rune
Leo Schirber
Martin,
Thanks for the honesty about your confusion concerning PSD estimation in Matlab (and elsewhere): it's widespread in my opinion. I will look over your 2 scripts, although one you tell me is in error. Here are some thoughts that lead up to a basic question. It'd be easier if I had an equation editor, but I'll try to get by without. My question is: what are consistent units for the magnitude squared of the Fourier transformed signal? I'm going to use the comm signal convention that the squared magnitude of x of t is power in Watts.
Assume first the continuous formulation with "infinite duration" time signals. (Later perhaps we can think of "real world" signals that are always finite in duration.) Assume that the instantaneous power in Watts of the signal x(t) is found from the square of the absolute value of x(t): i.e., P_inst(t) = abs(x(t))^2. A signal x(t) has finite energy E in Joules if the integral over all time of P_inst(t) is finite. If a signal doesn't have finite energy, I can say that it has finite power if the limit of 1/T times the integral from 0 to T of P_inst(t) converges. Call that converged value P in Watts.
For an energy signal with energy E, define the "energy spectral density" or ESD as the distribution of energy in frequency space. We can define the scaling for this function in different ways, but one simple way is just to set the ESD equal to absolute value the FT of x(t) squared:
ESD(f) = abs( FT(x(t)) )^2 = abs( X(f) )^2
where X(f) = FT{x(t)} = Integral_all_time{x(t)*exp[-j*2*pi*f*t)] dt }
Here f is a frequency variable in Hertz. With the above definition for the Fourier transform, I can write (from Parseval's theorem) that the energy E in both domains is the same.
E = Integral_all_time{abs(x(t))^2 dt} = Integral_all_frequency{abs(X(f)^2 df}
From this last expression I infer that the units for abs(X(f))^2 must be W-s^2, because the product of that quantity and Hz must be Joules. Would you agree with this?
Leo
dbd
> ...
> ... My question is: what are consistent units for the magnitude squared of the Fourier transformed signal? ...
The only consistent answer is "It depends." Matlab provides tools to
calculate the PSD. This is appropriate for stocastic signals that have
a defined power spectral density. A deterministic signal does not have
a power spectral density. It has a power spectrum. The power spectrum
has different units than the power spectral density. An energy signal
usefully characterized by other units yet. The real question is not
"What are the units for the PSD?" (magnitude squared per Hertz), but
"What are the correct units to describe my signal?" Deterministic
signals have power spectra. Energy signals have energy spectral
densities. The "real world" problem here is that signals can be sums
of different types of signals and the fft coefficients can represent
different types of signals for each coefficient. (Worse, an fft
coefficient can represent the sum of component of differing signal
types.)
> ...
> Assume first the continuous formulation with "infinite duration" time signals. (Later perhaps we can think of "real world" signals that are always finite in duration.) ...
The disconnect between arguments such as this and the calculations
Matlab performs is not about "ideal" versus "real" but rather about
continuous/infinite versus sampled/finite. These are different
theoretical mathematical cases. The common, incorrect practice seems
to be to analyze in the continuous/infinite domain and convert the
"ideal" case to "real world" by assuming that the fft (really DFT) of
samples of the continuous/infinite input will be samples of the
continuous/infinite Fourier transform of that input. Unfortunately,
that is never the case for signals with defined PSDs. The repeated
diversion of this question to continuous/infinite analysis is not only
not helpful, but is actively harmful and common. That is why the
confusion is widespread and continuing.
Outside of symbolic manipulations, Matlab processes sampled/finite
data and meets that requirement quite usefully when not confused with
inapplicable infinite/continuous expectations. The PSD calculation
does not determine if your signal is deterministic, random or
transient in nature. It is possible to test for such characteristics
in Matlab, but such algorithms will be limited in their success by the
nature of sampled/finite arithmetic.
Dale B. Dalrymple
Leo Schirber
Say we're given a finite stream of time-varying, real signal data. We could imagine that it's a member of a random (or stochastic) process, but that doesn't do a lot for us (unless of course we've manufactured the data based on a model with known random noise statistics, for example.) It's definitely not deterministic. We want to describe its frequency content in some way, shape, or form. We turn to Matlab, and try to determine what it is doing with its current set of spectral estimation routines.
Matlab adds to the mystery by not assigning units to every variable, but I can and will. Suppose that I'm a comm engineer, and I want to find either the ESD or the PSD estimates for my signal data, x(t), and x(t) is in units of volts or V. I use the common association of instantaneous signal power to the absolute squared magnitude of x(t). My questions are: how can I use Matlab to estimate (a) the ESD or (b) the PSD? These quantities are usually quoted in units of J/Hz for ESD and W/Hz for PSD.
Lest you just reply with a code snippet with a "periodogram" or "pwelch" call, let's review some of the definitions and the difficulties with the definitions of PSD that you find in books. If I have an energy signal, then it seems to me that I'm "safe" in following the crowd and using methods that basically find for me a smoothed magnitude squared FT of my data. The way I'd put it, what I'm asking for in this case is more properly called the ESD. I also believe that this is what is actually calculated with the periodogram and pwelch functions, even though it is referred to as a PSD.
However, if I imagine that I actually have a power signal (imagination is required perhaps because I have to imagine that the signal stream is extended "forever" and never "dies out"), then I think there is the need to determine a "time-averaged autocorrelation function". (See, for example, Proakis and Salehi "Contemporary Communication Systems using Matlab). In symbols, this time fuction can be defined as
r_xx(tau) = lim T -> inf { 1/T * Integral_from -T/2 to T/2: {x(t) x(t+tau) dt} }
Then one can take the FT of r_xx to find the theoretical PSD.
So, does Matlab do this in its functions like "xcorr"? I think not - Matlab has a tendency to use discrete analogs and not calculate limits - , but i's approach is not known to me, however, so I'm asking.
Leo
dbd
> ...
> Lest you just reply with a code snippet with a "periodogram" or "pwelch" call, let's review some of the definitions and the difficulties ...
No snippet has touched my hands. Please see the discussions of the
processing of sampled/finite data below.
> ... Matlab has a tendency to use discrete analogs and not calculate limits - , but i's approach is not known to me, however, so I'm asking.
>
> Leo
You are correct that Matlab calculates in the sampled/finite domain.
Your reversion to integral equations will continue to fail to give
useful insight about sample/finite domain implementation.
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 you have asked about,
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.
Dale B. Dalrymple
Leo Schirber
Dale,
Thanks. I've taken a look at the references. In particular, the journal article looks very good for getting into the details of analyzing data.
However, I still have to object to Matlab's (and the journal article's) use of the term PSD, as I've hinted at before. This term has a perfectly good definition and a long history. See for example http://www.cse.unt.edu/~rakl/class3020/Chapter8.pdf.
Matlab is not calculating a PSD as it as it is classically defined, as in the above article. Agreed? A yes or no answer is requested.
Leo
dbd
> ...
> Matlab is not calculating a PSD as it as it is classically defined, as in the above article. Agreed? A yes or no answer is requested.
>
> Leo
Matlab is not calculating PSD as PSD is 'classically' CALCULATED. If
that's a yes for you, then yes.
For a rectangular window, the psd of a sequence padded with zeros to
double the length will produce the same result as that produced with
correlation processing. Without zero padding, the DFT based psd
generates every other sample of the correlation calculated psd. In the
DFT based non-zero-padded signal processing path, the DFT coefficients
just before the calculation of the modulus contain all the information
of the original data. For the DFT of the zero padded sequence, every
other coefficient is identical to the coeffient DFT from the non-zero
padded sequence. In fact, the DFT can be inverted to regenerate the
original data. This characteristic is used in some analysis/synthesis
applications. The fft also provides a more efficient implementation
than the direct form correlation based calculation. Of course, the
correlation based calculation can be implemented with increased
efficiency for large data sequences by using ffts of the data sequence
zero padded to twice the original length.
Efficiency and the ability perform synthesis form intermediate data
(DFT coefficients) has lead many, including Matlab the calculate psd
and periodogram via fft.
Dale B. Dalrymple
dbd
> dbd <d...@ieee.org> wrote in message <db9fb79c-b802-4449-8b30-0b0254c63...@q40g2000prg.googlegroups.com>...
Leo Schirber
Dale, I think your "yes" means "no" to my question. You're saying that Matlab is calculating a PSD correctly, but not via a limit.
However, after reading my communications books, I'd agree with you with the following reservation: I think there's nothing really "wrong" with Matlab's PSD calculation, apart from perhaps a somewhat cavalier (in my opinion) accounting of the units of the quantities in the description. Here's how I'd describe PSD estimation, paying very careful attention to the units for each step of the way. I'd appreciate your comments or criticism, if you can get through all this muck. My reference is [1], Stremler, "Introduction to Communication Systems" 1982.
Let's start with a power signal x(t), a continuous function of time measured in volts [V] and defined for all time t in seconds. We make a sequence of N measurements of x(t) separated at equal intervals of time delta_t to create vector x[n]. Assume without loss of generality that the starting time for the measurements is t = 0, and call T = N*delta_T the "span" of the measurements in seconds. Let x[n] be a short-hand notation for the nth measurement at t_n = (n-1)*delta_t, i.e., x[n] = x(t_n) for n = 1,2, ... N.
The problem is to estimate the power spectral density (PSD) of x. Call the discrete estimate S(K), which is a shorthand for S(f_K), f a frequency variable in Hz. Assume for simplicity that x(t) is band-limited with cutoff frequency f_c in Hz, and that f_c < f_Nyquist = 1/(2*delta_t).
Recall the following definition from [1] of the PSD of x , S_x(f), or more simply, just S(f):
S(f) = limit_as_W_approaches_inf { (1/W)* | X_W(f) |^2 ) ... Units [W/Hz] (1)
Here X_W(f) is the Fourier transform of continuous function x(t) truncated at t = W. Explicitly,
X_W(f) = integral_from_0_to_W { x(t)*exp(-j*2*pi*f*t) dt } ... Units [V-s] (2)
Alternately, X_W(f) could be considered as the FT of the product of two time functions: a "rectangular-window" function (1 on [0,W], 0 elsewhere) and our original time function x(t). Hence we know that X_W(f) can be also be expressed as the convolution of the continuous FT of x(t) ( X(f) ) and a "sinc" ( sin(pi*f)/f ) function.
We make a first and rough estimate of S_x(f) from our finite data vector x[n] as follows. Define a "truncated estimate" of S(f) over [0,T] and denote it by S_T(f):
S_T(f) = (1/T) * | X_T(f) |^2 ... Units [W/Hz) (3)
We expect (actually hope is more accurate) for T "large" that S_T(f) will be a fair approximation to S(f), as we know by definition that the limit as T approaches infinity of (3) is (1).
Next we approximate (3) numerically by approximating the FT by a DFT. Consider X[K] as the (slightly nonstandard) expression for the DFT of x[n] as below:
X[K] = delta_t * sum_from_n=1_to_N { x[n]*exp(-j*2*pi*(n-1)*(K-1)/N) } for 1 <= K <= N
where X[K] is approximately equal to X(f_K), X(f) the continuous FT of x(t),
with f_K = (K-1)*delta_f, delta_f = 1/(N*delta_t), f_K in [Hz]
This expression is nonstandard only because an extra factor of delta_t appears in front of the sum: the reason for this factor is to conserve energy in the frequency and time domains. A heuristic explanation is that approximating the Fourier integral (integral over time of a voltage function) by a finite sum necessarily requires an extra factor of "delta_time" to appear out in front of the sum of voltages.
So using this expression for the DFT of x[n] - and agreeing with the association of K and frequency vector f_K as above for S_T - we can write for our rough and discrete approximation for S_x(f):
S_T[K] = (1/T) * | X[K] |^2 ; 1 <= K <= N ... Units [W/Hz] (4)
In (4) we have followed the somewhat confusing but common convention of allowing our frequencies to range between f(1)= 0 and f(N) = 2*f_Nyquist - delta_f, with the understanding that the frequencies above f_Nyquist are the "folded-over" negative frequencies, just as is conventionally done in X[K]. (A "one-sided" representation for S would add the squared magnitudes of the folded over values to the values between 2 and N/2, 1 being the DC term)
To refine this rough estimate, consider breaking up the length N vector into contiguous sub-sequences, say L sub-sequences of length M where N = L*M. We can then evaluate L expressions similar in form to (4) and find their average. This algorithm comprises "Bartlett's method" of spectral estimation.
For Welch's method - a variation on Bartlett's method - , a "window" is applied to the subsequence of length M, and the sub-sequences are also allowed to overlap. A window in this context is a real-valued time function on interval (a,b) that (generally) is maximized and somewhat flat near (a+b)/2 and drops to zero at t = a and t = b. Popular windows are named after Hamming, Hanning, Hann, Blackman, and others. Again, however, L expressions similar to (4) are found, and their average is computed to estimate S[K] and hence approximate S_x(f).
A couple of final comments:
(1) In Matlab, with the delta_t factor absent in the DFT, energy is not conserved in time and frequency domains. Admittedly, one can always normalize the squared magnitude of the FT result so that energy is conserved.
(2) The expression in (4) has the correct units of [W/Hz] if and only if the factor of 1/T is included; if omitted the units would be [J/Hz]. That was my point. That factor of 1/T is included in "Bartlett's method" for spectral estimation (see for example, Wikipedia on Bartlett's method), but not usually mentioned in describing "Welch's method" (in for example, Matlab). Admittedly, after the window is applied the energy would need to be adjusted also, so the factor would have to be adjusted from 1/T due to the windows effect.
Leo
dbd
>
> However, after reading my communications books, I'd agree with you with the following reservation: I think there's nothing really "wrong" with Matlab's PSD calculation, apart from perhaps a somewhat cavalier (in my opinion) accounting of the units of the quantities in the description. Here's how I'd describe PSD estimation, paying very careful attention to the units for each step of the way. I'd appreciate your comments or criticism, if you can get through all this muck. My reference is [1], Stremler, "Introduction to Communication Systems" 1982.
>
>
>
> Leo
If you want to understand what Matlab does, get a reference that
provides the discrete /finite forms of the laws and transforms you
are interested in. As soon as there is sampling the integral forms are
no longer applicable. They work where they belong, in theoretical
development and homework problems. Don't expect wikis to understand
the difference.
I took the time to explain once before:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dbd <d...@ieee.org> wrote in message <4e5bab06-88c7-424c-bc6a-
ec77d135e...@n26g2000yqh.googlegroups.com>...
...
The disconnect between arguments such as this and the calculations
Matlab performs is not about "ideal" versus "real" but rather about
continuous/infinite versus sampled/finite. These are different
theoretical mathematical cases. The common, incorrect practice seems
to be to analyze in the continuous/infinite domain and convert the
"ideal" case to "real world" by assuming that the fft (really DFT) of
samples of the continuous/infinite input will be samples of the
continuous/infinite Fourier transform of that input. Unfortunately,
that is never the case for signals with defined PSDs. The repeated
diversion of this question to continuous/infinite analysis is not only
not helpful, but is actively harmful and common. That is why the
confusion is widespread and continuing.
Outside of symbolic manipulations, Matlab processes sampled/finite
data and meets that requirement quite usefully when not confused with
inapplicable infinite/continuous expectations. The PSD calculation
does not determine if your signal is deterministic, random or
transient in nature. It is possible to test for such characteristics
in Matlab, but such algorithms will be limited in their success by the
nature of sampled/finite arithmetic.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Techniques such as Welch's can be used to calculate both the power
spectral density, energy spectral density and the power spectrum. You
must decide which correctly characterizes your signal and "fix the
units" yourself.
Dale B. Dalrymple
Leo Schirber
Dale,
I understand your main point, and realize you've made it before in this thread, and maybe a time or two before this. However, I disagree with your "philosophy" and here's why.
You tell me (and others) that there is a "dichotomy" between the infinite/continuous and the finite/sampled regimes. Well, yes. You disparage the attempts (often these attempts "fail" it's true) to explain the finite/sampled algorithms in terms of the infinite/continuous theory. Why would people continue on these incorrect paths when they apparently only lead to confusion and unhappiness?
The answer is: the definitions of the quantities or functions we want to estimate are given in terms of the infinite/continuous. The Fourier transform is defined as an integral over all time, and we actually want to estimate the FT of time functions that do not have finite duration. The PSD for a power signal x(t) is defined as the limit of the autocorrelation of a "windowed" x(t) as the window increases in duration without limit.
Yes, once we go to the computer we have to deal with the reality of sampled data and finite arrays. But that doesn't change the definitions. Algorithm developers have to "deal" with those definitions, and do the best they can to formulate approximations and "work-arounds". Successful algorithms can approximate the continuous definitions, with restrictions.
In my opinion, in order for a user to fully understand one of these algorithms involving infinite/continuous functions, the user needs to see how they are connected to the infinite/continuous definitions. Next a full disclosure of the assumptions, approximations, work-arounds, etc., that apply must be given. (Very helpful here is an accounting of the units of the quantities involved for an example calculation.) Only then can the user hope to understand the finite/sampled algorithm.
Finally, supplying only the finite/sampled details of an algorithm that attacks an infinite/continuous problem - no matter how precisely - without at least a "nod" to the infinite/continuous underlying ideas seems to me to be just plain unacceptable. No one will ever be able to figure out what is being done or why it is being done from that narrow description. Use of the algorithm is bound to lead to difficulties, because of a lack of understanding about the assumptions and where the approximations are.
Leo
dbd
>
> You tell me (and others) that there is a "dichotomy" between the infinite/continuous and the finite/sampled regimes. Well, yes. You disparage the attempts (often these attempts "fail" it's true) to explain the finite/sampled algorithms in terms of the infinite/continuous theory.
I think that to use one regime to approximate the other you need to
understand both. There exists a body of mathematics that deals with
the discrete finite regime that requires no basis in or requirement
for an analogy to the FT. I think you should try to link the regimes
by understanding both and comparing the analogous parts to understand
the differences. I chose not to participate in continuing discussion
of the infinite/continuous regime because I'm not convinced it is on
topic in comp.soft-sys.matlab and because I feel sure that such
discussions would only encourage some in their "failing attempts"
You'll have to find someone else to discuss that regime with. I have
repeatedly tried to get you to include the nature of the digital
regime in your analysis and you haven't.
> Why would people continue on these incorrect paths when they apparently only lead to confusion and unhappiness?
>
Because they were taught by educators who did not understand the
existence or nature of the two regimes.
> The answer is: the definitions of the quantities or functions we want to estimate are given in terms of the infinite/continuous. The Fourier transform is defined as an integral over all time, and we actually want to estimate the FT of time functions that do not have finite duration. ...
>
Agreed.
> Yes, once we go to the computer we have to deal with the reality of sampled data and finite arrays. But that doesn't change the definitions. Algorithm developers have to "deal" with those definitions, and do the best they can to formulate approximations and "work-arounds". Successful algorithms can approximate the continuous definitions, with restrictions.
>
> In my opinion, in order for a user to fully understand one of these algorithms involving infinite/continuous functions, the user needs to see how they are connected to the infinite/continuous definitions. Next a full disclosure of the assumptions, approximations, work-arounds, etc., that apply must be given. (Very helpful here is an accounting of the units of the quantities involved for an example calculation.)
You give lip service to the existence of the discrete/finite regime,
but your use of "full disclosure" and "work-arounds" shows an attitude
consistent with an attempt to explain discrete outputs without taking
the time to compare the fundamental definitions in the discrete regime
with their continuous counterparts.
>Only then can the user hope to understand the finite/sampled algorithm.
You cannot learn the characteristics of finite/sampled algorithms from
the definitions of physical parameters in the infinite/continuous
regime.
>
> Finally, supplying only the finite/sampled details of an algorithm that attacks an infinite/continuous problem - no matter how precisely - without at least a "nod" to the infinite/continuous underlying ideas seems to me to be just plain unacceptable. No one will ever be able to figure out what is being done or why it is being done from that narrow description. Use of the algorithm is bound to lead to difficulties, because of a lack of understanding about the assumptions and where the approximations are.
>
> Leo
My philosophy is that there are two mathematical regimes with equal
logical standing. If you want to use one to approximate the other, you
need to understand both so that you can compare and understand the
differences.
Trying to leap from a front end in one regime to a back end in the
output of the other is the usual method used in comp.soft-sys.dsp to
avoid the need to understand and come to terms with the nature of the
finite sampled regime.
Dale B. Dalrymple
Martin Kovac
I noticed that there has launched a heated debate on the PSD topic so I hope that you could help me. I try to find PSD of Gaussian monopulse by two methods-periodogram and FFT but I get different results in amplitude. I wonder why? Please, this is the matlab code:
% GAUSSIAN PULSE
fc = 2e9; % Center frequency
tc = gmonopuls('cutoff',fc); % Width of each pulse
fs=800*fc; % Sampling frequency
SegmentDuration = 2*tc;
N = fs*SegmentDuration;
t = -SegmentDuration:(1/fs):(SegmentDuration-1/fs); % Signal evaluation time
y = gmonopuls(t,fc);
figure();
plot(t,y)
figure();
periodogram(y, rectwin(length(y)), [], fs, 'oneside', 'power');
hgcf = gcf;
hgcf.Color = [1,1,1];
mean_power = mean(y.^2);
mean_power_dBm = 10*log10(1000*mean_power);
%-----------------------------------------------------------------------------------------------------------
rng default
fc = 2e9;
Fs = 1600e9; % Sampling frequency
t = 0:1/Fs:(4*tc-1/Fs); % Signal evaluation time
tc = gmonopuls('cutoff',fc); % Width of each pulse
x = pulstran(t,2*tc,@gmonopuls,fc);
figure();
plot(t,x)
N = length(x);
NFFT = 2^nextpow2(N);
xdft = fft(x,NFFT);
xdft = xdft(1:NFFT/2+1);
psdx = (1/(Fs*NFFT)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/NFFT:Fs/2;
figure();
plot(freq,10*log10(psdx))
grid on;
title('Periodogram Using FFT');
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');
Thank you for answers.