现代数字信号处理及应用仿真题答案

仿真作业

姓名：***

学号：S*********

4.17程序

clc;

clear;

for i=1:500

sigma_v1=0.27;

b(1)=-0.8458;

b(2)=0.9458;

a(1)=-(b(1)+b(2));

a(2)=b(1)*b(2);

datlen=500;

rand('state',sum(100*clock));

s=sqrt(sigma_v1)*randn(datlen,1);

x=filter(1,[1,a],s);

%%

sigma_v2=0.1;

u=x+sqrt(sigma_v2)*randn(datlen,1);

d=filter(1,[1,-b(1)],s);

%%

w0=[1;0];

w=w0;

M=length(w0);

N=length(u);

mu=0.005;

for n=M:N

ui=u(n:-1:n-M+1);

y(n)=w'*ui;

e(n)=d(n)-y(n);

w=w+mu.*conj(e(n)).*ui;

w1(n)=w(1);

w2(n)=w(2);

ee(:,i)=mean(e.^2,2);

end

end

ep=mean(ee');

plot（ep）;

xlabel('迭代次数');ylabel('MSE');title('学习曲线'); plot(w1);

hold;

plot(w2);

仿真结果：

步长0.015仿真结果

00.1

0.2

0.3

0.4

0.5

0.6

0.7

迭代次数M S E 学习曲线

步长0.025仿真结果

步长0.005仿真结果

4.18 程序

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;