clear all; close all; wp = 0.5; ws = 0.75; ap = 0.01; as = 60; %%%% Butter Worth Filter %%%%%% % [n, wn] = buttord(wp,ws,ap,as); % [b,a] = butter(n,wn); %%%% Chebysheb Filter 1 %%%%%% % [n,Wn] = cheb1ord(wp,ws,ap,as); % [b,a] = cheby1(n,wp,Wn,'low'); %%%% Chebysheb Filter 2 %%%%%% % [n,Wn] = cheb2ord(wp,ws,ap,as); % [b,a] = cheby2(n,as,Wn,'low'); %%%% Chebysheb Filter 2 %%%%%% [n,Wn] = ellipord(wp,ws,ap,as); [b,a] = ellip(n,ap,as,Wn,'low'); xlswrite('D:\Blog\iir\iira.xlsx',-a); xlswrite('D:\Blog\iir\iirb.xlsx',b); %[r,p,k] = residuez(b,a) %[b,a] = residuez(r,p,k) %%%%% Direct form TF to delayed parallel form TF.. %[Bm,Am,FIR]=tf2delparf(b,a); %%%%% Direct form TF to non delayed parallel form TF.. [Bm,Am,FIR]=tf2parf(b,a); xlswrite('D:\Blog\iir\iirAm.xlsx',-Am); xlswrite('D:\Blog\iir\iirBm.xlsx',Bm); %zplane(b,a) [H,F] = freqz(b,a,512,1); figure(1) plot(F,20*log10(abs(H))) xlabel('Normalized Frequency') ylabel('Magnitude') %%%%% Direct form TF to Cascaded form TF.. %[b,a] = eqtflength(hn',1); % make lengths equal [z,p,k] = tf2zp(b,a); [sos,g] = zp2sos(z,p,k); xlswrite('D:\Blog\iir\sos.xlsx',sos); %%%%%% Verification %%%%%%%%% ts=1/100e3; t = 0:ts:2e-3; x = sin(2*pi*22e3*t) .* sin(2*pi*20e3*t); xlswrite('D:\Blog\iir\xin.xlsx',x); %%% Functions for filtration.. %y = filter(b,a,x); %%Normal Filter function %y=delparfilt(Bm,Am,FIR,x); %%% delayed parallel form y=parfilt(Bm,Am,FIR,x); %%% non-delayed parallel form figure(2) plot(t,y); xlabel('Time') ylabel('Amplitude') %%%%%% This part is for verification and can be ignored............... lpf_out = xlsread('D:\Blog\iir\op.xlsx'); lpf_out1 = downsample( lpf_out , 2); xlswrite('D:\Blog\iir\op_float.xlsx',lpf_out1); %fm = max(f); lpf_out1 = xlsread('D:\Blog\iir\op1.xlsx'); lpfom = lpf_out1/max(lpf_out1); ym = y/max(y); figure(3) plot(1:201,[zeros(4,1);lpf_out1(1:197,1)],1:201,y) legend('FPGA Output','MATLAB Output'); xlabel('Time') ylabel('Amplitude') xlim([0 201]) RMSE = norm(ym - lpfom')/norm(ym); function [Bm,Am,FIR]=tf2parf(B,A); [r,p,k]=residuez(B,A); % partial fraction expansion [Bm,Am,FIR]=rpk2parf(r,p,k); % recombine to second-order sections end function [Bm,Am,FIR]=tf2delparf(B,A); [r,p,k]=residued(B,A); % perform partial fraction expansion to the delayed form [Bm,Am,FIR]=rpk2parf(r,p,k); % recombine to second-order sections end function y=parfilt(Bm,Am,FIR,x); y=zeros(size(x)); s=size(Am); for k=1:s y=y+filter(Bm(:,k),Am(:,k),x); end if length(FIR)>0, y=y+filter(FIR,1,x); end; end function y=delparfilt(Bm,Am,FIR,x); NFIR=length(FIR); if NFIR==1, % if the FIR coeff is zero then there is no FIR part if FIR==0, NFIR=0; end; end; y=zeros(size(x)); s=size(Am); for k=1:s(2), % computing the parallel second-order sections y(NFIR+1:end)=y(NFIR+1:end)+filter(Bm(:,k),Am(:,k),x(1:end-NFIR)); end; if NFIR>0, % computing the FIR part y=y+filter(FIR,1,x); end;