FIR等波纹滤波器优化设计.docx

FIR等波纹滤波器优化设计0,wp段位通带 ws,阻带 wp,ws过度带1是理想通带中能接受的振幅波动2是阻带中能接受的振幅波动通带波动dB值 Rp=-20log101-11+2阻带波动dB值 As=-20log1021+2 M阶数设计方法1. 窗函数法理想低通滤波器Hdejw=e-jn |w|wc0 wcwpi 对应时域hdn=12-wcwce-jn0wejnwdw=sin(wc(n-n0)(n-n0)无限长，非因果，无法实现hn=hd(n)RN(n)选取傅里叶级数的项数越多，引起的误差就越小另一种思路，不只简单地截取，对输出的傅里叶级数加权，就是用适当的窗函数形状进行加权2. 频域样本法先在频率域修改滤波器的预期频率特性，使它不太理想，接近于实际的需要和可能，然后直接利用它的傅里叶反变换作为滤波器系数。最优设计法 上面两种方法缺点：1不能精确地指定阻带频率ws和wp 2不能同时控制1和2 3 理想响应和实际响应之间的逼近误差在全频带区间不均匀，靠近边缘频率处误差较大，而在远离边缘处误差误差较小。（如果使误差均匀分布，就可以在满足相同技术指标的条件下，得到一个较低阶的滤波器。）A最小二乘法B最小最大波动法ew=WtwAw-D(w)=Wt(w)n=0Ld(n)cosnw-D(w)Wt(w)幅特性1和1.01 阻带0.011和0.001低通常取 Wtw=fx=12 0wwp1 wpw问题找到滤波器系数向量d(n)，使得在w=0,wp wp,ws频段内最大绝对误差ew为最小交替定理 Hejw=n=-MMhne-jnw H(ejw)=h(0)+n=1M2h(n)cosnw设w是在0,上的一个子集, 比如是0,wp和wp,ws的合集。则 Aw=n=0Ld(n)cosnw 是对D(w)的最佳逼近，当且仅当式的误差函数e（w）是等波动的，并且该误差函数e（w）在w中至少有L+2个极值点。即在w中，存在着 0w1w2若发现有新的极值点代替原有极值点，重复Step1，直到所有E(w)，且不变。代码：%主函数clear;wp=0.4*pi;ws=0.6*pi;jishuqi = 0; %用于测试M = 7; %滤波器设计指标,通带截止频率0.4pi,阻带截止频率0.6pi,阶数为7countwhile = 0; %迭代次数pold = ; %存储上一个计算的p值,用于测试迭代是否结束l=wp*(M+1)/(pi-ws+wp); %按照均分原则,确定通带内极值点频率的个数w1=wp/(l+1):wp/(l+1):wp; %通带内初始化极值点w2=ws:(pi-ws)/(M-l-1):pi; %阻带内初始化极值点w=w1,w2;%w = Initial_wset; %M+1 lengthwc = 0.5*pi; %理想低通滤波器截止频率%main loop%迭代开始while(1) jishuqi = jishuqi + 1; count = 0; H_MATRIX = ; E_MATRIX = ; a_k_Matrix,H_k_Matrix,H_d_Matrix,p = Remez_DSP(M,w,wc,wp); %根据滤波器阶数,极 %值点求出此时ak,Hk,Hd,p for index=1:length(w) %此步为了保留通带截止频率wp阻带截止频率ws不改变 if w(index) = wp wploc = index; end if w(index) = ws wcloc = index; endend %更新极值点 w_candidate = 0.0001:(w(wploc)-0.0001)/10000:w(wploc)-1/10000; %在通带内均分取10000个点,为了求出新的极值点 E_candidate = ; neww = ; for jj=1:length(w_candidate) E,H = Remez_H(w_candidate(jj),H_k_Matrix,wp,wc,M,a_k_Matrix,w); E_candidate = E_candidate E; End %求出10000个点所对应的E值 for jj=2:length(E_candidate)-1 if abs(E_candidate(jj)abs(E_candidate(jj-1) & abs(E_candidate(jj)abs(E_candidate(jj+1) temp1=; E_candidate(jj)=abs(E_candidate(jj); temp1=w_candidate(jj) E_candidate(jj); neww = neww temp1; end end %求出E的所有极值点 w_temp = ; neww=neww; neww=sortrows(neww,2); %对极值点根据对应E值排序,即误差最大的极值点 neww=neww; w1=neww(1,end-(M-3)/2:end); %根据l值选取最大的一些极值点 w_candidate = w(wcloc):(w(M+1)-w(wcloc)/10000:w(M+1)-1/10000; %在阻带内均分取10000个点,为了求出新的极值点 E_candidate = ; neww = ; for jj=1:length(w_candidate) E,H = Remez_H(w_candidate(jj),H_k_Matrix,wp,wc,M,a_k_Matrix,w); E_candidate = E_candidate E; end for jj=2:length(E_candidate)-1 if abs(E_candidate(jj)abs(E_candidate(jj-1) & abs(E_candidate(jj)abs(E_candidate(jj+1) temp1=; E_candidate(jj)=abs(E_candidate(jj); temp1=w_candidate(jj) E_candidate(jj); neww = neww temp1; end end %操作同通带内 w_temp = ; neww=neww; neww=sortrows(neww,2); neww=neww; w2=neww(1,end-(M-7)/2:end); w_temp = w1 w2 ws wp 0 pi; %更新新的w值,包括通带极值点,阻带极值%点,ws,wp,0,pi w = w_temp; w = sort(w); if p = pold %如果p值迭代后不改变,则退出迭代 break; end pold = p; countwhile = countwhile + 1;end%根据上一步求得w值求出a(0),a(1),a(M-1),并计算求出滤波器的w-H(w)PP = ;for j=1:M+1 temp = ; for k=1:M+1 if k=1 temp = temp 1; elseif k=M+1 temp = temp (-1)(j-1)/Remez_W(w(j),wp); else temp = temp cos(k-1)*w(j); end end PP = PP;temp; %构建书p255页的矩阵enda = inv(PP)*H_d_Matrix; a = a;PPPMatrix = ;for j=0:10000 w3 = pi*j/10000; PPP = Remez_P(a,w3,M); PPPMatrix = PPPMatrix PPP;endfor j=1:M+1 HH = Remez_P(a,w(j),M); H_MATRIX = H_MATRIX HH;end %求出ak,Hk,Hd,p由主函数调用function a_k_Matrix,H_k_Matrix,H_d_Matrix,p = Remez_DSP(M,w,wc,wp)%Function%index start with 0, so index+1 is used in Matlab as 1 standard.a_k_Matrix = ; %compute a_kfor k=0:M a_k=1; for i=0:M if(i=k) a_k=a_k*(1/(cos(w(i+1)-cos(w(k+1); end end a_k_Matrix = a_k_Matrix a_k;end%compute H_dH_d_Matrix = ;for k=0:M if abs(w(k+1) = wp H_d = 1; else H_d = 0; end H_d_Matrix = H_d_Matrix H_d;end%compute ppson = 0;for k=0:M-1 pson = pson + a_k_Matrix(k+1)*H_d_Matrix(k+1);endpmother = 0;for k=0:M pmother = pmother + (-1)(k)*a_k_Matrix(k+1)/Remez_W(w(k+1),wp);endp = pson/pmother;%compute H_kH_k_Matrix = ;for k=0:M H_k = H_d_Matrix(k+1)-(-1)(k)*p/Remez_W(w(k+1),wp); H_k_Matrix = H_k_Matrix H_k;end%求出E,H,由主函数调用function E,H = Remez_H(w,H_k_Matrix,wp,wc,M,a_k_Matrix,wo)%compute H H_son = 0; H_mother = 0; flag=0; for k2=0:M-1 if (cos(w)-cos(wo(k2+1)=0) H_son = H_son + (a_k_Matrix(k2+1)/(cos(w)-cos(wo(k2+1)*H_k_Matrix(k2+1); H_mother = H_mother + (a_k_Matrix(k2+1)/(cos(w)-cos(wo(k2+1); else flag=k2+1; break; end end if (flag=0) H = H_son/H_mother; else H=H_k_Matrix(flag); end %compute EE = Remez_W(w,wp)*(Remez_H_d(w,wc)-H);%求出矩阵的每一行,由主函数调用function PPP = Remez_P(a,w,M)PPP = 0;for j=1:M PPP = PPP + a(j)*cos(w*(