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

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