FIR等波纹滤波器优化设计

FIR等波纹滤波器优化设计[0,𝑤𝑝]段位通带[𝑤𝑠,π]阻带[𝑤𝑝,𝑤𝑠]过度带𝛿1是理想通带中能接受的振幅波动𝛿2是阻带中能接受的振幅波动通带波动dB值𝑅𝑝=−20𝑙𝑜𝑔101−𝛿11+𝛿2阻带波动dB值𝐴𝑠=−20𝑙𝑜𝑔10𝛿21+𝛿2M阶数设计方法1.窗函数法理想低通滤波器𝐻𝑑(𝑒𝑗𝑤)={𝑒−𝑗𝜔𝑛|𝑤|≤𝑤𝑐0𝑤𝑐𝑤≤𝑝𝑖对应时域ℎ𝑑(𝑛)=12π∫𝑒−𝑗𝑛0𝑤𝑒𝑗𝑛𝑤𝑑𝑤=sin(𝑤𝑐(𝑛−𝑛0))π(n−𝑛0)𝑤𝑐−𝑤𝑐无限长，非因果，无法实现h(n)=ℎ𝑑(𝑛)𝑅𝑁(𝑛)选取傅里叶级数的项数越多，引起的误差就越小另一种思路，不只简单地截取，对输出的傅里叶级数加权，就是用适当的窗函数形状进行加权2.频域样本法先在频率域修改滤波器的预期频率特性，使它不太理想，接近于实际的需要和可能，然后直接利用它的傅里叶反变换作为滤波器系数。最优设计法上面两种方法缺点：1不能精确地指定阻带频率𝑤𝑠和𝑤𝑝2不能同时控制𝛿1和𝛿23理想响应和实际响应之间的逼近误差在全频带区间不均匀，靠近边缘频率处误差较大，而在远离边缘处误差误差较小。（如果使误差均匀分布，就可以在满足相同技术指标的条件下，得到一个较低阶的滤波器。）A最小二乘法B最小最大波动法e(w)=𝑊𝑡(𝑤)[𝐴(𝑤)−𝐷(𝑤)]=𝑊𝑡(𝑤)∑𝑑(𝑛)[𝑐𝑜𝑠𝑛𝑤−𝐷(𝑤)]𝐿𝑛=0𝑊𝑡(w)幅特性1和1.01阻带0.011和0.001低通常取𝑊𝑡(w)=f(x)={𝛿1𝛿2⁄0≤𝑤≤𝑤𝑝1𝑤𝑝≤𝑤≤π问题找到滤波器系数向量d(n)，使得在w=[0,𝑤𝑝][𝑤𝑝,𝑤𝑠]频段内最大绝对误差|e(w)|为最小交替定理H(ejw)=∑ℎ(𝑛)𝑒−𝑗𝑛𝑤𝑀𝑛=−𝑀H(ejw)=h(0)+∑2ℎ(𝑛)𝑐𝑜𝑠𝑛𝑤𝑀𝑛=1设w是在[0,π]上的一个子集,比如是[0,𝑤𝑝]和[𝑤𝑝,𝑤𝑠]的合集。则A(w)=∑𝑑(𝑛)𝑐𝑜𝑠𝑛𝑤𝐿𝑛=0①是对D(w)的最佳逼近，当且仅当①式的误差函数e（w）是等波动的，并且该误差函数e（w）在w中至少有L+2个极值点。即在w中，存在着0≤𝑤1𝑤2⋯𝑤𝐿+2≤π其中包括𝑤𝑝和𝑤𝑠，在这些频率上误差正负交替，都达到同样的极值e(𝑤𝑖)=(−1)𝑖𝑒𝑚i=1,2,…,L+2②e(w)=𝑤𝑡(𝑤)[𝐴(𝑤)−𝐷(𝑤)]代入①𝑒(𝑤𝑖)=𝑊𝑡(𝑤)∑𝑑(𝑛)[𝑐𝑜𝑠𝑛𝑤−𝐷(𝑤)]𝐿𝑛=0代入②得到∑𝑑(𝑛)𝑐𝑜𝑠𝑛𝑤𝑖−(−1)𝑖𝑒𝑚𝑊𝑡(𝑤𝑖)=𝐷(𝑤𝑖)=𝐷𝑖𝐿𝑛=0L=4设计一个滤波器𝑤𝑝=0.2𝜋𝑅𝑝=0.25𝑑𝐵𝑤𝑠=0.3𝜋𝐴𝑠=50𝑑𝐵Remez交换算法Step1频率子集上均匀选取M个极值点𝑤0𝑤1…𝑤M计算ρ=∑𝑎𝑘𝐻𝑑(𝑒𝑗𝑤)𝑀𝑘=0∑(−1)𝑘𝑎𝑘𝑊(𝑤𝑘)⁄𝑀𝑘=0𝑎𝑘=∏1(𝑐𝑜𝑠𝑤𝑖−𝑐𝑜𝑠𝑤𝑘)𝑀𝑖=𝑜,𝑖≠𝑘Step2计算H(w)=∑[𝑎𝑘𝑐𝑜𝑠𝑤−𝑐𝑜𝑠𝑤𝑘]𝐻(𝑤𝑘)𝑀−1𝑘=0∑𝑎𝑘𝑐𝑜𝑠𝑤−𝑐𝑜𝑠𝑤𝑘𝑀−10其中H(𝑤𝑘)=𝐻𝑑(𝑤𝑘)−(−1)𝑘𝜌𝑊(𝑤𝑘)E(w)=W(w)[𝐻𝑑(w)−H(w)]在𝑤0,𝑤1,…,𝑤𝑀+1间检查是否有|E(w)|ρ若发现有新的极值点代替原有极值点，重复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+1lengthwc=0.5*pi;%理想低通滤波器截止频率%mainloop%%迭代开始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,pforindex=1:length(w)%此步为了保留通带截止频率wp阻带截止频率ws不改变ifw(index)==wpwploc=index;endifw(index)==wswcloc=index;endend%更新极值点w_candidate=0.0001:(w(wploc)-0.0001)/10000:w(wploc)-1/10000;%在通带内均分取10000个点,为了求出新的极值点E_candidate=[];neww=[];forjj=1:length(w_candidate)[E,H]=Remez_H(w_candidate(jj),H_k_Matrix,wp,wc,M,a_k_Matrix,w);E_candidate=[E_candidateE];End%求出10000个点所对应的E值forjj=2:length(E_candidate)-1ifabs(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=[newwtemp1];endend%求出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=[];forjj=1:length(w_candidate)[E,H]=Remez_H(w_candidate(jj),H_k_Matrix,wp,wc,M,a_k_Matrix,w);E_candidate=[E_candidateE];endforjj=2:length(E_candidate)-1ifabs(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=[newwtemp1];endend%操作同通带内w_temp=[];neww=neww';neww=sortrows(neww,2);neww=neww';w2=neww(1,end-(M-7)/2:end);w_temp=[w1w2wswp0pi];%更新新的w值,包括通带极值点,阻带极值%点,ws,wp,0,piw=w_temp;w=sort(w);ifp==pold%如果p值迭代后不改变,则退出迭代break;endpold=p;countwhile=countwhile+1;end%%根据上一步求得w值求出a(0),a(1),…,a(M-1),并计算求出滤波器的w-H(w)PP=[];forj=1:M+1temp=[];fork=1:M+1ifk==1temp=[temp1];elseifk==M+1temp=[temp(-1)^(j-1)/Remez_W(w(j),wp)];elsetemp=[tempcos((k-1)*w(j))];endendPP=[PP;temp];%构建书p255页的矩阵enda=inv(PP)*H_d_Matrix';a=a';PPPMatrix=[];forj=0:10000w3=pi*j/10000;PPP=Remez_P(a,w3,M);PPPMatrix=[PPPMatrixPPP];endforj=1:M+1HH=Remez_P(a,w(j),M);H_MATRIX=[H_MATRIXHH];end%%求出ak,Hk,Hd,p由主函数调用function[a_k_Matrix,H_k_Matrix,H_d_Matrix,p]=Remez_DSP(M,w,wc,wp)%Function%indexstartwith0,soindex+1isusedinMatlabas1standard.a_k_Matrix=[];%computea_kfork=0:Ma_k=1;fori=0:Mif(i~=k)a_k=a_k*(1/(cos(w(i+1))-cos(w(k+1))));endenda_k_Matrix=[a_k_Matrixa_k];end%computeH_dH_d_Matrix=[];fork=0:Mifabs(w(k+1))=wpH_d=1;elseH_d=0;endH_d_Matrix=[H_d_MatrixH_d];end%computeppson=0;fork=0:M-1pson=pson+a_k_Matrix(k+1)*H_d_Matrix(k+1);endpmother=0;fork=0:Mpmother=pmother+(-1)^(k)*a_k_Matrix(k+1)/Remez_W(w(k+1),wp);endp=pson/pmother;%computeH_kH_k_Matrix=[];fork=0:MH_k=H_d_Matrix(k+1)-(-1)^(k)*p/Remez_W(w(k+1),wp);H_k_Matrix=[H_k_MatrixH_k];end%%求出E,H,由主函数调用function[E,H]=Remez_H(w,H_k_Matrix,wp,wc,M,a_k_Matrix,wo)%computeHH_son=0;H_mother=0;flag=0;fork2=0:M-1if(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))));elseflag=k2+1;break;endendif(flag==0)H=H_son/H_mother;elseH=H_k_Matrix(flag);end%computeEE=Remez_W(w,wp)*(Remez_H_d(w,wc)-H);%%求出矩阵的每一行,由主函数调用functionPPP=Remez_P(a,w,M)PPP=0;forj