gusucode.com > 信号处理工具箱 - signal源码程序 > signal\signal\signal\private\crmz.m
function [h,a,delta,not_optimal,iext,HH,EE, ... M_str,HH_str,h_str,Lf,Lb,Ls,Lc,A] = ... crmz( L, sym, Lfft, indx_edges, PLOTS, TRACE ) %CRMZ Complex Remez multiple-exchange filter design algorithm. % Designs FIR filters with arbitary magnitude and % phase specifications; reduces to the Remez (Parks-McClellan) % algorithm in the real case. % Authors: L. Karam, J. McClellan % Revised: 22-Oct-96, D. Orofino % % Copyright (c) 1988-98 by The MathWorks, Inc. % $Revision: 1.2 $ $Date: 1998/07/14 21:58:24 $ global DES_CRMZ WT_CRMZ GRID_CRMZ TGRID_CRMZ IFGRD_CRMZ J = 1i; N1 = 0; if nargin==1 & isstruct(L), if isstr(L.plot) & strcmp(L.plot,'plot-result'), N2 = N1 + L.L - 1; H = L.HH .* exp(-J*2*pi*TGRID_CRMZ*(N1+N2)/2); crmz_plresult( H, L.EE, L.iext, L.iext, L.indx_edges, N1, N2, L.delta, L.sym ); drawnow end return end if TRACE, disp(' '); disp(' ***********************************************'); disp(' ***** Complex Remez *****'); disp(' ***********************************************'); end if nargin<2, sym = 0; end; N2 = N1 + L - 1; % Last index of impulse response is_odd = rem(L,2); if (is_odd) Lf = (L+1)/2; Ws = ['W(:,2:Lf)']; hr_str = ['[a([Lc:-1:2])/2; a(1); a([2:Lc])/2]']; hi_str = ['[-a([Lb:-1:(Lc+1)])/2; 0; a([(Lc+1):Lb])/2]']; ph_str = ['ones(length(TGRID_CRMZ),1)']; else Lf = L / 2; Ws = ['W']; hr_str = ['[a([Lc:-1:1])/2; a([1:Lc])/2]']; hi_str = ['[-a([Lb:-1:(Lc+1)])/2; a([(Lc+1):Lb])/2]']; ph_str = ['exp(-J*pi*TGRID_CRMZ)']; end; Lc = Lf; Ls = L - Lf; switch sym case 0, Lb = L; M_str = ['[cos(W) sin(' Ws ')]']; h_str = [hr_str '+J*' hi_str]; HH_str = ['HH = fftshift( fft(hc,Lfft) ).*' ph_str ';']; case 1, % Even-symmetry: Lb = Lc; M_str = ['cos(W)']; h_str = [hr_str]; HH_str = ['HH=fft(hc,Lfft); HH((Lfft/2+2):Lfft)=[]; HH=HH.* ' ph_str ';']; case 2, % Odd-symmetry: Lb = Ls; Lc = 0; M_str = ['sin(' Ws ')']; h_str = ['J*' hi_str]; HH_str = ['HH=fft(hc,Lfft); HH((Lfft/2+2):Lfft)=[]; HH=HH.* ' ph_str ';']; end %------------------------ A = DES_CRMZ .* exp(J*2*pi*GRID_CRMZ*(N1+N2)/2); if PLOTS % clf; plot(GRID_CRMZ*2,abs(DES_CRMZ),'--') hold on plot(GRID_CRMZ*2, angle(DES_CRMZ),'-') hold off title('Desired Filter Frequency Response') xlabel('normalized frequency') ylabel('magnitude (dashed) and phase (solid)') drawnow end %----------------------- if TRACE, disp('Calculating initial solution .....'), end vec_edges = TGRID_CRMZ(indx_edges); edges = reshape(vec_edges,2,length(vec_edges)/2)'; [fext, iext] = crmz_guess(edges, GRID_CRMZ, Lb); it = 0; delta = 0.0; delta_old = -1; no_stp = 1; exactTol = 1e-15; % tolerance for exactness (arbitrary) % if the maximum error relative to the absolute maximum of the % desired function is less than exactTol, the filter designed % is considered 'exact' and the iteration is terminated successfully. while (no_stp) it = it + 1; delta_old = abs(delta); fext = GRID_CRMZ(iext); W = 2*pi*fext*([0:(Lf-1)]+(~is_odd)*0.5); Mb = eval(M_str); M = [ Mb ((-1).^[0:Lb]')./WT_CRMZ(iext) ]; a = M \ A(iext); delta = a(Lb+1); h = eval(h_str); hc = crmz_rotate(zeropad(h,Lfft-L), -Ls); eval(HH_str); W = 2*pi*vec_edges*([0:(Lf-1)]+(~is_odd)*0.5); Mb = eval(M_str); HH(indx_edges) = Mb * a(1:Lb); EE = WT_CRMZ .* (A - HH(IFGRD_CRMZ)); EE(iext) = delta*( (-1) .^ [2:length(iext)+1] ); [jext,EEj] = crmz_find(EE,iext); if TRACE, fprintf('(%3.0f): ', it) fprintf('delta = %8.5f+j%8.5f, ', real(delta), imag(delta)) fprintf('|delta| = %5.2e, ', abs(delta)); fprintf('e_max = %5.2e\n', max(abs(EE))); end if PLOTS, crmz_plerror(EE,HH,iext,jext,indx_edges,GRID_CRMZ,delta,sym); drawnow end if all(iext == jext') | (max(abs(EE))/max(abs(A))<exactTol) no_stp = 0; end iext = jext'; end %%% end of Complex Remez algorithm %%%% %%%%% Assessing Optimality of the Solution %%%% e_max = max(abs(EE)); tlr = abs(delta)/100; %% Needed due to limited machine accuracy %% if (e_max <= (abs(delta)+tlr)) | (e_max/max(abs(A))<exactTol) not_optimal = 0; if TRACE, fprintf('Optimal Solution Obtained !\n') fprintf('(%3.0f): ', it) fprintf('delta = %8.5f+j%8.5f, ', real(delta), imag(delta)) fprintf('|delta| = %5.2e, ', abs(delta)); fprintf('e_max = %5.2e\n',max(abs(EE))); end else not_optimal = 1; if TRACE crmz_chckopt( EE, delta, tlr, Lfft ) fprintf('(%3.0f): ', it) fprintf('delta = %8.5f+j%8.5f, ', real(delta), imag(delta)) fprintf('|delta| = %5.2e, ', abs(delta)); fprintf('e_max = %5.2e\n',max(abs(EE))); end end %%%%%%%%%%%%% End Complex Remez Stage %%%%%%%%%%%%%%%%%%%%%%%% %======================================================== function crmz_chckopt( EE, delta, tlr, Lfft ) % Checks degree of suboptimality of solution not optimal % on desired frequency bands. % subint = find( abs(EE) <= (abs(delta)+tlr) ); rel_size = round((length(subint)/Lfft)*100); disp(['Optimal solution found on a frequency subset comprising ' ... sprintf('%3.0f percent \n',rel_size) 'of the original frequency set']) %======================================================== function [fnew,enew] = crmz_find( error, fold ) % Construct new subset of points using exchange rules % Authors: LJ Karam and JH McClellan % Reference: "Complex Chebyshev approximation for FIR filter design", % IEEE Trans. on Circuits and Systems II, March 1995. fold = fold(:); Nx = length( fold ); Ngrid = length( error ); if error(fold(1))~=0 sgn_error = error(fold(1))/abs(error(fold(1))); else sgn_error = 1; end error = real( conj(sgn_error)*error ); delta = min( abs( error(fold) ) ); %--- present value of delta up = sign(error(fold(1))) > 0; if up, fence = [ [1;fold(1:2:Nx)] [fold(1:2:Nx);Ngrid] ]; else fence = [ [1;fold(2:2:Nx)] [fold(2:2:Nx);Ngrid] ]; end Lf = length(fence(:,1)); emn = zeros(Lf,1); imn = emn; emx = emn; imx = emn; for i=1:Lf, [emn(i),imn(i)] = min( error(fence(i,1):fence(i,2) )); imn(i) = imn(i) + fence(i,1) -1; end imn(find((emn > -delta))) = []; fence = [ [1;imn] [imn;Ngrid] ]; Lf = length(fence(:,1)); for i=1:Lf [emx(i),imx(i)] = max( error(fence(i,1):fence(i,2) )); imx(i) = imx(i) + fence(i,1) -1; end imx(find((emx < delta))) = []; fnew = sort([imx;imn]); Nf = length(fnew); if ( Nf > Nx) if ( abs(error(fnew(1))) >= abs(error(fnew(Nf))) ) fnew = fnew(1:Nx); else fnew = fnew(Nf-Nx+1:Nf); end end fnew = fnew'; enew = error(fnew); %======================================================== function [fext,iext] = crmz_guess( edges, grid, nfcns ) %CRMZ_GUESS generate initial guess of "EXTREMAL FREQS" % for remez exchange algorithm % usage: % [Fext,Iext] = crmz_guess( Edges, Grid, Nfcns ) % % Edges: band-edges moved onto the grid (Revised edges) % Grid: frequencies on the fine grid % Nfcns: number of basis functions in the approx % Fext: initial guess of "extremal frequencies" % Iext: indices for the "extremal frequencies" % TOL = 5*eps; next = nfcns+1; Nbands = length(edges(:,2)); tt = edges; merged = tt; if Nbands > 1, jkl = find( abs(tt(1:(Nbands-1),2)-tt(2:Nbands,1)) > TOL ); merged = [ [tt(1,1); tt(jkl+1,1)] , [tt(jkl,2); tt(Nbands,2)] ]; end Nbands = length(merged(:,1)); bw = merged(:,2)-merged(:,1); if any(bw < 0), edges error('Internal error: obtained a negative bandwidth'); end percent_bw = bw / sum(bw); fext = zeros(next,1); n = 0; i = 1; while (n < next) & (i <= Nbands), nfreqs_i = min( next-n, ceil( percent_bw(i)*next ) ); if( nfreqs_i == 0 ) n=n+1; fext(n) = merged(i,1); else fext(n+[1:nfreqs_i]) = ... linspace(merged(i,1),merged(i,2),nfreqs_i); n = n+nfreqs_i; end i = i+1; end iext = 0*fext; for i=1:next [tt,iext(i)] = min( abs(fext(i)-grid) ); end if any(diff(iext) == 0), error('Internal error: two extremal frequencies at the same grid point') end fext = grid(iext); %======================================================== function rotated = crmz_rotate(x,num_places) %CRMZ_ROTATE circular shift of columns in a matrix % crmz_rotate(V,r) % circularly shifts the elements in the columns of V % by r places right (r>0); or r places left (r<0). % (Right is down; left is up.) % If the input is a row or column vector, the shift is % performed on the vector. % If the input is a signal matrix, each column is shifted % [M,N] = size(x); if M > 1, % rotate columns num_places = mod(num_places,M); % make num_places in range [0,M-1] rotated = [ x(M-num_places+1:M,:); x(1:M-num_places,:) ]; else % Assume N > 1 % rotate row vector num_places = mod(num_places,N); % make num_places in range [0,N-1] rotated = [ x(N-num_places+1:N) x(1:N-num_places) ]; end %======================================================== function crmz_plerror2( H, EE, iext, jext, indx_edges, N1, N2, delta, sym ) % crmz_plerror2 % plot error magnitude and real phase-rotated error. global DES_CRMZ WT_CRMZ GRID_CRMZ TGRID_CRMZ IFGRD_CRMZ J = 1i; EE1 = EE(iext(1)); EE_pl = []; grd_pl = []; i = 1; l = 0; while (i < (length(indx_edges)-1)) lo = l + 1; l = l + length([indx_edges(i):indx_edges(i+1)]); grd_pl = [grd_pl;GRID_CRMZ(lo:l);nan]; EE_pl = [EE_pl;EE(lo:l);nan]; i = i + 2; end; lo = l + 1; l = l + length([indx_edges(i):indx_edges(i+1)]); grd_pl = [grd_pl;GRID_CRMZ(lo:l)]; EE_pl = [EE_pl;EE(lo:l)]; EEE_pl = [real(EE_pl * conj(EE1/abs(EE1)) )-3*abs(delta),abs(EE_pl)]; EEE = [real(EE * conj(EE1/abs(EE1)) )-3*abs(delta), abs(EE)]; delr = real(delta); deli = imag(delta); axis('auto') plot( grd_pl*2, EEE_pl, '-', [GRID_CRMZ(1), 0.5]*2, abs(delta)*[1;1], '--',... [GRID_CRMZ(1), 0.5]*2, abs(delta)*[1 -1;1 -1]-3*abs(delta), ':',... GRID_CRMZ(jext)*2, EEE(jext,:), 'o', GRID_CRMZ(iext)*2, EEE(iext,:), '+') V = axis; axis([-1*(~sym) 1 V(3) V(4)]) %======================================================== function crmz_plresult( H, EE, iext, jext, indx_edges, N1, N2, delta, sym ) % crmz_plresult % % Generates subplots of final result with % 1. Magnitude response % 2. Passband Phase % 3. Error magnitude and Real rotated error % 4. Phase error in passband % Uses grp_delay.m global DES_CRMZ WT_CRMZ GRID_CRMZ TGRID_CRMZ IFGRD_CRMZ J = 1i; % Compute EE_pl for plotting error in bands: EE_pl = []; grd_pl = []; i = 1; l = 0; while (i < (length(indx_edges)-1)) lo = l + 1; l = l + length([indx_edges(i):indx_edges(i+1)]); grd_pl = [grd_pl;GRID_CRMZ(lo:l);nan]; EE_pl = [EE_pl;EE(lo:l);nan]; i = i + 2; end lo = l + 1; l = l + length([indx_edges(i):indx_edges(i+1)]); grd_pl = [grd_pl;GRID_CRMZ(lo:l)]; EE_pl = [EE_pl;EE(lo:l)]; E_pl = EE_pl .* exp(-2i*pi*grd_pl*(N1+N2)/2); % Needed for plot hfig = figure; subplot(111) axis('auto'); subplot(221) dbMax = 20*log(max(abs(H(:)))); %dbMin = 20*log(min(abs(H(:)))); %dbRange = dbMax - dbMin; dbRange = 80; plot(TGRID_CRMZ*2, db(abs(H(:)),dbRange,dbMax), '-') ph = gca; set(ph,'box','off'); xlabel('normalized frequency'); ylabel('magnitude (log scale)') subplot(222); if length(indx_edges) > 2 pas_edges = TGRID_CRMZ( indx_edges(3:4) ); else pas_edges = TGRID_CRMZ( indx_edges(1:2) ); end pas_edges = pas_edges'; pas_edges = pas_edges(:); i = 1; grd_pl = []; aH_pl = []; aD_pl = []; while (i < (length(pas_edges)-1)) psb = find((GRID_CRMZ >= pas_edges(i)) & (GRID_CRMZ <= pas_edges(i+1))); grd_pl = [grd_pl;GRID_CRMZ(psb);nan]; aH_pl = [aH_pl;angle(H(IFGRD_CRMZ((psb))));nan]; aD_pl = [aD_pl;angle(DES_CRMZ(psb));nan]; i = i + 2; end psb = find((GRID_CRMZ >= pas_edges(i)) & (GRID_CRMZ <= pas_edges(i+1))); grd_pl = [grd_pl;GRID_CRMZ(psb)]; aH_pl = [aH_pl;angle(H(IFGRD_CRMZ((psb))))]; aD_pl = [aD_pl;angle(DES_CRMZ(psb))]; plot(grd_pl*2, aH_pl, '-') hold on plot(grd_pl*2, aD_pl, '--') ph = gca; set(ph,'box','off'); xlabel('normalized frequency'), ylabel('passband phase (radians)') V = axis; V(1) = pas_edges(1); V(2) = pas_edges(length(pas_edges)); axis([2*V(1) 2*V(2) V(3) V(4)]); x = (V(1)+V(2))/(2.0); text(2*x,V(4),'( -- ideal )') hold off subplot(224) e_angle = mod(aD_pl,2*pi)-mod(aH_pl,2*pi); plot(grd_pl*2,e_angle) xlabel('normalized frequency'), ylabel('passband phase error') ph = gca; set(ph,'box','off'); V = axis; V(1) = pas_edges(1); V(2) = pas_edges(length(pas_edges)); axis([2*V(1) 2*V(2) V(3) V(4)]); subplot(223) crmz_plerror2( H, EE, iext, jext, indx_edges, N1, N2, delta, sym ) ph = gca; set(ph,'box','off'); axis('off') text(.95,-4.5*abs(delta),'1') text(-0.0108,-4.5*abs(delta),'0') if (~sym) text(-0.52*2,-4.5*abs(delta),'-1') str = sprintf('%5.4f',abs(delta)); text(-0.75*2,abs(delta),str,'fontsize',10); text(-0.75*2,-2*abs(delta),str,'fontsize',10); str = sprintf('%5.4f',-abs(delta)); text(-0.78*2,-4*abs(delta),str,'fontsize',10); text(-0.22*2,1.6*abs(delta),'Error magnitude','fontsize',10); text(-0.25*2,-1.5*abs(delta),'Real rotated error','fontsize',10); else str = sprintf('%5.4f',abs(delta)); text(-0.14*2,abs(delta),str,'fontsize',10); text(-0.14*2,-2*abs(delta),str,'fontsize',10); str = sprintf('%5.4f',-abs(delta)); text(-0.16*2,-4*abs(delta),str,'fontsize',10); text(0.155*2,1.6*abs(delta),'Error magnitude','fontsize',10); text(0.125*2,-1.5*abs(delta),'Real rotated error','fontsize',10); end xlabel('normalized frequency') %======================================================== function crmz_plerror(EE,HH,iext,jext,indx_edges,grd,delta,sym) %crmz_plerror % plot error magnitude, real phase rotated error, % error real part, error imaginary part. % clf; EE1 = EE(iext(1)); EE_pl = []; grd_pl = []; i = 1; l = 0; while (i < (length(indx_edges)-1)) lo = l + 1; l = l + length([indx_edges(i):indx_edges(i+1)]); grd_pl = [grd_pl;grd(lo:l);nan]; EE_pl = [EE_pl;EE(lo:l);nan]; i = i + 2; end; lo = l + 1; l = l + length([indx_edges(i):indx_edges(i+1)]); grd_pl = [grd_pl;grd(lo:l)]; EE_pl = [EE_pl;EE(lo:l)]; EEE_pl = [real(EE_pl * conj(EE1/abs(EE1)) )-3*abs(delta),... abs(EE_pl), real(EE_pl)-6*abs(delta),... imag(EE_pl)-9*abs(delta)]; %-- plottting offset EEE = [real(EE * conj(EE1/abs(EE1)) )-3*abs(delta),... abs(EE), real(EE)-6*abs(delta),... imag(EE)-9*abs(delta)]; %-- plotting offset delr = real(delta); deli = imag(delta); axis('auto') plot( grd_pl*2, EEE_pl, '-', [grd(1), 0.5]*2, abs(delta)*[1;1], '--',... [grd(1), 0.5]*2, abs(delta)*[1 -1;1 -1]-3*abs(delta), ':',... [grd(1), 0.5]*2, delr*[1 -1;1 -1]-6*abs(delta),':',... [grd(1), 0.5]*2, deli*[1 -1;1 -1]-9*abs(delta),':',... grd(jext)*2, EEE(jext,:), 'o', grd(iext)*2, EEE(iext,:), '+') V = axis; axis([-1*(~sym) 1 V(3) V(4)]) xlabel('normalized frequency') title(['Weighted Error: 1.magnitude 2.real-rotated 3.real part ' ... '4.imag. part']) %======================================================== function z_out = zeropad(x,num_of_zeros) %ZEROPAD % zeropad(x,L) will append L zeros to each column of x. % [M,N] = size(x); if M == 1, % input is a row vector z_out = [ x zeros(1,num_of_zeros) ]; else % input is a matrix or column z_out = [ x; zeros(num_of_zeros,N) ]; end %======================================================== function y = db( x, dBrange, dBmax ) %DB Convert an array to decibels % Y=DB( X, dbRANGE, dbMAX ) will compute 20 Log(X) % and then scale or clip the result so that % the minimum dB level is dbMAX-dbRANGE. % ex: db(X, 80, 0) gives the range 0 to -80 dB % Y = dB( X, dbRANGE ) defaults to dBmax = 0 % % NOTE: on some machines log(0) gives an error % and aborts this function (esp. microVAX) [M,N] = size(x); if dBrange <= 0, error('dBrange must be > 0'); end if nargin == 2, dBmax = 0; end y = abs(x); ymax = max(y)/10.0^(dBmax/20); if M == 1, %-- input is a row y = y/ymax; else y = x ./ ymax(ones(M,1),:); end thresh = 10.0^((dBmax-dBrange)/20); y = y.*(y>thresh) + thresh.*(y<=thresh); y = 20.0*log10(y); %========================================================