現代數字信號處理課后作業【完結篇】

感謝大家的陪伴，學無止境，

1、設各態歷經平穩隨機信號的觀測值 x ( n ) = x(n)= x(n)={1，-1，2，-2；n=0~3}；試用周期法和相關法求其功率譜估計值 P N ( w ) P_N(w) PN?(w)，驗證其一致性，

周期圖法：
①先求信號頻譜函式 X N ( e j w ) = ∑ n = 0 N ? 1 x ( n ) e ? j w n X_N(e^{jw})=\sum\limits_{n=0}^{N-1}x(n)e^{-jwn} XN?(ejw)=n=0∑N?1?x(n)e?jwn
②再求功率譜密度 P N ( w ) = 1 N ∣ X N ( e j w ) ∣ 2 P_N(w)=\dfrac{1}{N}|X_N(e^{jw})|^2 PN?(w)=N1?∣XN?(ejw)∣2

自相關法：
①首先在時域求序列自相關函式 r N ( m ) = 1 N x N ( m ) ? x N ? ( ? m ) = 1 N ∑ n = 0 N ? 1 x N ( n ) x N ( n ? m ) r_N(m)=\dfrac{1}{N}x_N(m)*x^*_N(-m)=\dfrac{1}{N}\sum\limits_{n=0}^{N-1}x_N(n)x_N(n-m) rN?(m)=N1?xN?(m)?xN??(?m)=N1?n=0∑N?1?xN?(n)xN?(n?m)
②再求其傅里葉變換得到功率譜密度 P N ( w ) = ∑ m = 1 ? N N ? 1 r N ( m ) e ? j w n P_N(w)=\sum\limits_{m=1-N}^{N-1}r_N(m)e^{-jwn} PN?(w)=m=1?N∑N?1?rN?(m)e?jwn

周期圖估計法：

x 4 ( n ) = { 1 , ? 1 , 2 , ? 2 } N = 4 x_4(n)=\{1,-1,2,-2\} \ \ \ \ \ N=4 x4?(n)={1,?1,2,?2} N=4

X 4 ( e j w ) = ∑ n = 0 3 x ( n ) e ? j w n = 1 ? e ? j w + 2 e ? j 2 w ? 2 e ? j 3 w X_4(e^{jw})=\sum\limits_{n=0}^{3}x(n)e^{-jwn}=1-e^{-jw}+2e^{-j2w}-2e^{-j3w} X4?(ejw)=n=0∑3?x(n)e?jwn=1?e?jw+2e?j2w?2e?j3w

= ( 1 ? c o s w + 2 c o s 2 w ? 2 c o s 3 w ) + j ( s i n w ? 2 s i n 2 w + 2 s i n 3 w ) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =(1-cosw+2cos2w-2cos3w)+j(sinw-2sin2w+2sin3w) =(1?cosw+2cos2w?2cos3w)+j(sinw?2sin2w+2sin3w)

P 4 ( w ) = 1 4 [ ( 1 ? c o s w + 2 c o s 2 w ? 2 c o s 3 w ) 2 + ( s i n w ? 2 s i n 2 w + 2 s i n 3 w ) 2 ] P_4(w)=\dfrac{1}{4}[(1-cosw+2cos2w-2cos3w)^2+(sinw-2sin2w+2sin3w)^2] P4?(w)=41?[(1?cosw+2cos2w?2cos3w)2+(sinw?2sin2w+2sin3w)2]

自相關估計法

由 于 x 4 ( n ) 為 實 序 列 ， 所 以 x ? ( ? n ) = x ( ? n ) 由于x_4(n)為實序列，所以x^*(-n)=x(-n) 由于x4?(n)為實序列，所以x?(?n)=x(?n)

r 4 ( m ) = x ( m ) ? x ( ? m ) = 1 4 ∑ n = 0 3 x 4 ( n ) x 4 ( n ? m ) r_4(m)=x(m)*x(-m)=\dfrac{1}{4}\sum\limits_{n=0}^{3}x_4(n)x_4(n-m) r4?(m)=x(m)?x(?m)=41?n=0∑3?x4?(n)x4?(n?m)

r 4 ( m ) = 1 4 { ? 2 , 4 , ? 7 , 10 , ? 7 , 4 , ? 2 } n = ? 3 , ? 2 , . . . , 2 , 3 r_4(m)=\dfrac{1}{4}\{-2,4,-7,10,-7,4,-2\}\ \ \ \ n=-3,-2,...,2,3 r4?(m)=41?{?2,4,?7,10,?7,4,?2} n=?3,?2,...,2,3

P 4 ( w ) = ∑ ? 3 3 r 4 ( m ) e ? j w m P_4(w)=\sum\limits_{-3}^{3}r_4(m)e^{-jwm} P4?(w)=?3∑3?r4?(m)e?jwm

= 1 4 ? ( ? 2 e j 3 w + 4 e j 2 w ? 7 e j w + 10 ? 7 e ? j w + 4 e ? j 2 w ? 2 e ? j 3 w ) \ \ \ \ \ \ \ \ \ \ \ \ =\dfrac{1}{4}\cdot(-2e^{j3w}+4e^{j2w}-7e^{jw}+10-7e^{-jw}+4e^{-j2w}-2e^{-j3w}) =41??(?2ej3w+4ej2w?7ejw+10?7e?jw+4e?j2w?2e?j3w)

= 1 4 [ ? 2 ( e j 3 w + e ? j 3 w ) + 4 ( e j 2 w + e ? j 2 w ) ? 7 ( e j w + e ? j w ) + 10 ] \ \ \ \ \ \ \ \ \ \ \ \ =\dfrac{1}{4}[-2(e^{j3w}+e^{-j3w})+4(e^{j2w}+e^{-j2w})-7(e^{jw}+e^{-jw})+10] =41?[?2(ej3w+e?j3w)+4(ej2w+e?j2w)?7(ejw+e?jw)+10]

= 1 4 ( ? 4 c o s 3 w + 9 c o s 2 w ? 14 c o s w + 10 ) \ \ \ \ \ \ \ \ \ \ \ \ =\dfrac{1}{4}(-4cos3w+9cos2w-14cosw+10) =41?(?4cos3w+9cos2w?14cosw+10)

Matlab繪圖驗證二者功率譜影像：

Matlab代碼：

clc;
clear;
close all;
 
xn=[1,-1,2,-2];
N=length(xn);
rm=xcorr(xn); %計算自相關

figure('name','功率譜估計');
B=xn;
A=1;
[H1,w]=freqz(B,A,'whole');
subplot(2,2,1);
plot(w/pi,1/N*(abs(H1)).^2,'r');
xlabel('\omega/\pi');ylabel('P_4(\omega)');grid on;title('周期圖估計');

B=rm/N;
A=1;
[H2,w]=freqz(B,A,'whole');
subplot(2,2,2);
plot(w/pi,abs(H2),'b');
xlabel('\omega/\pi');ylabel('P_4(\omega)');grid on;title('自相關估計');

P1=0.25*((1-cos(w)+2*cos(2*w)-2*cos(3*w)).^2+(-sin(w)+2*sin(2*w)-2*sin(3*w)).^2); %周期圖估計結果
P2=0.25*(-4*cos(3*w)+8*cos(2*w)-14*cos(w)+10); %自相關估計結果

subplot(2,2,3);plot(w/pi,P1,'r');
xlabel('\omega/\pi');ylabel('P_4(\omega)');grid on;title('周期圖估計');
subplot(2,2,4);plot(w/pi,P2,'b');
xlabel('\omega/\pi');ylabel('P_4(\omega)');grid on;title('自相關估計');