pzexpofn.m

非常好的数字处理教程

function pzexpofn(action,datastruct)
% This code borrows heavily from MAD polezero demonstration
    if nargin < 1
        action='init';
    end

	% Maybe 'apply to white noise' button
	% Send filter to pzfilterexplorer
	% How do you set the plots so that one zooms and the other doesn't????
	% Phase in radians, +- pi
    % Fix impulse response, zoom in...

	name = mfilename;
	figname = [name(1:end-2) '_fig'];
	f=findobj('Tag',figname);
    handles = get(f,'UserData');

    switch action
		case 'help'
			display_help(figname);
        case 'init'
			setdefaults;
			reset(handles.pz_plot);
			reset(handles.freq_plot);
			reset(handles.phase_plot);
			reset(handles.impulse_plot);
			handles.npoints = 256;
			makePZPlot(handles);
			handles = makePlots(handles);
			set(gcf,'KeyPressFcn','pzexpofn keypressed');
			%linkedzoom([handles.freq_plot, handles.phase_plot],'onx');
		case 'readinput'
			setdefaults;
			reset(handles.pz_plot);
			reset(handles.freq_plot);
			reset(handles.phase_plot);
			reset(handles.impulse_plot);
            handles.Fs = 44100;
			handles.npoints = 256;
			makePZPlot(handles);
			handles = makePlots(handles);
			set(gcf,'KeyPressFcn','pzexpofn keypressed');
            if isfield(datastruct, 'gain')
                set(handles.gain,'String',num2str(datastruct.gain));
            end
			z = datastruct.zeropos; 
			% Remove zeros from 0
			z = z(find(z));
			zc = conj(z);
            elem = find(z == zc);
            if (elem > 0)
                % Add just a bit so that the plot doesn't make them invisible
				mult = 0.00000;
                z = [z(1:elem-1); z(elem)+j*mult; z(elem)-j*mult; z(elem+1:end)];
				datastruct.zeropos = z;
            end
            p = datastruct.polepos;
			% Remove poles from 0
			p = p(find(p));
            pc = conj(p);
            elem = find(p == pc);
            if (elem > 0)
				mult = 0.00000;
                p = [p(1:elem-1); p(elem)+j*mult; p(elem)-j*mult; p(elem+1:end)];
                datastruct.polepos = p;
            end

			if isfield(datastruct,'zeropos')
                for i=1:length(datastruct.zeropos)/2,
                    x = real(datastruct.zeropos(i*2-1));
                    y = imag(datastruct.zeropos(i*2-1));
                    p1=text('Parent',handles.pz_plot,'Position',[x y],...
                        'String','o','ButtonDownFcn','pzexpofn selectzero',...
                        'EraseMode','xor','Tag','Zero',...
                        'HorizontalAlignment','center');
					% Don't create two where there should be one
					if y ~= 0
						p2=text('Parent',handles.pz_plot,'Position',[x -y],...
							'String','o','ButtonDownFcn','pzexpofn selectzero',...
							'EraseMode','xor','Tag','Zero',...
							'HorizontalAlignment','center');
						set(p1,'UserData',p2);
						set(p2,'UserData',p1);
					end
                end
                for i=1:length(datastruct.polepos)/2,
                    x = real(datastruct.polepos(i*2-1));
                    y = imag(datastruct.polepos(i*2-1));
                    p1=text('Parent',handles.pz_plot,'Position',[x y],...
                        'String','x','ButtonDownFcn','pzexpofn selectpole',...
                        'EraseMode','xor','Tag','Pole',...
                        'HorizontalAlignment','center');
					% Don't create two where there should be one
					if y ~= 0
						p2=text('Parent',handles.pz_plot,'Position',[x -y],...
							'String','x','ButtonDownFcn','pzexpofn selectpole',...
							'EraseMode','xor','Tag','Pole',...
							'HorizontalAlignment','center');
						set(p1,'UserData',p2);
						set(p2,'UserData',p1);
					end
                end
			end
			handles = updatePlots(handles);
		case 'db'
			axes(handles.freq_plot);
			if get(handles.dB,'Value')
				set(handles.freq_plot,'YLim',[-80 1]);
    			ylabel('Normalized Attenuation (dB)');
			else
				set(handles.freq_plot,'YLim',[0 1.05]);
    			ylabel('Normalized Attenuation');
			end
			pzexpofn 'plotfreqz'
		case 'unwrap'
			if get(handles.unwrap,'Value')
				set(handles.phase_plot,'YLimMode','auto');
			else
				set(handles.phase_plot,'YLim',[-pi pi]);
			end
			pzexpofn 'plotfreqz'
		case 'grid'
			if get(handles.grid,'Value')
				axes(handles.freq_plot);
				grid on;
				axes(handles.phase_plot);
				grid on;
				axes(handles.impulse_plot);
				grid on;
			else
				axes(handles.freq_plot);
				grid off;
				axes(handles.phase_plot);
				grid off;
				axes(handles.impulse_plot);
				grid off;
			end
		case 'logf'
			if get(handles.logf,'Value')
				set(handles.freq_plot,'XScale','log');
				set(handles.phase_plot,'XScale','log');
			else
				set(handles.freq_plot,'XScale','linear');
				set(handles.phase_plot,'XScale','linear');
			end
		case 'loadvowel'
			axes(handles.pz_plot);
            cla;
            delete(findobj(gcf,'Selected','on'))
            makePZPlot(handles);
            bws=[0.03 0.07 0.1];    %formant freqs from Rabiner & Juang
            switch datastruct
                case 'beet'
                    freqs=[270 2290 3010];
                case 'bit'
                    freqs=[390 1990 2550];
                case 'bet'
                    freqs=[530 1840 2480];
                case 'bat'
                    freqs=[660 1720 2410];
                case 'bart'
                    freqs=[730 1090 2440];
                case 'bort'
                    freqs=[570 840 2410];
                case 'but'
                    freqs=[440 1020 2240];
                case 'boot'
                    freqs=[300 870 2240];
                case 'bert'
                    freqs=[490 1350 1690];
                otherwise
            end 
            % convert freqs and bws to PZ locations
            angs=2*pi*freqs/44100;
            for fr=1:length(freqs)
                [x,y]=pol2cart(angs(fr),1-bws(fr));
                p1=text('Parent',handles.pz_plot,'Position',[x y],...
                    'String','x','ButtonDownFcn','pzexpofn selectpole',...
                    'EraseMode','xor','Tag','Pole','HorizontalAlignment','center');
                p2=text('Parent',handles.pz_plot,'Position',[x -y],...
                    'String','x','ButtonDownFcn','pzexpofn selectpole',...
                    'EraseMode','xor','Tag','Pole','HorizontalAlignment','center');
				set(p1,'UserData',p2);
				set(p2,'UserData',p1);
            end     
			pzexpofn 'plotfreqz'
		case 'addzero'
			handles = add_coef(handles,'zero');
			pzexpofn 'plotfreqz'
		case 'addpole'
			handles = add_coef(handles,'pole');
			pzexpofn 'plotfreqz'
		case 'selectpole'
			handles = select_coef(handles,'pole');
		case 'selectzero'
			handles = select_coef(handles,'zero');
		case 'movepole'
			handles = move_coef(handles);
			pzexpofn 'plotfreqz'
		case 'movezero'
			handles = move_coef(handles);
			pzexpofn 'plotfreqz'
		case 'releaseobject'
			set(gcf,'WindowButtonMotionFcn','');
			set(gcf,'WindowButtonUpFcn','');
			pzexpofn 'plotfreqz'
		case 'plotfreqz'
			handles = updatePlots(handles);
			pzexpofn 'timezoom';
		case 'keypressed'
			% if delete key is pressed, delete any selected object
			if double(get(gcf,'CurrentCharacter')) == 8
				delete(findobj(gcf,'Selected','on'))
			end
			pzexpofn 'plotfreqz'
		case 'convexpo'
			data.impulse = get(handles.implot(1),'YData');
			data.impulse = data.impulse(1:100);
			convexpogui(data);
		case 'pzfilterexpo'
			polepos=findobj(handles.pz_plot,'Tag','Pole');
			zeropos=findobj(handles.pz_plot,'Tag','Zero');
			pp = [];
			zz = [];
			for p=1:length(polepos)
				ppos=get(polepos(p),'Position');
				pp(p)=ppos(1)+ppos(2)*j;
			end
			for z=1:length(zeropos)
				zp=get(zeropos(z),'Position');
				zz(z)=zp(1)+zp(2)*j;
			end
			data.polepos = pp;
			data.zeropos = zz;
            data.gain = str2num(get(handles.gain,'String'));
			pzfilterexpogui(data);
		case 'timezoom'
			if isfield(handles,'imp')
                val = get(handles.timezoom,'Value');
                if val == 0,
                    val = 0.001;
                end
                set(handles.impulse_plot,'XLim',[0 val*100])
            end
		case 'unselectall'
			set(findobj(gcf,'Selected','on'),'Selected','off');
		case 'print'
			print_figure(f);
		case 'close'
			close_figure(f, figname(1:end-4));
			return;
	end
	set(f,'UserData',handles);

% --------------------------------------------------------------------
function makePZPlot(handles);
	axes(handles.pz_plot);
	set(handles.pz_plot, 'DrawMode','fast', 'NextPlot','replacechildren',...
		'Box','on', 'XLim', [-1.3 1.3], 'YLim', [-1.3 1.3]);

	% define a circle
	th = 0:pi/50:2*pi;
	xunit = [0.2:0.2:1.0]'*cos(th);
	yunit = [0.2:0.2:1.0]'*sin(th);
	rmax=1;
	tc = get(handles.pz_plot,'xcolor');

	% now really force points on x/y axes to lie on them exactly
	%inds = 1:(length(th)-1)/4:length(th);
	%xunit(inds(2:2:4)) = zeros(2,1);
	%yunit(inds(1:2:5)) = zeros(3,1);

	patch('xdata',xunit(5,:).*rmax,'ydata',yunit(5,:).*rmax, ...
		'edgecolor',tc,'facecolor',get(gca,'color'),...
		'ButtonDownFcn','pzexpofn unselectall');
	hold on;
	for i=1:4,
		plot(yunit(i,:).*rmax, xunit(i,:).*rmax, 'k:')
	end

	line([-1.1 1.1],[0 0]);
	line([0 0],[-1.1 1.1]);
	axis(rmax*[-1 1 -1.15 1.15]);
	axis image;
	axis off;

function handles = makePlots(handles);
	axes(handles.freq_plot);
	handles.freqs=linspace(0,1,handles.npoints);
	handles.freqline=line(handles.freqs,zeros(size(handles.freqs)));
	set(handles.freqline,'EraseMode','xor');
	axis normal;	% Strange hack.
	set(handles.freq_plot,'XLim',[0 0.5]);
	set(handles.freq_plot, 'DrawMode','fast', 'NextPlot','replacechildren',...
		'Box','off');
	set(handles.freq_plot,'XTick',[],'XTickLabel','');
    ylabel('Normalized Attenuation');
	grid;

	axes(handles.phase_plot);
	handles.phaseline=line(handles.freqs,zeros(size(handles.freqs)));
	set(handles.phaseline,'EraseMode','xor');
	axis normal;	% Strange hack.
	set(handles.phase_plot,'YLim',[-pi pi],'XLim',[0 0.5]);
	set(handles.phase_plot, 'DrawMode','fast', 'NextPlot','replacechildren',...
		'Box','off');
	xlabel('Normalized Frequency ({\varpi})');
    ylabel('Phase (radians)');
	set(handles.freq_plot,'XTick',get(handles.phase_plot,'XTick'));
	grid;

	
	axes(handles.impulse_plot);
	handles.imp = zeros(1,500);
	%handles.implot = plot(handles.imp);
	handles.implot = stem(handles.imp);
	set(handles.implot,'EraseMode','xor');
	set(handles.impulse_plot,'XLim',[-1 100]);
	set(handles.impulse_plot, 'DrawMode','fast', 'NextPlot','replacechildren');
	xlabel('Sample');
	grid;
	%zoom on;
	

% --------------------------------------------------------------------
function ud = updatePlots(ud)
    polepos=findobj(ud.pz_plot,'Tag','Pole');
    zeropos=findobj(ud.pz_plot,'Tag','Zero');
    len=max(length(polepos),length(zeropos));
    pp=zeros(1,len);  % pp=1e-10*rand(1,len);
    for p=1:length(polepos)
        ppos=get(polepos(p),'Position');
        pp(p)=ppos(1)+ppos(2)*i;
    end
    zz=zeros(1,len);  % zz=1e-10*rand(1,len);
    for z=1:length(zeropos)
        zp=get(zeropos(z),'Position');
        zz(z)=zp(1)+zp(2)*i;
    end
    ud.zz=zz;ud.pp=pp;
    [ud.num,ud.den]=zp2tf(zz',pp',1);
    [h,w]=freqz(ud.num,ud.den,ud.npoints,'whole');
    gain = str2num(get(ud.gain,'String'));
	if get(ud.dB,'Value')
		% normalise so max = 0 dB
		logmag=20*log10(abs(h));
		mag=logmag-max(logmag);
	else
		mag = gain*abs(h)./max(abs(h));
        %mag = gain*abs(h);
	end
    set(ud.freqline,'YData',mag);
    if get(ud.dB,'Value')
    else
        set(ud.freq_plot,'YLim',[0 max(mag)*1.1])
    end
	if get(ud.unwrap,'Value')
	    set(ud.phaseline,'YData',unwrap(angle(h)));
		set(ud.phase_plot,'YLim',[min(unwrap(angle(h)))-1 max(unwrap(angle(h)))+1]);
	else
	    set(ud.phaseline,'YData',angle(h));
	end
	imp = filter(ud.num, ud.den, [1, zeros(1,499)]);	
    set(ud.implot,'YData',imp(1:500));
% 	temp = get(ud.implot(2),'YData');
% 	temp(2:3:end) = imp(1:500);
%     set(ud.implot(2),'YData',temp);
	set(ud.impulse_plot,'YLim',[1.1*min(imp) 1.1*max(imp)]);
	


% --------------------------------------------------------------------
function ud = add_coef(ud,type)
	switch type
		case 'zero'
			string = 'o';
			callbak = 'pzexpofn selectzero';
			tag = 'Zero';
		case 'pole'
			string = 'x';
			callbak = 'pzexpofn selectpole';
			tag = 'Pole';
	end
	freq = rand*22050/2;
	ang = 2*pi*freq/44100;
	[x,y]=pol2cart(ang,1.0*rand);
	p1=text('Parent',ud.pz_plot,'Position',[x y],'String',string,...
		'ButtonDownFcn',callbak,'EraseMode','xor',...
		'Tag',tag,'HorizontalAlignment','center');
	p2=text('Parent',ud.pz_plot,'Position',[x -y],'String',string,...
		'ButtonDownFcn',callbak, 'EraseMode','xor',...
		'Tag',tag,'HorizontalAlignment','center');
	set(p1,'UserData',p2);
	set(p2,'UserData',p1);

function ud = select_coef(ud,type)
	switch type
		case 'zero'
			buttonfn = 'pzexpofn movezero';
		case 'pole'
			buttonfn = 'pzexpofn movepole';
	end
	s=get(gcbo,'UserData');
	set(gcf,'WindowButtonMotionFcn', buttonfn);
	set(gcf,'WindowButtonUpFcn','pzexpofn releaseobject');
	ud.currentcoef=gcbo;
	ud.partnercoef=s;
	set([findobj(gca,'Tag','Pole'); findobj(gca,'Tag','Zero')],'Selected','off');
	set([ud.currentcoef ud.partnercoef],'Selected','on');

function ud = move_coef(ud)
	cp=get(ud.pz_plot,'CurrentPoint');
	if (cp(1)*cp(1)+cp(3)*cp(3)) < 2.1
	  set(ud.currentcoef,'Position',[cp(1) cp(3)]);
	  set(ud.partnercoef,'Position',[cp(1) -cp(3)]);
	end

function showImpulse(ud)
    polepos=findobj(ud.pz_plot,'Tag','Pole');
    zeropos=findobj(ud.pz_plot,'Tag','Zero');
    len=max(length(polepos),length(zeropos));
    pp=zeros(1,len);  % pp=1e-10*rand(1,len);
    for p=1:length(polepos)
        ppos=get(polepos(p),'Position');
        pp(p)=ppos(1)+ppos(2)*i;
    end
    zz=zeros(1,len);  % zz=1e-10*rand(1,len);
    for z=1:length(zeropos)
        zp=get(zeropos(z),'Position');
        zz(z)=zp(1)+zp(2)*i;
    end
    [b_k,a_k]=zp2tf(zz',pp',1);