iirexpofn.m

非常好的数字处理教程

	handles.ylim_specplot_orig = get(handles.specplot,'YLim');
	handles.xlim_timeplot_orig = get(handles.timeplot,'XLim');
	handles.ylim_timeplot_orig = get(handles.timeplot,'YLim');
	handles.xlim_specplot = handles.xlim_specplot_orig;
	handles.ylim_specplot = handles.ylim_specplot_orig;
	handles.xlim_timeplot = handles.xlim_timeplot_orig;
	handles.ylim_timeplot = handles.ylim_timeplot_orig;
	handles.clim = get(handles.specplot,'CLim');


% --------------------------------------------------------------------
function updateTimePlot(handles)
	axes(handles.timeplot)
	plot(handles.t,handles.audiodata)
	set(handles.timeplot,'XLim',handles.xlim_timeplot);
	set(handles.timeplot,'YLim',handles.ylim_timeplot);
	xlabel('time (s)');
	grid;

function updateSpecPlot(handles)
	axes(handles.specplot);
	plotype = get(handles.fftdisplay,'String');
	switch lower(plotype{get(handles.fftdisplay,'Value')})
		case 'sonogram'
			% Check Interpolate and dB
			if (get(handles.interpolate,'Value'))
				if (get(handles.dBcheckbox,'Value'))
					pcolor(handles.st,handles.bin,20*log10(abs(handles.spectrum)));
				else
					pcolor(handles.st,handles.bin,abs(handles.spectrum));
				end
				shading interp;
			else
				if (get(handles.dBcheckbox,'Value'))
					imagesc(handles.t,handles.bin,20*log10(abs(handles.spectrum)));
				else
					imagesc(handles.t,handles.bin,abs(handles.spectrum));
				end
			end
			set(handles.specplot,'XTickLabel',['']);
			axis xy;
			ylabel('Frequency (Hz)');
			% Get Colormap
			contents = get(handles.colormap,'String');
			cmap = colormap(lower(contents{get(handles.colormap,'Value')}));
			% Handle Inverse
			if (get(handles.inverse,'Value'))
				colormap(flipud(cmap));
			else
				colormap(cmap);
			end
			set(handles.specplot,'XLim',handles.xlim_specplot);
			set(handles.specplot,'YLim',handles.ylim_specplot);
			set(handles.specplot,'CLim',handles.clim);
			h = colorbar('vert');
			if (get(handles.dBcheckbox,'Value'))
				set(get(h, 'YLabel'), 'String', 'dB');
			else
				set(get(h, 'YLabel'), 'String', 'Amplitude');
			end
			set(handles.specplot,'Units','Characters');
	        sonopos = get(handles.specplot,'Position')
	        set(h,'units','Characters', ...
		        'Position',[sonopos(1)+sonopos(3)+18 sonopos(2) 7 sonopos(4)])
            switch computer
                case {'GLNX86','MAC'}
                    %set(gca,'Position',handles.pos + [0 0.07 0 0])
                    set(gca,'Position',handles.pos)
                case 'PCWIN'
                    set(gca,'Position',handles.pos)
            end
			
		case 'waterfall'
			% Waterfall does not work on Linux, Mac, or Windows!
			%waterfall(handles.st,handles.bin,abs(handles.spectrum));
			%if (get(handles.dBcheckbox,'Value'))
				%me = surf(handles.st,handles.bin,20*log10(abs(handles.spectrum)));
                %else
				% perhaps downsample
				%me = surf(handles.st,handles.bin,abs(handles.spectrum));
				% This would be nice but unecessary
				%explore3d(gca);
            %end
			%zoom off;
			%xlabel('Time')
			%ylabel('Frequency')
			%if (get(handles.interpolate,'Value'))
			%	shading interp;
            %end
		end


% --------------------------------------------------------------------
function plotfreqz(handles)
	% Find a way to keep original xlims
	[a_k,b_k] = designfilter(handles);
	%[H, W] = freqz(b_k, a_k, 512, 'whole');
	[H, W] = freqz(b_k, a_k);

	if isfield(handles,'audiodata')
		Fs = handles.audiodata.Fs;
	else
		Fs = 44100;
	end

    f=findobj('Tag','freqzplot');
	if isempty(f)
	    f = figure('units','normal','paperunits','normal',...
	   		 'position',[.1,.1,.6,.6], 'Tag', 'freqzplot','MenuBar','none');
	    uicontrol('String','dB','Style','checkbox','Value', 0, ...
	   		'Position',[5 400 40 20],'Tag','dB',...
	   		'Callback','iirexpofn plotfreqz');
	    uicontrol('String','unwrap','Style','checkbox','Value', 1, ...
	   		'Position',[5 200 80 20],'Tag','unwrap',...
	   		'Callback','iirexpofn plotfreqz');
	end
	dB = findobj(f,'Tag','dB');
	unwrapbox = findobj(f,'Tag','unwrap');

	% Plot frequency response of H(w)
	freqplot = subplot(3,1,1);
	W = W.*Fs/(2*pi);
	if get(dB,'Value')
		plot(W,20*log10(abs(H)) - max(20*log10(abs(H))));
		axis([0.0 Fs/2 -200 10.1]);
		ylabel('Magnitude (dB)');
	else
		plot(W,abs(H)./max(abs(H)));
		axis([0.0 Fs/2 -0.1 1.1]);
		ylabel('Magnitude');
	end
	set(freqplot,'XTickLabel','');
	t_h = title('Frequency, Phase and Impulse Response');
	set(t_h, 'fontsize',15, 'fontweight','bold');
	grid;

	% Plot phase response of H(w)
	phasplot = subplot(3,1,2);
	if get(unwrapbox,'Value')
		plot(W,unwrap(angle(H)));
		axis([0.0 Fs/2 1.1*min(unwrap(angle(H))) 1.1*max(unwrap(angle(H)))]);
	else
		plot(W,angle(H));
		axis([0.0 Fs/2 -pi +pi]);
	end
	xlabel('Frequency (Hz)');
	ylabel('Phase (radians)');
	grid;
	linkedzoom([freqplot,phasplot],'onx');

	% Plot impulse response
    implot = subplot(3,1,3);
    impulse = filter(b_k, a_k, [1, zeros(1,99)]);
    plot(impulse);
    axis([0 length(impulse) -max(abs(impulse))*1.1 max(abs(impulse))*1.1]);
    xlabel('Sample');
    ylabel('Amplitude');
    grid;

function [a,b] = designfilter(handles)
	order = str2double(get(handles.filtorder,'String'));
	contents = get(handles.filtermenu,'String');
	typecontents = get(handles.type,'String');
	cut1 = str2double(get(handles.cut1,'String'));
	cut2 = str2double(get(handles.cut2,'String'));
	Rp = str2num(get(handles.rp,'String'));
	Rs = str2num(get(handles.rs,'String'));
	if isfield(handles,'audiodata')
		Fs = handles.audiodata.Fs;
	else
		Fs = 44100;
	end
	cut1 = cut1/(Fs/2);
	cut2 = cut2/(Fs/2);
	switch lower(typecontents{get(handles.type,'Value')})
		case 'butterworth'
			switch (lower(contents{get(handles.filtermenu,'Value')}))
				case 'allpass'
					b = 1;
					a = 1;
				case 'lowpass'
					[b,a] = butter(order, cut1);
				case 'highpass'
					[b,a] = butter(order, cut1,'high');
				case 'bandpass'
					[b,a] = butter(order, [cut1 cut2]);
				case 'notch'
					[b,a] = butter(order, [cut1 cut2], 'stop');
			end
		case 'elliptical'
			switch (lower(contents{get(handles.filtermenu,'Value')}))
				case 'allpass'
					b = 1;
					a = 1;
				case 'lowpass'
					[b,a] = ellip(order,Rp,Rs, cut1);
				case 'highpass'
					[b,a] = ellip(order,Rp,Rs, cut1,'high');
				case 'bandpass'
					[b,a] = ellip(order,Rp,Rs, [cut1 cut2]);
				case 'notch'
					[b,a] = ellip(order,Rp,Rs, [cut1 cut2], 'stop');
			end
		case 'chebyshev'
			switch (lower(contents{get(handles.filtermenu,'Value')}))
				case 'allpass'
					b = 1;
					a = 1;
				case 'lowpass'
					[b,a] = cheby1(order,Rp, cut1);
				case 'highpass'
					[b,a] = cheby1(order,Rp, cut1,'high');
				case 'bandpass'
					[b,a] = cheby1(order,Rp, [cut1 cut2]);
				case 'notch'
					[b,a] = cheby1(order,Rp,[cut1 cut2], 'stop');
			end
	end

function handles = apply_filter(handles)
	[a_k,b_k] = designfilter(handles);
	handles.audiodata2 = handles.audiodata;
	handles.audiodata.data = filter(b_k, a_k, handles.audiodata.data);
	% Store old data for undo
	handles.spectrum2 = handles.spectrum;
	handles.bin2 = handles.bin;
	handles.st2 = handles.st;
	set(handles.undo,'String','Undo');

function handles = undo(handles)
	hObject = handles.undo;
	if (strcmp(get(hObject,'String'),'Undo'))
		% Store old data for redo
		handles.spectrum3 = handles.spectrum;
		handles.bin3 = handles.bin;
		handles.st3 = handles.st;
		handles.audiodata3 = handles.audiodata;
		% Reload old data
		handles.spectrum = handles.spectrum2;
		handles.bin = handles.bin2;
		handles.st = handles.st2;
		handles.audiodata = handles.audiodata2;
		set(hObject,'String','Redo');
	else
		handles.spectrum = handles.spectrum3;
		handles.bin = handles.bin3;
		handles.st = handles.st3;
		handles.audiodata = handles.audiodata3;
		set(hObject,'String','Undo');
	end