现代数字信号处理及其应用作业(homework.m)

【核心代码】

%（1）产生N=512点的样本序列
data_len = 512;    %样本序列的长度
trials = 100;     %随机试验的次数
A=zeros(data_len,2);EA=zeros(data_len,1);
B=zeros(data_len,2);EB=zeros(data_len,1);
for m = 1: trials
a1 = -0.975;
a2 = 0.95;
sigma_v_2 =0.0731;
v = sqrt(sigma_v_2) * randn(data_len, 1, trials);%产生v(n)
u0 = [0 0];
num = 1;
den = [1 a1 a2];
Zi = filtic(num, den, u0);     %滤波器的初始条件
u = filter(num, den, v, Zi);  %产生样本序列u(n)
%（2）用LMS滤波器来估计w1和w2
mu1 = 0.05;
mu2 = 0.005;
w1 = zeros(2, data_len);
w2 = zeros(2, data_len);
e1 = zeros(data_len, 1);
e2 = zeros(data_len, 1);
d1 = zeros(data_len, 1);
d2 = zeros(data_len, 1);
%LMS迭代过程
for n =3 :data_len - 1
w1( :, n 1) = w1( :, n)   mu1 * u(n-1 : -1: n-2, : , m) * conj(e1(n));
w2( :, n 1) = w2( :, n)   mu2 * u(n-1 : -1: n-2, : , m) * conj(e2(n));
d1(n 1) = w1( : , n 1)' * u(n: -1: n-1, :, m);
d2(n 1) = w2( : , n 1)' * u(n: -1: n-1, :, m);
e1(n 1) = u(n 1, : ,m) - d1(n 1);
e2(n 1) = u(n 1, : ,m) - d2(n 1);
end
A = A   conj(w1)';
EA = EA  e1.^2;
B = B   conj(w2)';
EB = EB   e2.^2;
end
%剩余均方误差和失调参数
wopt=zeros(2,trials);
Jmin=zeros(1,trials);
sum_eig=zeros(trials,1);
for m=1:trials;
rm=xcorr(u(:,:,m),'biased');
R=[rm(512),rm(513);rm(511),rm(512)];
p=[rm(511);rm(510)];
wopt(:,m)=R\p;
[v,d]=eig(R);
Jmin(m)=rm(512)-p'*wopt(:,m);
sum_eig(m)=d(1,1) d(2,2);
end
sJmin=sum(Jmin)/trials;
e1_100trials_ave=sum(e1)/trials;
e2_100trials_ave=sum(e2)/trials;
Jex1=e1_100trials_ave-sJmin;
Jex2=e2_100trials_ave-sJmin;
sum_eig_100trials=sum(sum_eig)/100;
Jexfin=mu1*sJmin*(sum_eig_100trials/(2-mu1*sum_eig_100trials));
Jexfin2=mu2*sJmin*(sum_eig_100trials/(2-mu2*sum_eig_100trials));
M1=Jexfin/sJmin
M2=Jexfin2/sJmin
figure(1);
plot(A/trials);hold on;
plot(conj(w1)');
xlabel('迭代次数');ylabel('权向量');title('步长为0.05权向量收敛曲线');
figure(2);
plot(B/trials);hold on;
plot(conj(w2)');
xlabel('迭代次数');ylabel('权向量');title('步长为0.005权向量收敛曲线');
figure(3);
plot(EA/trials,'*');hold on;
plot(EB/trials,'-');
xlabel('迭代次数');ylabel('均方误差');title('步长分别为0.05和0.005学习曲线');