功率谱估计和频率估计

N = 256; n = 0:N-1;
v = randn(size(n)); x = filter(num,den,v);
p = 4; [a e] = aryule(x,p); [H1 w1] = freqz(sqrt(e),a); [H2 w2] = freqz(num,den);
plot(w1/pi,abs(H1).^2,'b-'); hold on;grid on; plot(w2/pi,abs(H2).^2,'r--'); legend('估计功率谱','理论功率谱'); xlabel('Frequency: (\omega/\pi)');ylabel('Amplitude');
= ω1 0= .4π ,ω2 0.4= 5π ,ω3 0.8π ,
a. 假定已知 x(n) 中包含 3 个复谐波，用 Pisarenko 谐波分解来估计频率，并分析估计的
精度。分别重复 20 次实现，平均后的估计精度提高了吗？估计的方差降低了吗？如果过高 地估计频率个数，会出现什么情况？如果过低估计频率个数，会出现什么情况?
w(n)
是方差为
σ
2 w
的白高斯噪声，
x(n)
是
AR()
过程，由单位方差的白噪声通过如下滤波
器所获得
H
(
z)
=
1
−
1.585
z
1
−1
+
0.96
z
−2
a. 画出 x(n) 和 y(n) 的理论功率谱。
b.
取
σ
2 w
= 0.5,1, 2,5 ，取 y(n) 的 N
= 100 个样本，采用
p
=
2的
MEM
plot(w3/pi,abs(H3).^2); H4 = H4 + H3; hold on; end
H4 = H4/cnt; hold on;grid on; plot(w3/pi,abs(H4).^2,'r.'); xlabel('Frequency: (\omega/\pi)');ylabel('Amplitude'); title('采用自相关方法');
plot(w7/pi,abs(H7).^2); H8 = H7 + H8; hold on; end
H8 = H8/cnt; hold on;grid on; plot(w7/pi,abs(H8).^2,'r.'); xlabel('Frequency: (\omega/\pi)');ylabel('Amplitude'); title('采用协方差方法');
b. 编写子函数来估计复谐波过程的功率，并用该函数来估计 a 中各频率估计的功率。 用真实频率来估计功率，又会出现什么结果。
c. 分别用 MUSIC 方法、特征向量法，最小范数法来重复 a 中的估计 20 次，比较不同方 法间的估计精度。
close all;clear;clc; %%%%%%%%%%%%% a %%%%%%%%%%%%% sos = [1 0 0 1 -0.5 0.5;1 0 0 1 0 0.5]; [num den] = sos2tf(sos)
%%%%%%%%%%%%% e %%%%%%%%%%%%% sos = [1 0 0 1 -1.585 0.96;1 0 0 1 -1.152 0.96]; [num den] = sos2tf(sos)
v = randn(size(n)); x = filter(num,den,v);
p = 4; [aa ee] = aryule(x,p); [H1 w1] = freqz(sqrt(ee),aa); [H2 w2] = freqz(num,den);
%%%%%%%%%%%%% d %%%%%%%%%%%%% %协方差方法 figure;p=4; H8 = 0;cnt = 20; for i = 1:cnt
v = randn(size(n)); x = filter(num,den,v);
[a e] = arcov(x,p); [H7 w7] = freqz(sqrt(e),a);
MEM
功率谱，重复
c
中的过程。会提高功
-1-
率谱估计精度吗？ 试验三、本试验主要验证频率估计。
3
∑ = x(n)
令 x(n) 是谐波过程
i =1
Aie jnωi
+
w(n)
，其中
w(n) 是单位方差的高斯白噪声。令
= A1 4= e jφ1 , A2 3= e jφ2 , A3 e jφ3 ， φi 是 在 [π , −π ] 间 均 匀 分 布 的 不 相 关 随 机 变 量 ， 取
功率谱并与真实功率谱相比。 b. 重复 a 中的计算 20 次，分别画出 20 次的重迭结果和平均结果。评论估计的方差并
说明怎样才能提高自相关方法估计功率谱的精度；
c. 分别取 p = 6,8,12 来重复 b 中的计算，描述模型阶数增加时会出现什么结果。
d. 分别采用协方差方法、修改的协方差方法来重复 b,c 中计算过程，说明对宽带 AR 过 程而言，哪种方法最好。
H4 = H4/cnt; hold on;grid on; plot(w3/pi,abs(H4).^2,'r.'); xlabel('Frequency: (\omega/\pi)');ylabel('Amplitude'); title('采用自相关方法');
%协方差方法 figure;p=4; H8 = 0;cnt = 20; for i = 1:cnt
2*pi*rand(1)))... + A3*exp(j*(n*w3 + 2*pi*rand(1))) + noise;
R = covar(x,M); [v,d] = eig(R);
-7-
Pmu = 0; for i = 1:(M-p)
Pmu = Pmu + abs(fft(v(:,M-i+1),Nfft)).^2; end Var = mean(real(diag(d(p+1:M,p+1:M)))); Pmu = 1./Pmu;
%自相关方法 figure; H4 = 0;cnt = 20; for i = 1:cnt
v = randn(size(n)); x = filter(num,den,v);
[a e] = aryule(x,p); [H3 w3] = freqz(sqrt(e),a);
-5-
plot(w3/pi,abs(H3).^2); H4 = H4 + H3; hold on; end
figure plot(w1/pi,abs(H1).^2,'b-'); hold on;grid on; plot(w2/pi,abs(H2).^2,'r--'); legend('估计功率谱','理论功率谱'); xlabel('Frequency: (\omega/\pi)');ylabel('Amplitude: (dB)');
%重复 c 的计算省略 %修改的协方差方法 figure; H10 = 0;cnt = 20; for i = 1:cnt
v = randn(size(n)); x = filter(num,den,v);
-6-
[a e] = armcov(x,p); [H9 w9] = freqz(sqrt(e),a);
实验四 功率谱估计
实验内容、步骤：
实验内容包括三个：
实验一、宽带 AR 过程 x(n) 是由单位方差的高斯白噪声通过滤波器
H
(z)
=
(1 −
0.5 z −1
+
1 0.5 z −2
)(1 +
0.5 z −2
)
a. 生成 x(n) 的 N = 256 个样本，取 p = 4 并用自相关方法来计算功率谱，画出估计的
plot(w9/pi,abs(H9).^2);
-4-
H10 = H9 + H10; hold on; end
H10 = H10/cnt; hold on;grid on; plot(w9/pi,abs(H10).^2,'r.'); xlabel('Frequency: (\omega/\pi)');ylabel('Amplitude'); title('采用修改的协方差方法');
plot(w5/pi,abs(H5).^2); H6 = H5 + H6; hold on; end
H6 = H6/cnt;
-3-
hold on;grid on; plot(w5/pi,abs(H6).^2,'r.'); xlabel('Frequency: (\omega/\pi)');ylabel('Amplitude'); end
%%%%%%%%%%%%% c %%%%%%%%%%%%% for p = [6 8 12]
figure; H6 = 0;cnt = 20; for i = 1:cnt
v = randn(size(n)); x = filter(num,den,v);
[a e] = aryule(x,p); [H5 w5] = freqz(sqrt(e),a);
plot(w9/pi,abs(H9).^2); H10 = H9 + H10; hold on; end