allpass5.m

计算全通滤波器的源程序,根据信道群延迟,计算滤波器系数,用以校正信号相位变化

clc;clear;close all;
f=10.23e+6;
fs=60.138e+6;
%%%%%%%%%%用最小二乘拟合相位曲线%%%%%%%%%%%%%%%%%
% ssamples1=[101	159	264	215	300	185	409	326	122	335	214];
% ssamples2=[43	101	203	155	241	126	349	267	61	276	152];
% ssamples3=[-15	40	147	99	173	70	284	213	-5	224	87];
% ssamples4=[-5	53	158	108	190	77	303	221	16	227	108];
% ssamples5=[	-97	-41	67	17	99	-11	204	131	-86	146	4];
% ssamples6=[	44	104	208	157	247	129	354	268	68	277	160];
% ssamples7=[	153	211	319	267	350	239	461	378	174	387	264];
% ssamples8=[	14	71	173	125	215	94	328	234	41	243	132];
% bfreq=[	0.18	0.22	0.38	0.42	0.42	0.58	0.62	0.78	0.82	0.98	1.023]*1e+6;
% dp=ssamples3.*pi./180;
% bf=bfreq;
% sw=2*pi.*bf./fs;
% %alf(1)=(sum(dp.*sw)-sum(sw.^2));
% [alf,S]= polyfit(sw,dp,1);

%%%%%%%%%  计算全通滤波器   %%%%%%%%%%%%%%%%
N=1024;
Nn=10;              %%%%%%%%%%  滤波器阶数  %%%%%%%%
M=9;
u=ones(N,N);
a=zeros(1,Nn+1);
a(1)=1;
n1=(0:1:N-1);
w=pi.*n1./N;
% alf=(1:1:50);
% beta=M-alf*pi;
% Theta=-(alf(1)*w.^2+beta(1)*w);
% Tau=2*alf(1)*w+beta(1);
% Theta=-(5*w.^2-w.^3);
% Theta=-(alf(1).*w.^2+alf(2).*w);
Theta=-((4/pi)*w.^2+6*w);
% Tau=10*w-3*w.^2;
Tau=(8/pi)*w+6;                      %给出期望群延迟
TauD=(-Tau./2+Nn/2);                 %分母群迟延
eTau=sum(Tau)/N;                     %群迟延平均值
Taus=(Tau-eTau)/2; %
% N2=round(eTau+Nn)+1;
% ea=[zeros(1,N2-1) 1];
%%%%%%%%%%  最小均方误差的求复倒谱方法  %%%%%%%%%
for n2=1:1:N
    for n3=0:1:N-1
     u(n3+1,n2)=n2*cos(n2*(pi*n3/N));
    end
end
W=[1.*ones(1,Nn)  1000.*ones(1,N-2*Nn) 1.*ones(1,Nn)];
bw=(W.*TauD)';

 A1=diag(W)*u;
 
d=(A1'*A1)\(A1'*bw);
d=d';
%%%%%%%%%     FFT   %%%%%%%%%%
%d=2.*ifft(TaoDe);
%%%%%%%%%%%%%%%%%%%%%%%
d=[0 d];

%%%%%%%%   求H(z)分母系数    %%%%%%%%%%%%%%
for n=1:1:Nn
    for k=1:1:n
        a(n+1)=a(n+1)+k*d(k+1)*a(n+1-k)/n;
    end
end

 den=a;
 num=fliplr(den);
% num=conv(a,ea);                   %%%%H(z)分母系数%%%
%   den=fliplr(num);      %%%%H(z)分子系数%%%%%%
% div=den(1);
% den=den./div;
% num=num./div;
%%%%%%%%%画图%%%%%%%%%
figure(1);
plot(w./pi,Theta,'r.');grid;
hold on;
[h,w]=freqz(num,den,1024);
phz=unwrap(angle(h));
plot(w./pi,phz);grid;
% axis([0,1,-30,40]);
title('理论期望的相位和计算仿真的展开相位特性')
figure(2);
freqz(num,den);
title('计算仿真的幅频特性和相频特性');
figure(3);
plot(w./pi,Tau,'r.');grid;
hold on;
grpdelay(num,den);
title('期望的群延迟和计算仿真的群延迟');