www.gusucode.com > 声音的处理有：LPC，FFT，共振峰，频谱源码程序 > siganlandsystemusingMatlab/SSUM/sample/samplexpofn.m

    function samplexpofn(action)

    if nargin < 1
        action='init';
    end

    % Option to show error?
    
	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;
			samplexpofn 'update';
		case 'update'
			update(handles);
			linkedzoom([handles.sampleplot, handles.interplot],'onx');
	end
	set(f,'UserData',handles);


% --------------------------------------------------------------------
function update(handles)
	Fs = str2double(get(handles.Fs,'String'));
	periods = str2double(get(handles.periods,'String'));
	freq1 = str2double(get(handles.freq1,'String'));
	amp1 = str2double(get(handles.amp1,'String'));
	phase1 = str2double(get(handles.phase1,'String'))*pi/180;
	offset = str2double(get(handles.offset,'String'));

% 	sample_t = [0:1/Fs:periods/freq1];
% 	t = [0:1/44100:periods/freq1];
	dur = 0.2;
	sample_t = [0:1/Fs:dur];
	aFs = 11025;
	dt = 0.05;
	t = [0:1/aFs:dur];
	
	plotdur = periods/freq1;

	contents = get(handles.wave,'String');
	switch lower(contents{get(handles.wave,'Value')})
		case 'sine'
			signal = offset + amp1*sin(2*pi*freq1*t + phase1);
			sample_signal = offset + amp1*sin(2*pi*freq1*sample_t + phase1);
		case 'triangle'
			signal = offset + amp1*sawtooth(2*pi*freq1*t + phase1,0.5);
			sample_signal = offset + amp1*sawtooth(2*pi*freq1*sample_t + phase1,0.5);
		case 'sawtooth'
			signal = offset + amp1*sawtooth(2*pi*freq1*t + phase1);
			sample_signal = offset + amp1*sawtooth(2*pi*freq1*sample_t + phase1);
		case 'square'
			signal = offset + amp1*square(2*pi*freq1*t + phase1);
			sample_signal = offset + amp1*square(2*pi*freq1*sample_t + phase1);
	end
	%Quantize samples
	nbits = str2num(get(handles.nbits,'String'));
	sample_signal = quantize(sample_signal,nbits);

	axes(handles.sampleplot);
	cla;
	plot(t-dt,signal);
	
	ylabel('Input');
	hold on;
	if (get(handles.samplebox,'Value'))
		stem(sample_t-dt,sample_signal,'k');
	end

	if (get(handles.gridbox,'Value'))
		grid on;
	else
		grid off;
	end
	% Set plot properties
	set(handles.sampleplot,'YLim',[-1.1 1.1]);
	set(handles.sampleplot,'XLim',[0 plotdur],'XTickLabel','');

	axes(handles.interplot);
	cla;
	if (get(handles.samplebox,'Value'))
		stem(sample_t-dt,sample_signal,'k');
	end
	hold on;

	% Now interpolate
	sample_times = sample_t';
	ts = linspace(0,dur,aFs*dur)';
	a = ts(:,ones(size(sample_times)));
	b = sample_times(:,ones(size(ts)))';
	y = sinc(Fs*(a-b))*sample_signal';
	h = plot(ts-dt,y);
	set(h,'LineWidth',2.0);

	if (get(handles.interpbox,'Value'))
		sinct = -0.1:1/aFs:0.1;
		sincfn = sinc(sinct*Fs);
		% Plot sincs at samples
		for i=1:length(sample_signal),
			plot(sinct+sample_t(i)-dt,sincfn*sample_signal(i));
		end
	end

	xlabel('Time (s)');
	ylabel('Output');
	if (get(handles.gridbox,'Value'))
		grid on;
	else
		grid off;
	end
	% Set plot properties
	set(handles.interplot,'YLim',[-1.1 1.1]);
	set(handles.interplot,'XLim',[0 plotdur]);


% --------------------------------------------------------------------

function signal = quantize(signal,nbits,maxbits)
	% Signal should be -1 <= x <= 1
	if (nbits == 1)
		ps = find(signal < 0);
		signal(ps) = -1;
		signal = ceil(signal);
	else
		signal = floor(signal.*2^(nbits-1))/2^(nbits-1);
		ps = find(signal == 0);
		signal(ps) = 1/2^(nbits-1);
	end