www.gusucode.com > signal 案例源码程序 matlab代码 > signal/PeriodogramExample.m

    %% Periodogram
% In general terms, one way of estimating the PSD of a process is to simply
% find the discrete-time Fourier transform of the samples of the process
% (usually done on a grid with an FFT) and appropriately scale the
% magnitude squared of the result. This estimate is called the
% _periodogram_.
%
% The periodogram estimate of the PSD of a signal $x_L(n)$ of length _L_ is
%
% $$ P_{xx}(f) =
% {1\over{LF_s}}\left|\sum_{n=0}^{L-1}x_L(n)e^{-j2\pi fn/F_s}\right|^2, $$
%
% where Fs is the sampling frequency.
%
% In practice, the actual computation of $P_{xx}(f)$ can be performed only
% at a finite number of frequency points, and usually employs an FFT. Most
% implementations of the periodogram method compute the $N$-point PSD
% estimate at the frequencies
%
% $$ f_k={{kF_s}\over N},\;\;\;k=0,1,\dots,N-1. $$
%
% In some cases, the computation of the periodogram via an FFT algorithm is
% more efficient if the number of frequencies is a power of two. Therefore
% it is not uncommon to pad the input signal with zeros to extend its
% length to a power of two.
%
% As an example of the periodogram, consider the following 1001-element
% signal |xn|, which consists of two sinusoids plus noise:

%%

fs = 1000;                % Sampling frequency
t = (0:fs)/fs;            % One second worth of samples
A = [1 2];                % Sinusoid amplitudes (row vector)
f = [150;140];            % Sinusoid frequencies (column vector)
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
% The three last lines are equivalent to
% xn = sin(2*pi*150*t) + 2*sin(2*pi*140*t) + 0.1*randn(size(t));

%%
% The periodogram estimate of the PSD can be computed using |periodogram|.
% In this case, the data vector is multiplied by a Hamming window to
% produce a modified periodogram.

[Pxx,F] = periodogram(xn,hamming(length(xn)),length(xn),fs);
plot(F,10*log10(Pxx))
xlabel('Hz')
ylabel('dB')
title('Modified Periodogram Power Spectral Density Estimate')

%%
% *Algorithm*
%
% Periodogram computes and scales the output of the FFT to produce the
% power vs. frequency plot as follows.
%
% # If the input signal is real-valued, the magnitude of the resulting FFT
% is symmetric with respect to zero frequency (DC). For an even-length FFT,
% only the first (1 + |nfft|/2) points are unique. Determine the number of
% unique values and keep only those unique points.
% # Take the squared magnitudes of the unique FFT values. Scale the squared
% magnitudes (except for DC) by $2/(F_sN)$, where _N_ is the length of
% signal prior to any zero padding. Scale the DC value by $1/(F_sN)$.
% # Create a frequency vector from the number of unique points, the nfft
% and the sampling frequency.
% # Plot the resulting magnitude squared FFT against the frequency.