2. 试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。滤波器设计参数为：Fp＝1.7KHz, Fr＝2.3KHz, Fs＝8KHz, Armin≥70dB。
解：
III. Matlab程序与结果
1）二带数字滤波器组的源代码：（源代码统一格式排版）
%二带数字滤波器组
clear all;
close all;
%滤波器设计参数
Fp=1700;
Fr=2300;
Fs=8000;
Omegap=tan(pi*Fp/Fs);
Omegar=tan(pi*Fr/Fs);
Omegac=sqrt(Omegap*Omegar);
k=Omegap/Omegar;
k1=sqrt(sqrt(1-k^2));
q0=0.5*(1-k1)/(1+k1);
q=q0+2*q0^5+15*q0^9+150*q0^13;
N=11;
N2=fix(N/4);
M=fix(N/2);
N1=M-N2;
%计算Omgai和Vi
for i=1:M
a=0;
for m=0:5
a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*i/N);
end
b=0;
for m=1:5
b=b+(-1)^m*q^(m^2)*cos(2*m*pi*i/N);
end
Omega(i)=2*q^0.25*a/(1+2*b);
V(i)=sqrt((1-k*Omega(i)^2)*(1-Omega(i)^2/k));
end
%计算alhpa(i)和beta(i)
for i=1:N1
alpha(i)=2*V(2*i-1)/(1+Omega(2*i-1)^2);
end
for i=1:N2
beta(i)=2*V(2*i)/(1+Omega(2*i)^2);
end
%计算A(i)和B(i)
for i=1:N1
A(i)=(1-alpha(i)*Omegac+Omegac^2)/(1+alpha(i)*Omegac+Omegac^2);
end
for i=1:N2
B(i)=(1-beta(i)*Omegac+Omegac^2)/(1+beta(i)*Omegac+Omegac^2);
end
w=0:0.0001:0.5;
HL=zeros(size(w));
HH=zeros(size(w));
for n=1:length(w)
z=exp(j*w(n)*2*pi);
H1=1;
for i=1:N1
H1=H1*(A(i)+z^(-2))/(1+A(i)*z^(-2)) ;
end
H2=1/z;
for i=1:N2
H2=H2*(B(i)+z^(-2))/(1+B(i)*z^(-2));
end
HL(n)=abs((H1+H2)/2);
HH(n)=abs((H1-H2)/2);
end
plot(w,HL,'g',w,HH,'m');
hold on;
xlabel('digital frequency');
ylabel('Amptitude');
title('二带数字滤波器组');
2）运行结果
3）四带数字滤波器组的源代码：
%四带滤波器组
clear all;
close all;
%滤波器器设计参数
Fp=1700;
Fr=2300;
Fs=8000;
Omegap=tan(pi*Fp/Fs);
Omegar=tan(pi*Fr/Fs);
Omegac=sqrt(Omegap*Omegar);
k=Omegap/Omegar;
k1=sqrt(sqrt(1-k^2));
q0=0.5*(1-k1)/(1+k1);
q=q0+2*q0^5+15*q0^9+150*q0^13;
N=11;
N2=fix(N/4);
M=fix(N/2);
N1=M-N2;
%计算Omega(i)和V(i)
for i=1:M
a=0;
for m=0:5
a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*i/N);
end
b=0;
for m=1:5
b=b+(-1)^m*q^(m^2)*cos(2*m*pi*i/N);
end
Omega(i)=2*q^0.25*a/(1+2*b);
V(i)=sqrt((1-k*Omega(i)^2)*(1-Omega(i)^2/k));
end
%计算alpha(i)和beta(i)
for i=1:N1
alpha(i)=2*V(2*i-1)/(1+Omega(2*i-1)^2);
end
for i=1:N2
beta(i)=2*V(2*i)/(1+Omega(2*i)^2);
end
%计算A(i)和B(i)
for i=1:N1
A(i)=(1-alpha(i)*Omegac+Omegac^2)/(1+alpha(i)*Omegac+Omegac^2);
end
for i=1:N2
B(i)=(1-beta(i)*Omegac+Omegac^2)/(1+beta(i)*Omegac+Omegac^2);
end
w=0:0.0001:0.5;
HLL=zeros(size(w));
HLH=zeros(size(w));
HHL=zeros(size(w));
HHHP=zeros(size(w));
for n=1:length(w)
z=exp(j*w(n)*2*pi);
H1=1;
for i=1:N1
H1=H1*(A(i)+z^(-2))/(1+A(i)*z^(-2)) ;
end
H21=1;
for i=1:N1
H21=H21*(A(i)+z^(-4))/(1+A(i)*z^(-4)) ;
end
H2=1/z;
for i=1:N2
H2=H2*(B(i)+z^(-2))/(1+B(i)*z^(-2));
end
H22=1/(z^2);
for i=1:N2
H22=H22*(B(i)+z^(-4))/(1+B(i)*z^(-4));
end
HL=((H1+H2)/2);
HH=((H1-H2)/2);
HLL(n)=abs((H21+H22)/2*HL);
HLH(n)=abs((H21-H22)/2*HL);
HHH(n)=abs((H21+H22)/2*HH);
HHL(n)=abs((H21-H22)/2*HH);
end
plot(w,HLL,'b',w,HLH,'y',w,HHL,'r',w,HHH,'k')
hold on
xlabel('digital frequency');
ylabel('Amptitude');
title('四带数字滤波器组');
4）运行结果
3. 根据《现代数字信号处理》（姚天任等，华中理工大学出版社，2001）第四章附录提供的数据(pp.352-353)，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：
1)Levinson算法 2）Burg算法 3）ARMA模型法 4）MUSIC算法
1）Levinson算法的MATLAB源程序如下：
其中，参数X为输入序列，p为AR模型的阶数，函数调用形式为：Levinson（X，p）。这是一个Levinson.m文件:
Levinson.m文件
function S=Levinson(X,p)
N=length(X);
for m=0:N-1
R(m+1)=sum(conj(X([1:N-m])).*X([m+1:N]))/N;
end
a=-R(2)/R(1);
sigma2=(1-abs(a)^2)*R(1);
gamma=-a;
for k=2:p
sigma2(k)=R(1)+sum(a.*conj(R([2:k])));
D=R(k+1)+sum(a.*R([k:-1:2]));
gamma(k)=D/sigma2(k);
sigma2(k)=(1-abs(gamma(k))^2)*sigma2(k);
a=[a-gamma(k)*conj(fliplr(a)),-gamma(k)];
end
sigma2=real(sigma2);
f=linspace(-0.5,0.5,512);
for k=1:512
S(k)=10*log10(sigma2(end)/(abs(1+sum(a.*exp(-j*2*pi*f(k)*[1:p])))^2));
end
plot(f,S);
title(['Levinson: ',int2str(p),' 阶']);
xlabel('归一化频率')；
ylabel('相对谱/dB');
close all; clear all; p=[10 25 40 55]; x=randn(1,1000);%独立的高斯分布随机数 for k=1:4 subplot(2,2,k); Levinson(x,p(k)); end
3）运行结果
略
II） Burg算法
函数调用形式为：Burg（X，p）。其中参数X为输入序列，p为AR模型的阶数。
<pre name="code" class="plain">function [S,A,sigma2]=Burg(X,p)
N=length(X);
ef=X;
eb=X;
sigma2=sum(X*X')/N;
A=[];
for k=1:p
efp=ef(k+1:end);
ebp=eb(k:end-1);
gamma(k)=2*efp*ebp'/(efp*efp'+ebp*ebp');
sigma2(k+1)=(1-abs(gamma(k))^2)*sigma2(k);
ef(k+1:end)=efp-gamma(k)*ebp;
eb(k+1:end)=ebp-gamma(k)'*efp;
A=[A-gamma(k)*conj(fliplr(A)),-gamma(k)];
end
f=linspace(-0.5,0.5,512);
for k=1:512
S(k)=10*log10(sigma2(end)/(abs(1+sum(A.*exp(-j*2*pi*f(k)*[1:p])))^2));
end
plot(f,S);
title(['Burg: ',int2str(p),' 阶']);
xlabel('归一化频率'), ylabel('相对谱/dB');
分别对应于10阶、25阶、40阶、55阶的Burg算法的执行可以调用Burg函数实现。.m文件如下：
close all; clear all; p=[10 25 40 55]; x=randn(1,1000);%独立的高斯分布随机数 for k=1:4 subplot(2,2,k); Burg(x,p(k)); End
3）运行结果
略
III） ARMA模型法
函数调用形式为：ARMA（X，p，q）。其中参数X为输入序列，p为AR模型的阶数，q为MA模型的阶数。
function S=ARMA(X,p,q)
N=length(X);
M=N-1;
r=xcorr(X,'unbiased');
for k=1:p
R(:,k)=(r([N+q-k+1:N+M-k])).';
end
a=(-pinv(R)*(r([N+q+1:N+M]).')).';
Y=filter([1,a],[1],X);
pp=5*q;
[S,A,K]=Burg(Y,pp);
P=K(end);
[S,A,K]=Burg(A,q);
b=A;
f=linspace(-0.5,0.5,512);
for k=1:512
S(k)=10*log10(P*(abs(1+b*(exp(-j*2*pi*f(k)*[1:q]).'))/abs(1+a*(exp(-j*2*pi*f(k)*[1:p]).')))^2);
end
plot(f,S);
title(['ARMA: (',int2str(p),',',int2str(q),') 阶']);
xlabel('归一化频率');
ylabel('相对谱/dB');
分别对应于ARMA（30，10）、ARMA（30，1）、ARMA（40，10）、ARMA（50，10）的ARMA算法的执行可以用.m文件实现
close all; clear all; x=randn(1,1000);%独立的高斯分布随机数 subplot(2,2,1);ARMA(x,30,10); subplot(2,2,2);ARMA(x,30,1); subplot(2,2,3);ARMA(x,40,10); subplot(2,2,4);ARMA(x,50,10);
3)运行结果
略
IV） MUSIC算法
函数调用形式为：MUSIC（X，p）。其中参数X为输入序列，p为AR模型的阶数，
function S=MUSIC(X,p)
N=length(X);
r=xcorr(X,'biased');
clear R
for k=1:N
R(:,k)=(r([N-(k-1):2*N-k])).';
end
[V,D] = eig(R);
f=linspace(-0.5,0.5,512);
S0=fft(V(:,p+1:end),512);
for k=1:512
S(k)=10*log10(1/(S0(k,:)*S0(k,:)'));
end
S=[S(257:512) S(1:256)];
plot(f,S);
title(['MUSIC: ',int2str(p),' 维']);
xlabel('归一化频率');
ylabel('相对谱/dB');
分别对应于10维、25维、40维、55维的MUSIC算法的执行可以用.m文件实现：
clear all; close all; x=randn(1,1000);%独立的高斯分布随机数 p=[10 25 40 55]; for k=1:4 subplot(2,2,k); MUSIC(x,p(k)); end
4.图1为均衡带限信号所引起失真的横向或格型自适应均衡器（其中横向FIR系统长度M＝11），系统输入是取值为±1的随机序列x（n），其均值为零；参考d（n）＝x（n－7）；信道具有脉冲响应：
式中W用来控制信道的幅度失真（W＝2～4，例如，取W＝2.9,3.1,3.3,3.5），而且信道受到均值为零、方差为 （例如，取 ＝0.001，相当于信噪比为30dB）的高斯白噪声v（n）的干扰，试比较基于下列五种算法自适应均衡器在不同信道失真、不同噪声干扰下的收敛情况（对应于每一种情况，在同一坐标下画出其学习曲线），并分析其结果。
1>横向/格-梯型结构LMS算法 2>横向/格-梯型结构RLS算法
I）横向结构LMS算法
close all;
clear all;
w=[2.9,3.1,3.3,3.5];%信道的幅度失真
M=11;%横向滤波器的长度
sigma2=0.001; %噪声方差
T=7; %参考信号延时
N=400; % 训练次数（信号样本数）
iteration=500; %迭代次数
L=iteration+T+M-1; %单个输入信号长度
u=0.025; %迭代步长
value=zeros(length(w),L-M+1-T);
for ww=1:length(w)
%信道具有的脉冲响应
h=ones(1,3);
h(1)=0.5*(1+cos(2*pi/w(ww)));
h(3)=h(1);
e2=zeros(N,L-M+1-T);%误差信号的平方
for n=1:N
rand('seed',n*N);
X=sign(2*rand(1,L)-1); %产生一个随机信号序列
D=zeros(size(X));
for mm=T+1:L
D(mm)=X(mm-T);%对应的参考信号
end
U=conv(X,h);%信号经过信道的输出
randn('seed',n*N);
V=randn(size(U)).*sqrt(sigma2); %产生高斯噪声
R=U+V; %自适应滤波器输入信号
W=zeros(M,1); %滤波器参数的初始值为0
for m=1:iteration
r=R(T+m:T+m+M-1)';
y=r'*W;
e=D(T+m+M-1)-y; % 误差信号
e2(n,m)=e.^2;
W=W+2*u*e*r; %滤波器参数迭代
end
end
value(ww,:)=mean(e2); %均方误差值
end
semilogy(value(1,:),'y--');hold on;
semilogy(value(2,:),'k-.');hold on;
semilogy(value(3,:),'m:');hold on;
semilogy(value(4,:),'b-');hold on;
grid on;
xlabel('the number of iteration');
ylabel('mean square error');
title('LMS Arithmetic');
legend(‘w1=2.9’,‘w2=3.1’,‘w3=3.3’,‘w4=3.5’);
II）横向结构RLS算法
close all;
clear all;
w=[2.9,3.1,3.3,3.5];
M=11;
sigma2=0.001; %噪声方差
T=7;%参考信号延迟
N=400;%训练次数
iteration=500; %迭代次数
L=iteration+T+M-1;%输入信号长度
lamda=0.999; %遗忘因子
value=zeros(length(w),L-M+1-T);
for ww=1:length(w)
h=ones(1,3);
h(1)=0.5*(1+cos(2*pi/w(ww)));
h(3)=h(1);
e2=zeros(N,L-M+1-T);
for n=1:N
rand('seed',n*N);
X=sign(2*rand(1,L)-1);
D=zeros(size(X));
for mm=T+1:L
D(mm)=X(mm-T);
end
U=conv(X,h);
randn('seed',n*N);
V=randn(size(U)).*sqrt(sigma2);
R=U+V; %自适应滤波器输入信号
W=zeros(M,1);
P=eye(M)/0.004;
for m=1:iteration
r=R(T+m+M-1:-1:T+m)';
e=D(T+m+M-1)-r'*W;
e2(n,m)=e.^2;
K=P*r/(lamda+r'*P*r);
P=(eye(M)-K*r')*P/lamda;
W=W+K.*e;
for mm=2:M-1
r(mm)=r(mm+1);
end
r(1)=R(T+m+M);
end
end
value(ww,:)=mean(e2);
end
semilogy(value(1,:),'y--');hold on;
semilogy(value(2,:),'k-.');hold on;
semilogy(value(3,:),'m:');hold on;
semilogy(value(4,:),'b-');hold on;
grid on;
xlabel('the number of iteration');
ylabel('mean square error');
title('RLS Arithmetic');
III）格-梯型结构LMS算法
close all;
clear all;
W=[2.9,3.1,3.3,3.5];%信道的幅度失真
w=0.9999;
T=7;%延迟7个周期
sigma2=0.001;%噪声方差
L=500;%数据长度
M=11;%横向滤波器的长度
N=400;%次数
e2=zeros(length(W),L);
for ww=1:length(W)
n=[1,2,3];
h(n)=1/2*(1+cos(2*pi*(n-2)/W(ww)));
for nn=1:N
km=ones(L+2,M+1);
vm=ones(L+2,M+1);
beta=ones(L+2,M+1);
fm=ones(L+2,M+1);
gm=ones(L+2,M+1);
e=ones(L+2,M+2);
rand('seed',N*nn);
X=sign(2*rand(1,L+T)-1);
U=conv(X,h);
randn('seed',N*nn);
V=randn(size(U))*sqrt(sigma2);
R=U+V;
for m=1:M+1
km(2,m)=0; %n=2
vm(2,m)=0.8;
beta(3,m)=0;
fm(2,m)=0;
gm(2,m)=0;
gm(1,m)=0;
end
for n=3:L+2 %m=1
e(n,2)=X(n-2);
fm(n,1)=R(n-2+T);
gm(n,1)=R(n-2+T);
end
for n=3:L+2
for m=2:M+1
vm(n,m)=w*vm(n-1,m)+(abs(fm(n,m-1)))^2+(abs(gm(n-1,m-1)))^2; km(n,m)=km(n-1,m)+(fm(n-1,m-1)*conj(gm(n-1,m))+conj(gm(n-2,m-1))*fm(n-1,m))/vm(n-1,m);
fm(n,m)=fm(n,m-1)-km(n,m)*gm(n-1,m-1);
gm(n,m)=gm(n-1,m-1)-conj(km(n,m))*fm(n,m-1);
e(n,m+1)=e(n,m)-gm(n,m)*beta(n,m);
beta(n+1,m)=beta(n,m)+2*e(n,m+1)*conj(gm(n,m))/vm(n,m);
end
e2(ww,n-2)=e2(ww,n-2)+e(n,M+2)*e(n,M+2);
end
end
end
e2=e2/N;
semilogy(1:L,e2(1,:),'y-');hold on;
semilogy(1:L,e2(2,:),'b-.');hold on;
semilogy(1:L,e2(3,:),'k:');hold on;
semilogy(1:L,e2(4,:),'m--');hold on;
grid on;
xlabel('data');
ylabel('mean square error');
title('格-梯型结构LMS算法');
IV）格-梯型结构RLS算法
close all;
clear all;W=[2.9,3.1,3.3,3.5];%信道的幅度失真
lamda=0.999; %遗忘因子
T=7;%延迟7个周期
L=500; %数据长度
M=11;%横向滤波器的长度
N=400; %迭代次数
sigma2=0.001;%噪声方差
e2=zeros(length(W),L);
for i=1:length(W)
h=ones(1,3);
h(1)=0.5*(1+cos(2*pi/W(i)));
h(3)=h(1);
alpha=zeros(L+2,M);
k=zeros(L+2,M);
kf=zeros(L+2,M);
kb=zeros(L+2,M);
fm=zeros(L+2,M);
gm=zeros(L+2,M);
Ef=zeros(L+2,M);
Eb=zeros(L+2,M);
delta=zeros(L+2,M);
kasi=zeros(L+2,M);
e=zeros(L+2,M+1);
for j=1:N
rand('seed',N*j);
X=sign(2*rand(1,L+T)-1);
D=zeros(size(X));
U=conv(X,h);
randn('seed',N*j);
V=randn(size(U)).*sqrt(sigma2);
R=U+V; %自适应滤波器输入信号
for m=1:M
alpha(1,m)=1;
k(1,m)=0;
Eb(1,m)=0.9;
Eb(2,m)=0.9;
Ef(2,m)=0.9;
delta(1,m)=0;
fm(2,m)=0;
gm(2,m)=0;
gm(1,m)=0;
e(2,m)=0;
end
for n=3:L+2
alpha(n-1,1)=1;
e(n,1)=X(n-2);
fm(n,1)=R(n-2+T);
gm(n,1)=R(n-2+T);
Ef(n,1)=lamda*Ef(n-1,1)+(abs(R(n-2+T)))^2;
Eb(n,1)=Ef(n,1);
end
for n=3:L+2
for m=1:M-1
k(n-1,m+1)=lamda*k(n-2,m+1)+alpha(n-2,m)*fm(n-1,m)*conj(gm(n-2,m));
kf(n-1,m+1)=-k(n-1,m+1)/Eb(n-2,m);
kb(n-1,m+1)=-conj(k(n-1,m+1))/Ef(n-1,m);
fm(n,m+1)=fm(n,m)+kf(n-1,m+1)*gm(n-1,m);
gm(n,m+1)=gm(n-1,m)+kb(n-1,m+1)*fm(n,m);
Ef(n-1,m+1)=Ef(n-1,m)-(abs(k(n-1,m+1)))^2/Eb(n-2,m);
Eb(n-1,m+1)=Eb(n-2,m)-(abs(k(n-1,m+1)))^2/Ef(n-1,m);
alpha(n-1,m+1)=alpha(n-1,m)-(alpha(n-1,m))^2* (abs(gm(n-1,m)))^2/Eb(n-1,m);
end
for m=1:M
delta(n-1,m)=lamda*delta(n-2,m)+alpha(n-1,m)* conj(gm(n-1,m))*e(n-1,m);
kasi(n-1,m)=-delta(n-1,m)/Eb(n-1,m);
e(n,m+1)=e(n,m)+kasi(n-1,m)*gm(n,m);
end
e2(i,n-2)=e2(i,n-2)+e(n,M+1)^2;
end
end
end
e2=e2/N;
clf;
semilogy(1:L,e2(1,:),'y-');hold on;
semilogy(1:L,e2(2,:),'b-.');hold on;
semilogy(1:L,e2(3,:),'c:');hold on;
semilogy(1:L,e2(4,:),'k--');hold on;
grid on;
xlabel('data');
ylabel('mean square error');
title('格-梯型结构RLS算法');
几种算法结果比较分析：
(1)在不同的信道幅度失真下，信道幅度失真越小，收敛时的均方误差越小。横向LMS算法对信道幅度失真十分敏感，当w=3.5时，横向LMS不仅稳态误差较大，其收敛速度也明显变慢。
(2)从收敛速度上看，横向RLS和格型RLS快于横向LMS和格型LMS。
(3)从稳态误差上看，横向RLS和格型RLS普遍低于横向LMS和格型LMS。
(4)从计算量上看,RLS大于LMS,格型大于横向结构。