对于 AR 模型法，由于其对于给定的序列，它隐含着数据和自相关函数的外推，使其可 能的数据长度超过给定的长度，因而分辨率可以不受反比于数据长度 N 的限制，得到很高的 分辨率。图中也可以看出，对于 Berg 法可以很清晰的分辨 0.21 和 0.22 这两个靠得很近的 频率。由于 AR 模型是一个有理分式，它估计出的谱要比周期图法平滑。
分辨率为 2/N ，足以把 0.15 和 0.16 的两条谱线分开。方差性能不是很好，所以谱线不是
很光滑。
2．用 Welch 方法对上述周期图法作改 进，例如，可选每段 64 点，重叠 32
点，用 Hamming 窗，画出所求功率谱曲线 。
5 0 -5 -10 Bd) -15 w P -20 -25 -30 -35
多少为最好，还要在实践中对所得的结果做多次比较后予以确定。
3．试比较用周期图法和模型法估计出的功率谱 的分辨率及其他性能。
周期图法物理概念明确，并且可以用 FFT 算法直接计算，计算量相对较小。但是周期图
法的分辨率较低，它正比于 2/ N （N 指所用数据的长度），这是受窗函数的影响，使得真
正谱在窗口主瓣内的功率向边瓣泄漏，降低了分辨率。周期图法的方差性能也不好，它不是 真实功率谱的一致估计。周期图的平滑和平均法主要是用来改善方差性能，但它们都和所用 的窗函数紧密相关，往往在方差得到改善的同时，分辨率和偏差却增大。这些在上面两题功 率谱估计的结果中可以得到证明：周期图法的分辨率相对要低于模型法的。
AR 模型 自相关法和 Berg 算法 %homework010301 clear all;clc; load x(n) x; N=4096; Order=33; w=-0.5:1/N:0.5-1/N; xpsd=abs(pyulear(x,order,N,1)); pmax=abs(max(xpsd)); xpsd=xpsd/pmax; xpsd=10*log10(xpsd+0.000001); figure(1) plot(w,fftshift(xpsd));grid on;
5
0
-5
-10
-15 Bd )w -20 P
-25
-30
-35
-40
AR 模型自相关法估计 x(n)的功率谱
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 w/2pi
5
0
-5
--4150
-15
Bd) w
-20
P -25
-30
-35
-40
AR 模型 Burg 法估计 x(n)的功率谱
Байду номын сангаас
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 w/2pi
2．试利用 FPE 或 AIC 准则来确定本信号所需要的 AR 模型的阶次。
用 FPE 和 AIC 准则计算出的 AR 模型都是 33 阶。采用估算出的阶次再做功率谱估计如下图所 示。图中出现了虚假的峰值，原因是阶次有点大了。
FPE 准则和 AIC 准则计算 AR 模型阶次 %homework 010302 clear all; clc; load x(n).mat; N=256; rx=xcorr(x,N,'biased'); R=rx(N+1:2*N+1); E=zeros(1,N); for order=1:N; [a,E(order),k]=levinson(R,order); end for i=1:length(x)
周期图法 %homework010201 clear;clc; load x(n) x; N=length(x); M=4096; w=-0.5:1/M:0.5-1/M; px = periodogram(x,boxcar(N),M,1); xpsd=10*log10(px/max(px)+0.000001); figure(1) plot(w,fftshift(xpsd));grid on; ylabel('P(w)/dB'); xlabel('w/(2pi)'); title('周期图法估计 x(n)的功率谱');
)see -600
ged esah
-800
P -1000
-1200
-1400
-160 00
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Normalized Frequency (x2pi rad/sample)
（二）经典功率谱估计
j
1．用周期图估计 X 'n的功率谱，画出 P per (e ) 曲 线；
Bd) -20 w P -30
-40
-50
-60 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 w/2pi
经调整，当阶次取 25 的时候能够达到比较好的估计质量（图形不再给出）。这里的 FPE
准则和 AIC 准则仅为阶次的选择提供了依据，对所研究的某一个具体的信号，究竟阶次取为
Welch 方法 %homework010202 clear;clc;
load x(n) x; M=4096; w=-0.5:1/M:0.5-1/M; px = pwelch(x,hamming(65),32,M,'twosided');%Welch,Hamming xpsd=10*log10(px/max(px)+0.000001); figure(1) plot(w,fftshift(xpsd));grid on; ylabel('P(w)/dB');xlabel('w/(pi)'); title('welch 估计 x(n)的功率谱');
welch 估 计 x(n)功 率 谱
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 w/(pi)
Welch 方法估计功率谱如下，由于分段并且有重叠，方差性能得到了很大的改善，谱线 变得平滑了，但是分辨率明显降低。可以看出，能够将-0.16 与 0.15、0.16 个频率分开，但 是没有将 0.15 和 0.16 区分开，这是由于分段使得数据长度变小所致。
附：代码如下
滤波器设计： %homework0101 clear all;clc; load test.mat x f=[0 0.4 0.5 1]; A=[1 1 0 0]; weigh=[10 1]; b=remez(32,f,A,weigh); [h,f]=freqz(b,1,256,1); h1=abs(h); h2=20*log10(h1); figure(1) plot(f,180/pi*unwrap(angle(h)));grid; xlabel('Normalized Frequency (x2pi rad/sample)'); ylabel('Phase(degrees)') figure(2) plot(f,h2);grid; xlabel('Normalized Frequency (x2pi rad/sample)'); ylabel('Magnitude(dB)') y=conv(b,x); x=[y,zeros(1,96)]; save x(n) x;
说明：在画功率谱曲线时，频率分点取 409 6，最大值按 0dB 归一化。 ^ 周 期 图 法 估 计 x(n)功 率 谱 10
0
-10
Bd)w -20 P -30
-40
-50
-60 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 w/(2pi)
利用 Matlab 中的 periodogram 函数做功率谱估 计，功率谱如上图所示。本题中取 N=4096 点 FFT 计算。(程序见附录)从图中可以看出，能够很明显的区分开 3 个不同频率的谱线，
^
3．用自相关法估计 X 'n的功率谱，令 M=63，画出 P per (e j ) 曲线；
自相关法估计功率谱如下图：-0.16 处频率已经很难与噪声区分开，并且 0.15 和 0.16 频 率也没有很好的区分开，可见，BT 法起到了平滑曲线的作用，但由于对数据进行截短，降低 了分辨率。
-40
自 相 关 法 估 计 x(n)功 率 谱 10
10
0
-10
-20
)
Bd
eduing
-30 -40
a
M
-5 0
X: 0.25 Y: -24.3
-6 0
-7 0
0
-8 0
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Normalized Frequency (x2pi rad/sample)
0
-200
-400
由下图，功率谱的谱线十分平滑，但分辨率都不好，没能区分开 0.15 和 0.16 两个频率成 分。因为自相关法等效于对前向预测的误差序列前后加窗，因而使得该方法分辨率降低。由 于所设计滤波器的阶数比较低，数据通过滤波器后的长度只有 160 点，在后面补零造成了分 辨率很低。调整模型的阶次也可以相应提高分辨率，当模型阶数在 25 左右就能把 0.15 和 0.16 分开，但是不符合题目要求。
0
-10
Bd) -20 wtb P -30
-40
-50
-60 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 w/(2pi)
（三）用模型法做功率谱估计
1．用 AR 模型估计 X'n的功率谱，阶次 p 可在 6－15 之间，由自己给定，给