aliasexpofn.m

非常好的数字处理教程

% The following taken in large part from MAD.

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

% clean up freq labels for log freq and linear freq
% Nyquist legend disappers after 'dB'
% 'Zoom to window' button?

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.timeplot);
	case 'loadsound'
		handles.audiodata = load_audiodata;
		if ~isfield(handles.audiodata, 'filenamepath')
			return;
		end
		set(handles.downsampleFactor,'Value',1);
		set(handles.effectiveFs,'String',num2str(handles.audiodata.Fs));
		if (size(handles.audiodata.data,2) > 1)
			handles.audiodata.data = to_mono(handles.audiodata.data);
		end
		handles = makeTimePlot(handles);
		handles = makeSpecPlot(handles);
		place_header(f, handles);
		set(handles.freqzoom,'Value',1);
		set(handles.ampzoom,'Value',0.5);
	case 'readinput'
		handles.audiodata = datastruct;
		clear datastruct;
		handles.audiodata.filenamepath = '';
		handles.audiodata.nBits = 16;
		handles.audiodata.channels = size(handles.audiodata.data,2);
		handles.original = handles.audiodata;
		set(handles.downsampleFactor,'Value',1);
		set(handles.effectiveFs,'String',num2str(handles.audiodata.Fs));
		handles = makeTimePlot(handles);
		handles = makeSpecPlot(handles);
		place_header(f, handles);
	case 'playsound'
		if isfield(handles, 'audiodata')
			audiodata = handles.audiodata;
			if isfield(handles,'filtercoefs')
				a = handles.filtercoefs.a;
				b = handles.filtercoefs.b;
				audiodata.data = filter(b,a,audiodata.data);
			end
			if (max(abs(audiodata.data)) > 1.0)
				audiodata.data = normalize(audiodata.data);
			end
			play_audiodata(audiodata, handles.play);
		end
	case 'normalize'
		if isfield(handles, 'audiodata')
			handles.audiodata.data = normalize(handles.audiodata.data);
			handles = makeTimePlot(handles);
			handles = makeSpecPlot(handles);
			aliasexpofn 'freqzoom';
			aliasexpofn 'ampzoom';
		end
	case 'logfreq'
		if (get(handles.logfreq,'Value'))
			set(handles.specplot,'XScale','log');
		else
			set(handles.specplot,'XScale','linear');
		end
	case {'db'}
		if isfield(handles, 'audiodata')
			updateSpecPlot(handles);
			if (get(handles.dB,'Value'))
				axes(handles.specplot)
				ylabel('Log Magnitude (dB)');
			else
				axes(handles.specplot)
				ylabel('Amplitude');
			end
			aliasexpofn 'freqzoom';
			aliasexpofn 'ampzoom';
		end
	case 'entire'
		if isfield(handles, 'audiodata')
			handles = makeSpecPlot(handles);
			aliasexpofn 'freqzoom';
			aliasexpofn 'ampzoom';
		end
	case 'stem'
		if isfield(handles, 'audiodata')
			handles = makeSpecPlot(handles);
			aliasexpofn 'freqzoom';
			aliasexpofn 'ampzoom';
		end
	case {'fftsize','window','smoothfft'}
		if isfield(handles, 'audiodata')
			contents = get(handles.fftsize,'String');
			fftsize = contents{get(handles.fftsize,'Value')};
			if strncmp(fftsize,'Entire',6)
				fftsize = length(handles.audiodata.data);
			else
				fftsize = str2num(fftsize);
			end
			wintime = fftsize/handles.audiodata.Fs;
			handles.right = handles.left + wintime;
			if handles.right > handles.maxtime
				handles.right = handles.maxtime;
				handles.left = handles.maxtime - wintime;
			end
			set(handles.l,'XData',[handles.left handles.left]);
			set(handles.r,'XData',[handles.right handles.right]);
			handles = makeSpecPlot(handles);
			aliasexpofn 'freqzoom';
			aliasexpofn 'ampzoom';
		end
	case 'left'
		set(gcf,'WindowButtonMotionFcn','aliasexpofn moveleft');
		set(gcf,'WindowButtonUpFcn','aliasexpofn release');
	case 'right'
		set(gcf,'WindowButtonMotionFcn','aliasexpofn moveright');
		set(gcf,'WindowButtonUpFcn','aliasexpofn release');
	case 'moveleft'
		handles = moveleft(handles);
		updateSpecPlot(handles);
	case 'moveright'
		handles = moveright(handles);
		updateSpecPlot(handles);
	case 'release'
		set(gcf,'WindowButtonMotionFcn','');
		set(gcf,'WindowButtonUpFcn','');
		handles = makeSpecPlot(handles);
	case 'downsample'
		contents = get(handles.downsampleFactor,'String');
		factorn = str2num(contents{get(handles.downsampleFactor,'Value')});
		if isfield(handles.audiodata,'data')
			filename = handles.audiodata.filenamepath;
			if ~isempty(filename)
				handles.audiodata = load_audiodata(filename);
			else
				handles.audiodata = handles.original;
			end
			if (size(handles.audiodata.data,2) > 1)
				handles.audiodata.data = to_mono(handles.audiodata.data);
			end
			if get(handles.lowpass,'Value') & factorn > 1
				nyq = 1/factorn;
				[b,a] = ellip(10,0.5,60,nyq*0.97);
				handles.audiodata.data = filter(b,a,handles.audiodata.data);
				handles.filtercoefs.a = a;
				handles.filtercoefs.b = b;
			end
			lengthorig = length(handles.audiodata.data);
			% decimate
			audiodatatemp = handles.audiodata.data(1:factorn:end);
			
			% Upsample to original sampling rate
% 			handles.audiodata.data = upsample(audiodatatemp, factorn);
			
			% Zeroth-order hold to go back to original number of samples
			for i=1:size(audiodatatemp,1)
				handles.audiodata.data((i-1)*factorn+1:i*factorn) = audiodatatemp(i);
			end
			handles.audiodata.data = handles.audiodata.data(1:lengthorig);
			
			set(handles.effectiveFs,'String', ...
				num2str(floor(handles.audiodata.Fs/factorn)));
			updateTimePlot(handles);
			handles = makeSpecPlot(handles);
			aliasexpofn 'freqzoom';
			aliasexpofn 'ampzoom';
		end
	case 'freqzoom'
		if isfield(handles,'audiodata')
			val = get(handles.freqzoom,'Value');
			if val == 0,
				val = 0.001;
			end
			Fs = handles.audiodata.Fs;
			if get(handles.entire,'Value')
				set(handles.specplot,'XLim',[-val*Fs/2 val*Fs/2])
			else
				set(handles.specplot,'XLim',[0 val*Fs/2])
			end
		end
	case 'ampzoom'
		if isfield(handles,'audiodata')
			val = get(handles.ampzoom,'Value');
			if val == 0,
				val = 0.001;
			end
			if get(handles.dB,'Value'),
				set(handles.specplot,'YLim',[-100 (2*val*100-40)])
			else
				set(handles.specplot,'YLim',[0 val*200])
			end
		end
	case 'sonogram'
		if isfield(handles,'audiodata'),
			audiodata = handles.audiodata;
			switch action
				case 'sonogram'
					sonoexpogui(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 = makeTimePlot(handles);
	contents = get(handles.fftsize,'String');
	fftsize = contents{get(handles.fftsize,'Value')};
	if strncmp(fftsize,'Entire',6)
		fftsize = length(handles.audiodata.data);
	else
		fftsize = str2num(fftsize);
	end

	axes(handles.timeplot);
	cla;
	handles.t = [0:1/handles.audiodata.Fs:(length(handles.audiodata.data)-1)/...
		handles.audiodata.Fs];
	maxtime = length(handles.t)/handles.audiodata.Fs;
	handles.maxtime = maxtime;
	handles.winleft = 0;
	handles.winright = maxtime;
	handles.timedata = plot(handles.t,handles.audiodata.data);
	xlabel('time (s)');
	set(handles.timeplot,'NextPlot','replacechildren','Drawmode','fast',...
		'XLim',[0 maxtime], 'YLim', [-1.1 1.1]);
	grid on;
	hold on;
	handles.left = 0.1*maxtime;
	handles.delta = fftsize/handles.audiodata.Fs;
	handles.right = handles.left + handles.delta;
	handles.l=line([handles.left handles.left],get(gca,'YLim'),'Color','r',...
		'ButtonDownFcn','aliasexpofn left','LineWidth', 2,'EraseMode','xor');
	handles.r=line([handles.right handles.right],get(gca,'YLim'),'Color','r',...
		'ButtonDownFcn','aliasexpofn right','LineWidth', 2,'EraseMode','xor');
	handles.step=0.2*maxtime;


function handles = makeSpecPlot(handles);
contents = get(handles.fftsize,'String');
fftsize = contents{get(handles.fftsize,'Value')};
if strncmp(fftsize,'Entire',6)
	fftsize = length(handles.audiodata.data);
else
	fftsize = str2num(fftsize);
end
contents = get(handles.Window,'String');
shape = contents{get(handles.Window,'Value')};

axes(handles.specplot);
cla;
% Get x-position of samples
startpos = floor(handles.left*handles.audiodata.Fs+1);
if (startpos+fftsize) > length(handles.audiodata.data)
	startpos = length(handles.audiodata.data)-fftsize;
end
if startpos <= 0
	handles.audiodata.data(end:end+abs(startpos)) = 0;
	startpos = 1;
end
signal = handles.audiodata.data(startpos:startpos+fftsize-1).*...
	get_windowdata(fftsize, shape);
fftsize = max(2048,fftsize);
f = fftshift(fft(signal,fftsize));
if (get(handles.dB,'Value'))
	f = 20*log10(1e-10+abs(f));
else
	f = abs(f);
end
freqs = linspace(-handles.audiodata.Fs/2,handles.audiodata.Fs/2,fftsize);

if (get(handles.smoothfft,'Value'))
	f=interp1(linspace(-handles.audiodata.Fs/2,handles.audiodata.Fs/2, ...
		length(f)),f,freqs);
end

if get(handles.stem,'Value')
	handles.sp = stem(freqs,f);
else
	handles.sp = plot(freqs,f);
end

contents = get(handles.downsampleFactor,'String');
factorn = str2num(contents{get(handles.downsampleFactor,'Value')});
if get(handles.entire,'Value')
	set(handles.specplot, ...
		'XLim',[-handles.audiodata.Fs/2 handles.audiodata.Fs/2], ...
		'YLim',[-100 60]);
	handles.nyquist = line([-(handles.audiodata.Fs/factorn)/2 ...
		-(handles.audiodata.Fs/factorn)/2], [-100 200], ...
		'Color','r','LineWidth', 2,'EraseMode','xor');
	handles.nyquist = line([(handles.audiodata.Fs/factorn)/2 ...
		(handles.audiodata.Fs/factorn)/2], [-100 200], ...
		'Color','r','LineWidth', 2,'EraseMode','xor');
else
	set(handles.specplot,'XLim',[0 handles.audiodata.Fs/2], ...
		'YLim',[-100 60]);
	handles.nyquist = line([(handles.audiodata.Fs/factorn)/2 ...
		(handles.audiodata.Fs/factorn)/2], ...
		[-100 200],'Color','r','LineWidth', 2,'EraseMode','xor');
end
if (get(handles.logfreq,'Value'))
	set(handles.specplot,'XScale','log');
else
	set(handles.specplot,'XScale','linear');
end
axis xy; grid on;
xlabel('Frequency (kHz)');
ylabel('Log Magnitude (dB)');
legend([handles.nyquist],'Nyquist Frequency');
set(handles.sp, 'EraseMode','background','LineWidth',1);


% --------------------------------------------------------------------
function updateTimePlot(handles)
	axes(handles.timeplot)
	set(handles.timedata,'YData',handles.audiodata.data);

function updateSpecPlot(handles)
contents = get(handles.fftsize,'String');
fftsize = contents{get(handles.fftsize,'Value')};
if strncmp(fftsize,'Entire',6)
	fftsize = length(handles.audiodata.data);
else
	fftsize = str2num(fftsize);
end
contents = get(handles.Window,'String');
shape = contents{get(handles.Window,'Value')};
% Get x-position of samples
startpos = floor(handles.left*handles.audiodata.Fs+1);
signal = handles.audiodata.data(startpos:startpos+fftsize-1).*...
	get_windowdata(fftsize, shape);
fftsize = max(2048,fftsize);
f = fftshift(fft(signal,fftsize));
if (get(handles.dB,'Value'))
	f = 20*log10(1e-10+abs(f));
else
	f = abs(f);
end
freqs = linspace(-handles.audiodata.Fs/2,handles.audiodata.Fs/2,fftsize);
if (get(handles.smoothfft,'Value'))
	f=interp1(linspace(-handles.audiodata.Fs/2,handles.audiodata.Fs/2,length(f)),f,freqs);
end
if (get(handles.logfreq,'Value'))
	set(handles.specplot,'XScale','log');
else
	set(handles.specplot,'XScale','linear');
end
if get(handles.stem,'Value')
	set(handles.sp(1),'YData',f);
	temp = get(handles.sp(2),'YData');
	temp(2:3:end) = f;
	set(handles.sp(2),'YData',temp);
else
	set(handles.sp, 'YData', f);
end


% --------------------------------------------------------------------
function ud = moveleft(ud)
contents = get(ud.fftsize,'String');
fftsize = contents{get(ud.fftsize,'Value')};
if strncmp(fftsize,'Entire',6)
	fftsize = length(ud.audiodata.data);
else
	fftsize = str2num(fftsize);
end
currentPoint=get(gca,'CurrentPoint');
ud.left=currentPoint(1);
wintime = fftsize/ud.audiodata.Fs;

if ud.left <= ud.winleft
	ud.left=0;
end
if ud.left >= (ud.winright-wintime)
	ud.left=ud.winright-wintime;
end

set(ud.l,'XData',[ud.left ud.left]);
set(ud.r,'XData',[ud.left+wintime ud.left+wintime]);

function ud = moveright(ud)
contents = get(ud.fftsize,'String');
fftsize = contents{get(ud.fftsize,'Value')};
if strncmp(fftsize,'Entire',6)
	fftsize = length(ud.audiodata.data);
else
	fftsize = str2num(fftsize);
end
currentPoint=get(gca,'CurrentPoint');
ud.right=currentPoint(1);
wintime = fftsize/ud.audiodata.Fs;

if ud.right <= (ud.winleft+wintime)
	ud.right=wintime;
end
if ud.right >= (ud.winright)
	ud.right=ud.winright;
end
ud.left = ud.right - wintime;

set(ud.r,'XData',[ud.right ud.right]);
set(ud.l,'XData',[ud.right-wintime ud.right-wintime]);