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.