功率谱估计和频率估计

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
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)');
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);
wk.baidu.com
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');
e. 把宽带 AR 过程改为下列窄带 AR 过程，
H
(z)
=
(1 − 1.585 z −1
+
0.96 z −2
1 )(1 − 1.152 z −1
+
0.96 z −2
)
重复 a,b,c,d 中的所有分析。
实验二、本实验是验证最大熵方法的功率谱估计。 对随机过程 y= (n) x(n) + w(n) ，
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;
p = 3; M = p + 3;
cnt = 20; vect = 0; vale = 0; Vsum = 0; Psum = 0;
for m = 1:cnt noise = randn(size(n)) + j*randn(size(n)); x = A1*exp(j*(n*w1 + 2*pi*rand(1))) + A2*exp(j*(n*w2 +
%%%%%%%%%%%%% 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);
-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('采用修改的协方差方法');
v = randn(size(n)); x = filter(num,den,v);
[a e] = arcov(x,p); [H7 w7] = freqz(sqrt(e),a);
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('采用协方差方法');
= ω1 0= .4π ,ω2 0.4= 5π ,ω3 0.8π ,
a. 假定已知 x(n) 中包含 3 个复谐波，用 Pisarenko 谐波分解来估计频率，并分析估计的
精度。分别重复 20 次实现，平均后的估计精度提高了吗？估计的方差降低了吗？如果过高 地估计频率个数，会出现什么情况？如果过低估计频率个数，会出现什么情况?
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('采用协方差方法');
%%%%%%%%%%%%% 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);
%重复 c 的计算省略 %修改的协方差方法 figure; H10 = 0;cnt = 20; for i = 1:cnt
v = randn(size(n)); x = filter(num,den,v);
[a e] = armcov(x,p); [H9 w9] = freqz(sqrt(e),a);
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)
close all;clear all;clc; %%%%%%%%%%%%% a %%%%%%%%%%%%% N = 500; n = 0:N-1; Nfft = 1024; w = linspace(0,2*pi,Nfft);
A1 = 5; A2 = 4; A3 = 3; w1 = 0.4*pi; w2 = 0.55*pi; w3 = 0.8*pi;
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 是 在 [π , −π ] 间 均 匀 分 布 的 不 相 关 随 机 变 量 ， 取
plot(w9/pi,abs(H9).^2); 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(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('采用自相关方法');
%%%%%%%%%%%%% 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);
方法由
y(n) 来
估计 x(n) 的功率谱，看看噪声对功率谱估计的精度有多大影响。
c. 改 p = 5 ，再重复 b 中的过程，分析所观测的结果；
d. 由于自相关序列为 ry= (k) rx (k) + σ w2δ (k) ，如果在计算 MEM 功率谱前从自相关值
ry
(0)
中减去
σ
2 ω
，用修改后的自相关序列来估计
实验四 功率谱估计
实验内容、步骤：
实验内容包括三个：
实验一、宽带 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 并用自相关方法来计算功率谱，画出估计的
%重复 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);
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
-2-
%%%%%%%%%%%%% b %%%%%%%%%%%%% %自相关方法 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);
功率谱并与真实功率谱相比。 b. 重复 a 中的计算 20 次，分别画出 20 次的重迭结果和平均结果。评论估计的方差并
说明怎样才能提高自相关方法估计功率谱的精度；
c. 分别取 p = 6,8,12 来重复 b 中的计算，描述模型阶数增加时会出现什么结果。
d. 分别采用协方差方法、修改的协方差方法来重复 b,c 中计算过程，说明对宽带 AR 过 程而言，哪种方法最好。
%自相关方法 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
w(n)
是方差为
σ
2 w
的白高斯噪声，
x(n)
是
AR(2)
过程，由单位方差的白噪声通过如下滤波
器所获得
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