采用iirgrpdelay函数设计全通滤波器校正系统相位

函数目标：根据描述的群延时设计全通滤波器对指定区间的相位进行校正。
iirgrpdelay函数共有以下7种用法：
1、[num,den] = iirgrpdelay(n,f,edges,a)
2、[num,den] = iirgrpdelay(n,f,edges,a,w)
3、[num,den] = iirgrpdelay(n,f,edges,a,w,radius)
4、[num,den] = iirgrpdelay(n,f,edges,a,w,radius,p)
5、[num,den] = iirgrpdelay(n,f,edges,a,w,radius,p,dens)
6、[num,den] = iirgrpdelay(n,f,edges,a,w,radius,p,dens,initden)
7、[num,den] = iirgrpdelay(n,f,edges,a,w,radius,p,dens,initden,tau)
拿第2个[num,den] = iirgrpdelay(n,f,edges,a) 来说：
输出：num是IIR滤波器分子，den是IIR滤波器分母
输入：
n是滤波器阶数，f是角频率范围矢量如0.3:0.001:0.6
edges是通带范围如[0.3,0.6]，a是描述的群延迟差如max(g)-g
下面是例子：

clc;clear all;close all;
load B_S;S=B_S;
load B_G;G=B_G;
H = dfilt.df2sos(B_S,B_G);
F=0.2:0.001:0.9;
g = grpdelay(H,F,2);   % Equalize the passband（g给0.2~0.9的群时延）
Gd = max(g)-g;          %Gd为需要补偿的群时延量
% Design the allpass delay equalizer
[num,den] = iirgrpdelay(20, F,[0.2 0.9], Gd);      %返回20阶全通滤波器，num和den是系数。
[GdA,w] = grpdelay(num,den);                    %计算全通滤波器的群时延
% 画图
He1=dfilt.df2(num,den);
He_all=dfilt.cascade(H,He1);
grpdelay(H)
grpdelay(He1)
grpdelay(He_all)  

实验结果：



可以看见校正后的系统相频响应接近一条直线，近似线性，校正效果非常好。