IIR+全通滤波器实现相位平衡_matlab仿真

http://blog.csdn.net/shenziheng1/article/details/53409743

1.全通系统实现相位均衡

本文从逆向角度出发，基于搜索最优化策略对全通系统相位平衡问题进行了初步研究。主要思想：在频域中对它在单位圆上抽样N个点，按照所要求的频率特性，找出一个全通滤波器，使其满足一定的相位特性。最后利用最优化（例如遗传算法）策略进行迭代，找到满足要求的最优解。

适用于低通系统相位补偿的二阶全通滤波系统如下：

2.设计满足特定相位条件的全通滤波器

2.1 产生实验信号，组份0.01 0.02（保留） 0.6、0.7（低通滤波滤掉）

2.2 设计低通滤波器，进行滤波，并观察滤波后的信号是否发生了相位失真

低通滤波器作用下，高频成分被过滤掉，仅剩下低频成分0.01、0.02。

我们再看看滤波后，幅值和相位发生了什么变化？

我们发现，幅度图像有点变形！转过来看一下相位图像！！！问题出现了，发生了相位失真！就是说，经过低通滤波时0.01和0.02两个频率成分“相对步调”不一致，这就是群延时的意义，说明该滤波器的群延时不是常数。

3.构造特定全通滤波器相位，并利用遗传算法进行参数搜索设计

经过全通滤波器后，发现相位失真被补偿回来了！！！0.01/0.02步调又一致了。

最后，我们分析一下，相位补偿过程中，我们做的工作：

这张图才是相位均衡的核心，黑丝线条是低通巴特沃兹滤波器引进的非线性相位；蓝色曲线是我们利用全通滤波器构造的非线性相位，两者通过系统级联完成了线性相位的设计！

 

 

function [ delt ] = myfit( X  )
%X(1)  a
%X(2)  b
global AP_Phase
[h, w] = freqz([X(2),X(1),1],[1,X(1),X(2)]);
AP_Phase_temp =unwrap(angle(h))*180/pi;
delt = sum(abs(AP_Phase-AP_Phase_temp));  
end

　　

 

clc;clear all;close all;  
%%  生成模拟模拟信号  
% Fs = ;   
% t = ;   
% %归一化频率设计0.01 0.02 0.06  0.07   
% y = ;    
% yreal = ;  
%% 设计一个低通滤波器*********************************************%%  
fp=15/500;%通带截止频率 0.03  
fs=40/500;%阻带截止频率 0.06  
wp=3;     %通带最大衰减   
ws=30;    %阻带最大衰减   
[n,fc]=buttord(fp,fs,wp,ws);   
[b1,a1]=butter(n,fc);    
[h1,w]=freqz(b1,a1);  
LP_phase = unwrap(angle(h1))*180/pi;  %低通滤波器的相位函数  
%% 设计全通滤波器 进行相位补偿 ************************************%%  
Ideal_Phase = zeros(512,1);  
Y0 = min(min(LP_phase));  
for i=1:1:512  
    Ideal_Phase(i,1) = (Y0*i)/512;  
end  
%全通滤波器需要补偿的相位  
AP_Phase = zeros(512,1);  
for j=1:1:512   
    AP_Phase(j,1) = Ideal_Phase(j,1) - LP_phase(j,1);  
end  
global AP_Phase
%% 利用遗传算法优化设计二阶全通滤波器 ******************************%%  
%滤波器原型为：H(z)=（b+a*z^(-1)+z^(-2)）/(1+a*z^(-1)+b*z^(-2));  
%目标函数：minimum -> delt = sum(abs(AP_Phase-AP_Phase_temp));  
%动态范围：a->[-100,100]  b->[-100,100]  
options = gaoptimset(); %遗传算法
options.Generations = 200000; %最大迭代数设为20000  
UP = [50,50]; LOW = [-50,-50];%变量上下限  
[X,FVAL] = ga(@myfit, 2 ,[], [],[],[],LOW,UP,[],options);  
[h3,w]=freqz([X(2),X(1),1],[1,X(1),X(2)]);  
AP_Phase_Lab = unwrap(angle(h3))*180/pi;  %低通滤波器的相位函数  
%% 相位校正后的效果 **********************************************%%  
figure;
subplot(3,1,1)
plot(LP_phase, 'k')
hold on;
plot(AP_Phase);
hold on;
plot(Ideal_Phase, 'r')
subplot(3,1,2)
plot(AP_Phase_Lab,'b')
subplot(3,1,3)
plot(LP_phase+AP_Phase_Lab);