ch11.m

数字信号处理中功率谱的重要算法。代码中对比了经典谱估计直接法和间接法之间的关系

% The demonstration of relationship between Periodogram and
% Blackman-Turkey.
clear all;
N_fft=8;% length of signal.
x=[ 1 4 3 6 2 4 7 8];
y = [x  zeros(1,N_fft)];% signal with N_fft 0.

% PSD estimated by Periodogram.
subplot(431);
stem(x);xlim([1 16]);
title('Signal x')
X = fft(x,N_fft);
Pxx = X.* conj(X);
subplot(432);
stem(Pxx);xlim([1 16]);
title('|X_N(k)|^2');
subplot(433);
stem(Pxx/N_fft);xlim([1 16]);
title('|X_N(k)|^2/N');

subplot(434);
stem(y);xlim([1 16]);
title('x+0')

% PSD estimated by Periodogram with N_fft zeros at the end of signal x.
Y = fft(y,2*N_fft);
Pyy = Y.* conj(Y);
subplot(435);
stem(Pyy);xlim([1 16]);title('|X_2_N(k)|^2');
subplot(436);
stem(Pyy/N_fft);xlim([1 16]);title('|X_2_N(k)|^2/N');

% Comparison of totla power.
E_x=sum(x.*x);
E_X=sum(Pxx)/N_fft;
E_Y=sum(Pyy)/(2*N_fft);

% PSD estimated by Blackman-Turkey method.
% 2N_fft points. (two side)
Rx=xcorr(x)/N_fft;
Rx_fft=[Rx(N_fft:2*N_fft-1) Rx(1:N_fft)];
RX=fft(Rx_fft,2*N_fft);
subplot(437);
stem((-(N_fft-1):(N_fft-1)),Rx);xlim([-(N_fft-1) (N_fft-1)]);title('r(m)');
subplot(438);
stem(Rx_fft);xlim([1 16]);title('r_0(m)');
subplot(439);
stem(abs(RX));xlim([1 16]);title('P_B_T^2^N');

% PSD estimated by Blackman-Turkey method.
% N_fft points. (one side)
Rx_half=Rx(N_fft:2*N_fft-1);
RX_half=fft(Rx_half,N_fft);
subplot(4,3,10);
stem((-(N_fft-1):(N_fft-1)),Rx);xlim([-(N_fft-1) (N_fft-1)]);title('r(m)');
subplot(4,3,11);
stem(Rx_half);xlim([1 16]);;title('r(m):m\geq0');
subplot(4,3,12);
stem(2*real(RX_half)-Rx_half(1));xlim([1 16]);title('P_B_T^N');