www.gusucode.com > sigdemos 工具箱matlab源码程序 > sigdemos/CZTExample.m
function CZTExample(action,s)
%CZTExample Example showing the FFT and CZT in the Signal Processing Toolbox.
%
% This example shows how to compare the transform of a filter using two
% different algorithms. They are
% FFT - the Fast Fourier Transform
% CZT - the Chirp-Z Transform
%
% The upper plot shows the z-plane, and the lower plot shows the transform of
% a bandpass elliptic digital filter on the points within the wedge shown on
% the upper plot. The filter has a passband from .4 to .7 of the Nyquist
% frequency (Nyquist = Fs/2).
%
% The FFT computes the Z-transform on equally spaced points around the unit circle.
%
% The CZT computes the Z-transform on a spiral or "chirp" contour. The
% contour is defined by initial frequency Fmin and radius R1, and final
% frequency Fmax and radius R2. Fs is the sampling frequency.
%
% Fmin and Fmax define a "wedge" of the unit circle.
%
% Npoints is the number of Z-transform points computed on the unit circle in
% the wedge defined by Fmin and Fmax.
%
% With FFT, the length of the transform is Npoints*Fs/(Fmax-Fmin), which
% computes Npoints points in the range Fmin to Fmax. If you are interested in
% a small frequency range, the CZT is more efficient because it "zooms in" on
% the range you are interested in and only computes Npoints.
%
% See also FILTDEMO, MODDEMO.
% Copyright 1988-2012 The MathWorks, Inc.
% Possible actions:
% initialize
% Fs
% Wmin
% Wmax
% points
% radius1
% radius2
% design
% button callbacks:
% radio
% infobtn
% close
if nargin<1,
action='initialize';
end;
if nargin <= 1,
feval(action);
else
feval(action,s);
end
function design
% evaluate fft or czt
set(gcf,'Pointer','watch');
ud=get(gcf,'UserData');
hndlList=ud.hndlList;
b=ud.filter.num;
a=ud.filter.den;
zplaneHndl = findobj(hndlList,'Tag','z-Plane');
responseHndl = findobj(hndlList,'Tag','response');
FsHndl = findobj(hndlList,'Tag','Fs');
Fs = get(FsHndl,'UserData');
WminHndl = findobj(hndlList,'Tag','Fmin');
Wmin = get(WminHndl,'UserData');
WmaxHndl = findobj(hndlList,'Tag','Fmax');
Wmax = get(WmaxHndl,'UserData');
pointsHndl = findobj(hndlList,'Tag','points');
Npoints = get(pointsHndl,'UserData');
radius1Hndl = findobj(hndlList,'Tag','Radius1');
R1 = get(radius1Hndl,'UserData');
radius2Hndl = findobj(hndlList,'Tag','Radius2');
R2 = get(radius2Hndl,'UserData');
btn1Hndl = findobj(hndlList,'Tag','FFT');
btn2Hndl = findobj(hndlList,'Tag','CZT');
closeHndl = findobj(hndlList,'Tag','Close');
set(gcf,'NextPlot','add')
if (get(btn1Hndl,'Value')==1), % FFT button checked?
n = (ceil(Npoints*Fs / max(Wmax-Wmin,eps)));
if n>2048, n = 2^nextpow2(n); end
ww = (0:min(n,400)-1)'/min(n,400)*2*pi;
axes(zplaneHndl)
cla(zplaneHndl)
[domainHndl,h2,h3] = zplane(exp(j*ww),[],zplaneHndl);
if length(h3)>1, delete(h3(2:length(h3))), end
rcolor = get(gcf,'DefaultAxesColorOrder');
rcolor = rcolor(min(2,size(rcolor,1)),:);
domainHndl(2) = line('XData',[0 cos(Wmin*2*pi/Fs)],...
'Color',rcolor,'YData',[0 sin(Wmin*2*pi/Fs)]);
domainHndl(3) = line('XData',[0 cos(Wmax*2*pi/Fs)],...
'Color',rcolor,'YData',[0 sin(Wmax*2*pi/Fs)]);
set(closeHndl,'UserData',domainHndl)
%title('Domain of FFT')
title('z-Plane')
axes(responseHndl)
F = fft(b,n)./fft(a,n);
nreplicas=ceil(Wmax/Fs)-floor(Wmin/Fs);
f=linspace(floor(Wmin/Fs)*Fs,ceil(Wmax/Fs)*Fs-1/n,nreplicas*n);
F=repmat(F,1,nreplicas);
subplot(212);
plot(f,20*log10(abs(F)),'.');
set(gca,'Tag','response');
if n<=Npoints, %(Wmin == 0)&(Wmax>=Fs*(1-1/n))
title(sprintf('%g point FFT of elliptic bandpass filter',n));
else
title(sprintf('Close-up of %g point FFT of elliptic bandpass filter',n));
end
elseif (get(btn2Hndl,'Value')==1),% CZT button checked.
M = Npoints;
A = R1*exp(j*Wmin*2*pi/Fs);
W = ( R2/R1 )^(-1/(M-1)) * exp(-j*(Wmax-Wmin)*2*pi/Fs/(M-1)) ;
axes(zplaneHndl)
z = A*W.^(-(0:M-1)');
cla(zplaneHndl)
[domainHndl,h2,h3] = zplane(z,[],zplaneHndl);
if length(h3)>1, delete(h3(2:length(h3))), end
rcolor = get(gcf,'DefaultAxesColorOrder');
rcolor = rcolor(min(2,size(rcolor,1)),:);
domainHndl(2) = line('XData',[0 R1*cos(Wmin*2*pi/Fs)],...
'Color',rcolor,'YData',[0 R1*sin(Wmin*2*pi/Fs)]);
text(R1/2*cos(Wmin*2*pi/Fs),R1/2*sin(Wmin*2*pi/Fs),' R1',...
'VerticalAlignment','top');
domainHndl(3) = line('XData',[0 R2*cos(Wmax*2*pi/Fs)],...
'Color',rcolor,'YData',[0 R2*sin(Wmax*2*pi/Fs)]);
text(R2/2*cos(Wmax*2*pi/Fs),R2/2*sin(Wmax*2*pi/Fs),'R2 ',...
'HorizontalAlignment','right','VerticalAlignment','top');
set(closeHndl,'UserData',domainHndl)
%title('Domain of CZT')
title('z-Plane')
axes(responseHndl)
w = unwrap(angle(z));
w = linspace(Wmin,Wmax,M)*2*pi/Fs;
Ga = czt(a,M,W,A);
Gb = czt(b,M,W,A);
% If any elements of Ga or Gb are zero set to eps
zeroIndxs_a=find(Ga==0);
zeroIndxs_b=find(Gb==0);
Ga(zeroIndxs_a)=eps;
Gb(zeroIndxs_b)=eps;
F = Gb./Ga;
cla(responseHndl)
subplot(212);
plot(w*Fs/2/pi,20*log10(abs(F)),'.')
set(gca,'Tag','response');
title(sprintf('%g point CZT of elliptic bandpass filter',M));
end
xlabel('Frequency')
ylabel('Magnitude (dB)')
set(gca,'XLim',[Wmin Wmax])
ylim = get(gca,'YLim');
set(gca,'YLim',[max(-350,ylim(1)) ylim(2)])
set(gcf,'Pointer','arrow');
function initialize
shh = get(0,'ShowHiddenHandles');
set(0,'ShowHiddenHandles','on')
figNumber=figure( ...
'Name','CZT and FFT Demo', ...
'HandleVisibility','callback',...
'IntegerHandle','off',...
'NumberTitle','off');
%[b,a]=ellip(9,3,50,[.4 .7]); % design filter - store in btn2 userdata
% inline filter coeffs for speed:
num = [ 1.036215553331465e-02
2.103525287321029e-02
4.618180706244246e-02
6.626949942636884e-02
1.047645705817928e-01
1.135461439917620e-01
1.182161372812089e-01
9.184711304310156e-02
5.839803125783760e-02
-5.115907697472721e-13
-5.839803125883236e-02
-9.184711304391158e-02
-1.182161372818484e-01
-1.135461439921812e-01
-1.047645705820628e-01
-6.626949942649851e-02
-4.618180706250907e-02
-2.103525287323005e-02
-1.036215553332198e-02]';
den = [ 1.000000000000000e+00
2.616784951166920e+00
9.010794557412463e+00
1.664491445555174e+01
3.429175523852171e+01
4.935300363227517e+01
7.526259700870287e+01
8.775364685258239e+01
1.063920847920586e+02
1.018212388254725e+02
1.008239318769620e+02
7.874570918388088e+01
6.397972011971305e+01
3.962383394162555e+01
2.603920255580061e+01
1.187892654523706e+01
6.065579909513929e+00
1.633745653505110e+00
5.883332413893858e-01]';
ud.filter.den=den;
ud.filter.num=num;
%==================================
% Set up the image axes
axes( ...
'Units','normalized', ...
'Position',[0.10 0.35 0.60 0.6], ...
'XTick',[],'YTick',[], ...
'Box','on');
set(figNumber,'DefaultAxesPosition',[0.10 0.1 0.60 0.80])
zplaneHndl = subplot(2,1,1);
set(zplaneHndl,'Tag','z-Plane');
set(gca, ...
'Units','normalized', ...
'XTick',[],'YTick',[], ...
'Box','on');
responseHndl = subplot(2,1,2);
set(gca, ...
'Units','normalized', ...
'XTick',[],'YTick',[], ...
'Box','on');
set(responseHndl,'Tag','response');
%====================================
% Information for all buttons (and menus)
labelColor=[0.8 0.8 0.8];
yInitPos=0.90;
menutop=0.95;
btnTop = 0.6;
top=0.75;
left=0.785;
btnWid=0.175;
btnHt=0.06;
textHeight = 0.05;
textWidth = 0.07;
% Spacing between the button and the next command's label
spacing=0.019;
%====================================
% The CONSOLE frame
frmBorder=0.019; frmBottom=0.04;
frmHeight = 0.92; frmWidth = btnWid;
yPos=frmBottom-frmBorder;
frmPos=[left-frmBorder yPos frmWidth+2*frmBorder frmHeight+2*frmBorder];
h=uicontrol( ...
'Style','frame', ...
'Units','normalized', ...
'Position',frmPos, ...
'BackgroundColor',[0.5 0.5 0.5]);
%====================================
% fft radio button
btnTop = menutop-spacing;
btnNumber=1;
yPos=btnTop-(btnNumber-1)*(btnHt+spacing);
labelStr='FFT';
callbackStr='CZTExample(''radio'',1);';
% Generic button information
btnPos=[left yPos-btnHt btnWid btnHt];
btn1Hndl=uicontrol( ...
'Style','radiobutton', ...
'Units','normalized', ...
'Position',btnPos, ...
'String',labelStr, ...
'Value',0, ...
'Callback',callbackStr,...
'Tag',labelStr,...
'TooltipString','Fast Fourier Transform');
%====================================
% czt radio button
btnTop = menutop-spacing;
btnNumber=2;
yPos=btnTop-(btnNumber-1)*(btnHt+spacing);
labelStr='CZT';
callbackStr='CZTExample(''radio'',2);';
% Generic button information
btnPos=[left yPos-btnHt btnWid btnHt];
btn2Hndl=uicontrol( ...
'Style','radiobutton', ...
'Units','normalized', ...
'Position',btnPos, ...
'String',labelStr, ...
'Value',1, ...
'Callback',callbackStr,...
'Tag',labelStr,...
'TooltipString','Chirp z-Transform');
yPos = yPos - spacing;
%===================================
% Sampling Frequency
top = yPos - btnHt - spacing;
labelWidth = frmWidth-textWidth-.01;
labelBottom=top-textHeight;
labelLeft = left;
labelRight = left+btnWid;
labelPos = [labelLeft labelBottom labelWidth textHeight];
h = uicontrol( ...
'Style','text', ...
'Units','normalized', ...
'Position',labelPos, ...
'HorizontalAlignment','left', ...
'String','Fs', ...
'Interruptible','off', ...
'BackgroundColor',[0.5 0.5 0.5], ...
'ForegroundColor','white',...
'TooltipString','Sampling Frequency');
% Text field
textPos = [labelRight-textWidth labelBottom textWidth textHeight];
callbackStr = 'CZTExample(''Fs'')';
FsHndl = uicontrol( ...
'Style','edit', ...
'Units','normalized', ...
'Position',textPos, ...
'HorizontalAlignment','right', ...
'BackgroundColor','white', ...
'ForegroundColor','black', ...
'String','1000','UserData',1000, ...
'Callback',callbackStr,...
'Tag','Fs');
%===================================
% Wmin frequency (1) label and text field
labelBottom=top-2*textHeight-spacing;
labelLeft = left;
labelPos = [labelLeft labelBottom labelWidth textHeight];
h = uicontrol( ...
'Style','text', ...
'Units','normalized', ...
'Position',labelPos, ...
'String','Fmin', ...
'HorizontalAlignment','left', ...
'Interruptible','off', ...
'BackgroundColor',[0.5 0.5 0.5], ...
'ForegroundColor','white',...
'TooltipString','Lower Frequency Bound');
% Text field
textPos = [labelRight-textWidth labelBottom textWidth textHeight];
callbackStr = 'CZTExample(''Wmin'')';
WminHndl = uicontrol( ...
'Style','edit', ...
'Units','normalized', ...
'Position',textPos, ...
'HorizontalAlignment','center', ...
'BackgroundColor','white', ...
'ForegroundColor','black', ...
'String','200','UserData',200, ...
'Callback',callbackStr,...
'Tag','Fmin');
%===================================
% Wmax frequency label and text field
labelBottom=top-3*textHeight-2*spacing;
labelLeft = left;
labelPos = [labelLeft labelBottom labelWidth textHeight];
h = uicontrol( ...
'Style','text', ...
'Units','normalized', ...
'Position',labelPos, ...
'String','Fmax', ...
'HorizontalAlignment','left', ...
'Interruptible','off', ...
'BackgroundColor',[0.5 0.5 0.5], ...
'ForegroundColor','white',...
'TooltipString','Upper Frequency Bound');
% Text field
textPos = [labelRight-textWidth labelBottom textWidth textHeight];
callbackStr = 'CZTExample(''Wmax'')';
WmaxHndl = uicontrol( ...
'Style','edit', ...
'Units','normalized', ...
'Position',textPos, ...
'HorizontalAlignment','center', ...
'BackgroundColor','white', ...
'ForegroundColor','black', ...
'String','350','UserData',350, ...
'Callback',callbackStr,...
'Tag','Fmax');
%===================================
% Number of points label and text field
labelBottom=top-4*textHeight-3*spacing;
labelLeft = left;
labelPos = [labelLeft labelBottom labelWidth textHeight];
h = uicontrol( ...
'Style','text', ...
'Units','normalized', ...
'Position',labelPos, ...
'String','Npoints', ...
'HorizontalAlignment','left', ...
'Interruptible','off', ...
'BackgroundColor',[0.5 0.5 0.5], ...
'ForegroundColor','white',...
'TooltipString','Number of Points of Interest');
% Text field
textPos = [labelRight-textWidth labelBottom textWidth textHeight];
callbackStr = 'CZTExample(''points'')';
pointsHndl = uicontrol( ...
'Style','edit', ...
'Units','normalized', ...
'Position',textPos, ...
'HorizontalAlignment','center', ...
'BackgroundColor','white', ...
'ForegroundColor','black', ...
'String','300','UserData',300, ...
'Callback',callbackStr,...
'Tag','points');
%===================================
% Radius (1) label and text field
labelBottom=top-5*textHeight-4*spacing;
labelLeft = left;
labelPos = [labelLeft labelBottom labelWidth textHeight];
h = uicontrol( ...
'Style','text', ...
'Units','normalized', ...
'Position',labelPos, ...
'String','R1', ...
'HorizontalAlignment','left', ...
'Interruptible','off', ...
'BackgroundColor',[0.5 0.5 0.5], ...
'ForegroundColor','white',...
'TooltipString','Radius 1');
% Text field
textPos = [labelRight-textWidth labelBottom textWidth textHeight];
callbackStr = 'CZTExample(''radius1'')';
radius1Hndl = uicontrol( ...
'Style','edit', ...
'Units','normalized', ...
'Position',textPos, ...
'HorizontalAlignment','center', ...
'BackgroundColor','white', ...
'ForegroundColor','black', ...
'String','1','UserData',1, ...
'Tag','Radius1',...
'Callback',callbackStr);
%===================================
% Radius (2) label and text field
labelBottom=top-6*textHeight-5*spacing;
labelLeft = left;
labelPos = [labelLeft labelBottom labelWidth textHeight];
h = uicontrol( ...
'Style','text', ...
'Units','normalized', ...
'Position',labelPos, ...
'String','R2', ...
'HorizontalAlignment','left', ...
'Interruptible','off', ...
'BackgroundColor',[0.5 0.5 0.5], ...
'ForegroundColor','white',...
'TooltipString','Radius 2');
% Text field
textPos = [labelRight-textWidth labelBottom textWidth textHeight];
callbackStr = 'CZTExample(''radius2'')';
radius2Hndl = uicontrol( ...
'Style','edit', ...
'Units','normalized', ...
'Position',textPos, ...
'HorizontalAlignment','center', ...
'BackgroundColor','white', ...
'ForegroundColor','black', ...
'String','1','UserData',1, ...
'Tag','Radius2',...
'Callback',callbackStr);
%====================================
% The INFO button
labelStr='Info';
callbackStr='CZTExample(''infobtn'')';
helpHndl=uicontrol( ...
'Style','pushbutton', ...
'Units','normalized', ...
'Position',[left frmBottom+btnHt+spacing btnWid btnHt], ...
'String',labelStr, ...
'Callback',callbackStr,...
'Tag',labelStr);
%====================================
% The CLOSE button
labelStr='Close';
callbackStr='close(gcf)';
closeHndl=uicontrol( ...
'Style','pushbutton', ...
'Units','normalized', ...
'Position',[left frmBottom btnWid btnHt], ...
'String',labelStr, ...
'Callback',callbackStr,...
'Tag',labelStr);
ud.hndlList=[zplaneHndl responseHndl FsHndl WminHndl WmaxHndl pointsHndl ...
radius1Hndl radius2Hndl btn1Hndl btn2Hndl helpHndl closeHndl];
set(figNumber, ...
'Visible','on', ...
'UserData',ud);
design
set(0,'ShowHiddenHandles',shh);
function Fs
ud = get(gcf,'UserData');
hWmax = findobj(ud.hndlList,'Tag','Fmax');
Wmax = get(hWmax,'UserData');
FsHndl = findobj(ud.hndlList,'Tag','Fs');
v = get(FsHndl,'UserData');
s = get(FsHndl,'String');
vv = eval(s,num2str(v));
if vv<Wmax*2, vv = v; end
vv = round(vv*100)/100;
set(FsHndl,'UserData',vv,'String',num2str(double(vv)))
design
function Wmin
ud = get(gcf,'UserData');
hWmax = findobj(ud.hndlList,'Tag','Fmax');
Wmax = get(hWmax,'UserData');
WminHndl = findobj(ud.hndlList,'Tag','Fmin');
v = get(WminHndl,'UserData');
s = get(WminHndl,'String');
vv = eval(s,num2str(v));
if vv>=Wmax, vv = v; end
vv = round(vv*100)/100;
set(WminHndl,'UserData',vv,'String',num2str(vv))
design
function Wmax
ud = get(gcf,'UserData');
hWmin = findobj(ud.hndlList,'Tag','Fmin');
Wmin = get(hWmin,'UserData');
WmaxHndl = findobj(ud.hndlList,'Tag','Fmax');
v = get(WmaxHndl,'UserData');
s = get(WmaxHndl,'String');
vv = eval(s,num2str(v));
if vv<=Wmin, vv = v; end
vv = round(vv*100)/100;
set(WmaxHndl,'UserData',vv,'String',num2str(vv))
design
function points
ud = get(gcf,'UserData');
pointsHndl = findobj(ud.hndlList,'Tag','points');
v = get(pointsHndl,'UserData');
s = get(pointsHndl,'String');
vv = eval(s,num2str(v));
if vv<length(ud.filter.den), vv = v; end
vv = round(vv*100)/100;
set(pointsHndl,'UserData',vv,'String',num2str(vv))
design
function radius1
ud = get(gcf,'UserData');
radius1Hndl = findobj(ud.hndlList,'Tag','Radius1');
v = get(radius1Hndl,'UserData');
s = get(radius1Hndl,'String');
vv = eval(s,num2str(v));
if vv<=0 | vv>10, vv = v; end
vv = round(vv*1000)/1000;
set(radius1Hndl,'UserData',vv,'String',num2str(vv))
design
function radius2
ud = get(gcf,'UserData');
radius2Hndl = findobj(ud.hndlList,'Tag','Radius2');
v = get(radius2Hndl,'UserData');
s = get(radius2Hndl,'String');
vv = eval(s,num2str(v));
if vv<=0 | vv > 10, vv = v; end
vv = round(vv*1000)/1000;
set(radius2Hndl,'UserData',vv,'String',num2str(vv))
design
function radio(s)
hFig=gcf;
ud=get(hFig,'UserData');
axHndl=gca;
rbtnHndl=findobj(ud.hndlList,'Style','radiobutton');
set(rbtnHndl,'Value',0) % Disable all the buttons
set(rbtnHndl(s),'Value',1) % Enable selected button
radiusHndl(1)=findobj(ud.hndlList,'Style','edit','Tag','Radius1');
radiusHndl(2)=findobj(ud.hndlList,'Style','edit','Tag','Radius2');
if strcmp('CZT',get(rbtnHndl(s),'String')),
set(radiusHndl,'Enable','on','BackgroundColor','white');
else
parentHndl=get(radiusHndl(1),'Parent');
bgcolor=get(parentHndl,'Color');
set(radiusHndl,'Enable','off','BackgroundColor',bgcolor);
end
drawnow
design
function infobtn
set(gcf,'Pointer','arrow')
ttlStr = get(gcf,'Name');
hlpStr= [...
'This example shows how to compare the transform of a filter using two '
' different algorithms. They are '
' FFT - the Fast Fourier Transform '
' CZT - the Chirp-Z Transform '
' '
'The upper plot shows the z-plane, and the lower plot shows the transform of'
'a bandpass elliptic digital filter on the points within the wedge shown on '
'the upper plot. The filter has a passband from .4 to .7 of the Nyquist '
'frequency (Nyquist = Fs/2). '
' '
'The FFT computes the Z-transform on equally spaced points around the unit '
'circle. '
' '
'The CZT computes the Z-transform on a spiral or "chirp" contour. The '
'contour is defined by initial frequency Fmin and radius R1, and final '
'frequency Fmax and radius R2. Fs is the sampling frequency. '
' '
'Fmin and Fmax define a "wedge" of the unit circle. '
' '
'Npoints is the number of Z-transform points computed on the unit circle in '
'the wedge defined by Fmin and Fmax. '
' '
'With FFT, the length of the transform is Npoints*Fs/(Fmax-Fmin), which '
'computes Npoints points in the range Fmin to Fmax. If you are interested in'
'a small frequency range, the CZT is more efficient because it "zooms in" on'
'the range you are interested in and only computes Npoints. '
];
myFig = gcf;
helpwin(hlpStr,ttlStr);
% [EOF]