featurexpofn.m

非常好的数字处理教程

function featurexpofn(action,datastruct)
    if nargin < 1
        action='init';
    end


	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.specplot);
			reset(handles.featureplot1);
			reset(handles.featureplot2);
		case 'loadsound'
			audiodata = load_audiodata;
            if isempty(audiodata)
                return
            end
			handles.fileopen = true;
			handles.channels = size(audiodata.data,2);
			if size(audiodata.data,2) > 1,
				audiodata.data = to_mono(audiodata.data);
			end
			handles.audiodata = audiodata;
			handles.t = [0:1/handles.audiodata.Fs:(length(handles.audiodata.data)-1)/handles.audiodata.Fs];
            contents = get(handles.fftsize,'String');
            fftsize = str2double(contents{get(handles.fftsize,'Value')});
            contents = get(handles.window,'String');
            shape = contents{get(handles.window,'Value')};
            [handles.spectrum, handles.bin, handles.st] = ...
                spectrogram(handles.audiodata, fftsize, shape);
            set(handles.freqzoom,'Value',1.0);
			handles = makeSpecPlot(handles);
			handles.featuretype = {{'zero'}, {'rms'}, {'spec_centroid'}, ...
					{'spec_rolloff'}, {'harmonicity'}, {'pitch'}};
			handles.featuredata = {};
			handles = makeFeatureData(handles);
			handles = drawFeaturePlot(handles,1);
			handles = drawFeaturePlot(handles,2);
			timeticks = max(handles.ftimes);
			set(handles.specplot,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',[0 timeticks], 'XTickLabel', '');
			set(handles.featureplot1,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',[0 timeticks], 'XTickLabel', '');
			set(handles.featureplot2,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',[0 timeticks]);
            axes(handles.featureplot2);
            xlabel('Time (s)');
			linkedzoom([handles.specplot, handles.featureplot1, handles.featureplot2],'onx');
			place_header(f,handles);
            featurexpofn 'freqzoom';
		case 'readinput'
			% This will read in analysis data and load accompanying sound
            setdefaults;
			reset(handles.specplot);
			reset(handles.featureplot1);
			reset(handles.featureplot2);
			audiodata = datastruct;
			clear datastruct;
            handles.fileopen = true;
			if size(audiodata.data,2) > 1,
				audiodata.data = to_mono(audiodata.data);
			end
			handles.audiodata = audiodata;
			handles.t = [0:1/handles.audiodata.Fs:(length(handles.audiodata.data)-1)/handles.audiodata.Fs];
            contents = get(handles.fftsize,'String');
            fftsize = str2double(contents{get(handles.fftsize,'Value')});
            contents = get(handles.window,'String');
            shape = contents{get(handles.window,'Value')};
            [handles.spectrum, handles.bin, handles.st] = ...
                spectrogram(handles.audiodata, fftsize, shape);
            set(handles.freqzoom,'Value',1.0);
			handles = makeSpecPlot(handles);
			handles.featuretype = {{'zero'}, {'rms'}, {'spec_centroid'}, ...
					{'spec_rolloff'}, {'harmonicity'}, {'pitch'}};
			handles.featuredata = {};
			handles = makeFeatureData(handles);
			handles = drawFeaturePlot(handles,1);
			handles = drawFeaturePlot(handles,2);
			timeticks = max(handles.ftimes);
			set(handles.specplot,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',[0 timeticks], 'XTickLabel', '');
			set(handles.featureplot1,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',[0 timeticks], 'XTickLabel', '');
			set(handles.featureplot2,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',[0 timeticks]);
            axes(handles.featureplot2);
            xlabel('Time (s)');
			linkedzoom([handles.specplot, handles.featureplot1, handles.featureplot2],'onx');
			place_header(f,handles);
            featurexpofn 'freqzoom';
		case 'playsound'
			if isfield(handles,'audiodata'),
                times = get(handles.featureplot1,'XLim');
                audiodata = handles.audiodata;
                samples = floor(times*audiodata.Fs);
                if (samples(1) <= 0)
                    samples(1) = 1;
                end
                if (samples(2) > length(audiodata.data))
                    samples(2) = length(audiodata.data);
                end
                audiodata.data = audiodata.data(samples(1):samples(2));
                if (max(abs(audiodata.data)) > 1.0)
                    audiodata.data = normalize(audiodata.data);
                end
                play_audiodata(audiodata, handles.play);
            end
		case {'db','inverse','interpolate','colormap'}
			if isfield(handles,'audiodata'),
                xlim = get(handles.specplot,'XLim');
				updateSpecPlot(handles);
                set(handles.specplot,'XLim',xlim);
                featurexpofn 'freqzoom';
			end
		case {'fftsize','window'}
			if isfield(handles,'audiodata'),
				contents = get(handles.fftsize,'String');
                fftsize = str2double(contents{get(handles.fftsize,'Value')});
                contents = get(handles.window,'String');
                shape = contents{get(handles.window,'Value')};

                [handles.spectrum, handles.bin, handles.st] = ...
                    spectrogram(handles.audiodata, fftsize, shape);
                xlim = get(handles.specplot,'XLim');
				updateSpecPlot(handles);
                featurexpofn 'freqzoom';
				handles.featuredata = {};
				handles = makeFeatureData(handles);
				handles = drawFeaturePlot(handles,1);
                handles = drawFeaturePlot(handles,2);
				timeticks = max(handles.ftimes);
				set(handles.specplot,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',xlim, 'XTickLabel', '');
			    set(handles.featureplot1,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',xlim, 'XTickLabel', '');
			    set(handles.featureplot2,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',xlim);
                axes(handles.featureplot2);
                xlabel('Time (s)');
			end
		case {'feature1_update','feature2_update'}
			switch action
				case 'feature1_update'
                    plotnum = 1;
				case 'feature2_update'
                    plotnum = 2;
			end
			handles = drawFeaturePlot(handles, plotnum);
			timeticks = max(handles.ftimes);
            xlim = get(handles.specplot,'XLim');
			set(handles.specplot,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',xlim, 'XTickLabel', '');
			set(handles.featureplot1,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',xlim, 'XTickLabel', '');
			set(handles.featureplot2,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',xlim);
            axes(handles.featureplot2);
            xlabel('Time (s)');
        case 'freqzoom'
			if isfield(handles,'audiodata')
				val = get(handles.freqzoom,'Value');
				if val == 0,
					val = 0.001;
				end
				Fs = handles.audiodata.Fs;
				set(handles.specplot,'YLim',[0 val*Fs/2])
			end
        case 'zoomreset'
            if isfield(handles,'audiodata')
                timeticks = max(handles.ftimes);
                set(handles.specplot,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',[0 max(timeticks)], 'XTickLabel', '');
			    set(handles.featureplot1,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',[0 max(timeticks)], 'XTickLabel', '');
			    set(handles.featureplot2,'XTick',[0:timeticks*0.1:timeticks], ...
					'XLim',[0 max(timeticks)]);
            end
        case {'fourier', 'alias', 'sonogram', 'firfilter', 'pzfilter','sws'}
            if isfield(handles,'audiodata'),
                audiodata = handles.audiodata;
			    switch action
				    case 'fourier'
				        fourierexpogui(audiodata);
                    case 'alias'
                        aliasexpogui(audiodata);
                    case 'sonogram'
                        sonoexpogui(audiodata);
                    case 'firfilter'
                        firexpogui(audiodata);
                    case 'pzfilter'
                        pzfilterexpogui(audiodata);
                    case 'sws'
                        swsexpogui(audiodata);
			    end
            end
		case 'print'
			print_figure(f);
		case 'close'
			close_figure(f,figname(1:end-4));
            return;
	end
	set(f,'UserData',handles);

% --------------------------------------------------------------------
function handles = drawFeaturePlot(handles, plotnum);
	switch plotnum
		case 1
			contents = get(handles.feature1,'String');
			feature = contents{get(handles.feature1,'Value')};
			axes(handles.featureplot1);
		case 2
			contents = get(handles.feature2,'String');
			feature = contents{get(handles.feature2,'Value')};
			axes(handles.featureplot2);
	end

    contents = get(handles.fftsize,'String');
    window_size = str2double(contents{get(handles.fftsize,'Value')});
	contents = get(handles.window,'String');
    shape = contents{get(handles.window,'Value')};
    windowdata = get_windowdata(window_size,shape);
	num_samples = size(handles.audiodata.data,1);
	% What if window_skip used in outside analysis isn't cola?
	window_skip = window_size/2; 	% COLA
    num_windows = floor(num_samples/window_skip)-1;

	handles.ftimes = (1:length(handles.featuredata))*window_skip/handles.audiodata.Fs;
	switch (lower(feature))
		case 'zero crossings'
			fstring = 'zero';
			ylabel('Zero Crossings');
		case 'rms'
			fstring = 'rms';
			ylabel('RMS');
		case 'spectral centroid'
			fstring = 'spec_centroid';
			ylabel('Spectral Centroid (Hz)');
		case 'spectral rolloff'
			fstring = 'spec_rolloff';
			ylabel('Spectral Rolloff');
		case 'harmonicity'
			fstring = 'harmonicity';
			ylabel('Harmonicity');
		case 'pitch'
			fstring = 'pitch';
			ylabel('Pitch (Hz)');
		case 'time series'
			ylabel('Amplitude');
			plot(handles.t,handles.audiodata.data);
			set(gca,'XLim',[0 max(handles.t)]);
			set(gca,'YLim',[-1.1 1.1]);
			grid;
			return;
	end

	col = find(strcmp([handles.featuretype{:}],fstring));
	data = [handles.featuredata{:,col}];
	stairs(handles.ftimes,data);
	grid;

function handles = makeSpecPlot(handles);
	axes(handles.specplot);
	handles.pos = get(gca,'Position'); % Save axes position
	if (get(handles.dB,'Value'))
		imagesc(handles.t,handles.bin,20*log10(abs(handles.spectrum)));
	else
		imagesc(handles.t,handles.bin,abs(handles.spectrum));
	end
	%shading interp;

	axis xy;
	ylabel('Frequency (Hz)');
	colormap('Jet');
	set(handles.specplot,'XTickLabel',['']);
	h = colorbar('vert');
	if (get(handles.dB,'Value'))
		set(get(h, 'YLabel'), 'String', 'dB');
	else
		set(get(h, 'YLabel'), 'String', 'Amplitude');
	end
	set(h,'Position',[0.7 0.6162 0.0404327 0.3651])
	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


function updateSpecPlot(handles)
	axes(handles.specplot);
	% Check Interpolate and dB
    if (get(handles.interpolate,'Value'))
		if (get(handles.dB,'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.dB,'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

	h = colorbar('vert');
	if (get(handles.dB,'Value'))
		set(get(h, 'YLabel'), 'String', 'dB');
	else
		set(get(h, 'YLabel'), 'String', 'Amplitude');
	end
	set(h,'Position',[0.7 0.6162 0.0404327 0.3651])
    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


function handles = makeFeatureData(handles);
    contents = get(handles.fftsize,'String');
    window_size = str2double(contents{get(handles.fftsize,'Value')});
	contents = get(handles.window,'String');
    shape = contents{get(handles.window,'Value')};
    windowdata = get_windowdata(window_size,shape);
	num_samples = size(handles.audiodata.data,1);
	window_skip = window_size/2; 	% COLA
    num_windows = floor(num_samples/window_skip)-1;

	Fs = handles.audiodata.Fs;
	[b,a] = butter(10, 900/(Fs/2)); 
	% Calculate all datatypes if necessary
	if (isempty(handles.featuredata) | handles.feature_window_size ~= window_size)
		h = waitbar(0,'Analyzing');
		for i = 1:num_windows,
			start_sample = (i-1)*window_skip + 1;
			end_sample = start_sample + window_size-1;
			frame = handles.audiodata.data(start_sample:end_sample).*windowdata;
			for j = 1:length(handles.featuretype),
				switch (lower(char(handles.featuretype{j})))
					case 'zero'
						handles.featuredata{i,j} = zero_crossings(frame);
					case 'rms'
						handles.featuredata{i,j} = rms(frame);
					case 'spec_centroid'
						handles.featuredata{i,j} = spec_centroid(frame,Fs);
					case 'spec_rolloff'
						handles.featuredata{i,j} = spec_rolloff(frame,Fs);
					case 'harmonicity'
						handles.featuredata{i,j} = 0;
					case 'pitch'
						% low-pass filter frame < 900 Hz
						frame = filter(b,a,frame);
						handles.featuredata{i,j} = pitch1(frame, Fs);
				end
			end
			waitbar(i/num_windows,h);
		end
		close(h);
		handles.feature_window_size = window_size;
	end