# 現代數字訊號處理大型作業

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源程式如下：

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演算法

``````<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');``````

``````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模型法

``````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');``````

`````` 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演算法

``````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）；通道具有脈衝響應：

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,格型大於橫向結構。