crmz.m

有关matlab的电子书籍有一定的帮助希望有用

   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 frequecies 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);

%========================================================