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

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

% Analysis size text entry
% Synthesis window, size, skip parameters
% Option to apply specific amplitude envelope to synthesis
%	such as RMS of original
% Options to substitute particular bins

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.timeplot1);
		reset(handles.timeplot2);
		reset(handles.resplot);
		set(handles.excite_text,'Visible','off');
		set(handles.excite_freq,'Visible','off');
		linkedzoom([handles.timeplot1,handles.timeplot2, handles.resplot],'onxy');
		cola_check(handles);
	case 'loadsound'
		audiodata = load_audiodata;
		if isempty(audiodata)
			return
		end
		if (max(abs(audiodata.data)) > 1.0)
			audiodata.data = normalize(audiodata.data);
		end
		handles.audiodata = audiodata;
		makeTimePlot(handles.audiodata, handles.timeplot1);
		place_header(f,handles);
	case 'readinput'
		handles.audiodata = datastruct;
		clear datastruct;
		handles.audiodata.filenamepath = '';
		handles.audiodata.nBits = 16;
		handles.audiodata.channels = size(handles.audiodata.data,2);
		setdefaults;
		reset(handles.timeplot1);
		reset(handles.timeplot2);
		reset(handles.resplot);
		makeTimePlot(handles.audiodata, handles.timeplot1);
		place_header(f,handles);
	case 'analyze'
		if isfield(handles,'audiodata')
			windowsize = str2double(get(handles.analwindowsize,'String'));
			windowskip = str2double(get(handles.analwindowskip,'String'));
			windowshapecontents = get(handles.analwindowshape,'String');
			windowshape = windowshapecontents{get(handles.analwindowshape,'Value')};
			order = str2num(get(handles.order,'String'));
			[handles.analysisdata, handles.residual] = analyze(handles.audiodata, windowsize, ...
				windowskip, windowshape, order);
			% Get and plot residual
% 			handles.residual = get_residual(handles.analysisdata, ...
% 				handles.audiodata, windowsize, windowskip, windowshape);
			[m,n] = size(handles.audiodata.data);
			handles.residual.data = handles.residual.data(1:m,1:n);
			makeTimePlot(handles.residual, handles.resplot);
		end
	case 'synthesize'
		if isfield(handles,'analysisdata')
			windowsize = str2double(get(handles.analwindowsize,'String'));
			windowskip = str2double(get(handles.analwindowskip,'String'));
			windowshapecontents = get(handles.analwindowshape,'String');
			windowshape = windowshapecontents{get(handles.analwindowshape,'Value')};
			
			excitecontents = get(handles.excite_menu,'String');
			switch lower(excitecontents{get(handles.excite_menu,'Value')})
				case 'residual'
					excitation = handles.residual;
				case 'noise'
					[m,n] = size(handles.residual.data);
					excitation.data = 1-2*rand(m,n);
					excitation.Fs = handles.residual.Fs;
				case 'pulses'
					[m,n] = size(handles.residual.data);
					freq = str2num(get(handles.excite_freq,'String'));
					per = round(handles.residual.Fs/freq);
					excitation.data = zeros(m,n);
					excitation.data(1:per:end,1:n) = 1;
					excitation.Fs = handles.residual.Fs;
				case 'other'
					% Load soundfile
					excitation = load_audiodata;
					[m,n] = size(handles.residual.data);
					[me,ne] = size(excitation.data);
					if (me > m)
						excitation.data = excitation.data(1:m,:);
					elseif (me < m)
						while me < m
							excitation.data = [excitation.data; excitation.data];
							[me,ne] = size(excitation.data);
						end
					end
					if (ne == 2 & n == 1)
						excitation.data = to_mono(excitation.data);
					elseif (ne == 1 & n == 2)
						excitation.data = [excitation.data, ...
							excitation.data];
					end
			end
			handles.synthdata = synthesize(handles.analysisdata, excitation, ...
				windowsize, windowskip, windowshape);
			[m,n] = size(handles.audiodata.data);
			handles.synthdata.data = handles.synthdata.data(1:m,1:n);
			%if length(handles.synthdata.data) > length(handles.audiodata.data)
			%	handles.synthdata.data = ...
			%		handles.synthdata.data(1:length(handles.audiodata.data),:);
			%end
			if (max(abs(handles.synthdata.data)) > 1.0)
				handles.synthdata.data = normalize(handles.synthdata.data);
			end
			makeTimePlot(handles.synthdata, handles.timeplot2);
		end
	case 'normalize'
		if isfield(handles,'synthdata')
			handles.synthdata.data = normalize(handles.synthdata.data);
			xlim = get(handles.timeplot1,'XLim');
			makeTimePlot(handles.synthdata, handles.timeplot2);
			set(handles.timeplot2,'XLim',xlim);
		end
	case 'norm_residual'
		if isfield(handles,'residual')
			handles.residual.data = normalize(handles.residual.data);
			xlim = get(handles.resplot,'XLim');
			makeTimePlot(handles.residual, handles.resplot);
			set(handles.resplot,'XLim',xlim);
		end
	case {'playsound1','playsound2'}
		audiodata = {};
		times = get(handles.timeplot1,'XLim');
		switch action
			case 'playsound1'
				if isfield(handles,'audiodata')
					audiodata = handles.audiodata;
					playbutton = handles.play1;
				end
			otherwise
				if isfield(handles,'synthdata')
					audiodata = handles.synthdata;
					playbutton = handles.play2;
				end
		end
		if ~isempty(audiodata)
			samples = round(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));
			play_audiodata(audiodata, playbutton);
		end
	case 'play_residual'
		audiodata = {};
		times = get(handles.timeplot1,'XLim');
		if isfield(handles,'residual')
			audiodata = handles.residual;
			playbutton = handles.play_residual;
		end
		if ~isempty(audiodata)
			samples = round(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));
			play_audiodata(audiodata, playbutton);
		end
	case {'orig_fourier','orig_sonogram'}
		if isfield(handles,'audiodata'),
			audiodata = handles.audiodata;
			if size(audiodata.data,2) > 1
				audiodata.data = to_mono(audiodata.data);
			end
			switch action
				case 'orig_fourier'
					fourierexpogui(audiodata);
				case 'orig_sonogram'
					sonoexpogui(audiodata);
			end
		end
	case {'synth_fourier','synth_sonogram'}
		if isfield(handles,'synthdata'),
			audiodata = handles.synthdata;
			if size(audiodata.data,2) > 1
				audiodata.data = to_mono(audiodata.data);
			end
			switch action
				case 'synth_fourier'
					fourierexpogui(audiodata);
				case 'synth_sonogram'
					sonoexpogui(audiodata);
			end
		end
	case 'saveresynth'
		if isfield(handles, 'synthdata')
			save_audiodata(handles.synthdata);
		end
	case 'colacheck'
		cola_check(handles);
	case 'print'
		print_figure(f);
	case 'close'
		close_figure(f,figname(1:end-4));
		return;
end
set(f,'UserData',handles);


function makeTimePlot(audiodata, axesnum);
if max(max(audiodata.data)) > 1
	audiodata.data = normalize(audiodata.data);
end
axes(axesnum)
t = [0:1/audiodata.Fs:(length(audiodata.data)-1)/audiodata.Fs];

plot(t,audiodata.data)
maxtime = length(t)/audiodata.Fs;
set(axesnum,'XLim',[0 maxtime]);
set(axesnum,'YLim',[-1.0 1.0]);
grid;
xlabel('time (s)');


% --------------------------------------------------------------------
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 [analysis,residual] = analyze(audiodata, windowsize, windowskip, windowshape, order);
% Perform LPC analysis of audiodata
analysis.Fs = audiodata.Fs;
analysis.windowsize = windowsize;
%analysis.windowskip = windowskip;
windowskip = windowsize/4;
analysis.windowskip = windowskip;

analysis.windowshape = windowshape;
windowdata = get_windowdata(windowsize,windowshape);
num_channels = min(size(audiodata.data));
num_windows = floor(length(audiodata.data)/windowskip);

if (num_windows-1)*windowskip+windowsize > size(audiodata.data,1)
	rem_samples = ((num_windows-1)*windowskip+windowsize) - size(audiodata.data,1);
	audiodata.data = [audiodata.data; zeros(rem_samples,1)];
end

residual.data = zeros(size(audiodata.data,1),num_channels);
residual.Fs = audiodata.Fs;
Zf = zeros(1,order);

%bat = figure;
h = waitbar(0, 'Analyzing ...');
for i=1:num_windows, % For each window
	for k = 1:num_channels % For each channel
		start_sample = (i-1)*windowskip + 1;
		end_sample = start_sample + windowsize-1;
		sig = audiodata.data(start_sample:end_sample,k);
% 		R = xcorr(frame,'coeff');
% 		R = R(windowsize:end);
% 		R = R(1:order+1);
% 		Rt = toeplitz(R(1:end-1));
% 		alpha = inv(Rt)*R(2:end);
% 		analysis.data{i,k} = [1 -alpha'];
		a = lpc(windowdata.*sig,order);
		analysis.data{i,k} = a;
		
		% now find residual in first parts of window
		[estimate,Zf] = filter([0 -a(2:end)],1,sig,Zf);
% 		b = residual.data(start_sample+windowskip:end_sample-windowskip);
% 		c = sig(windowskip:end-windowskip-1);
% 		whos
		
% 		residual.data(start_sample+windowskip:end_sample-windowskip+1) = ...
% 			sig(windowskip:end-windowskip) - estimate(windowskip:end-windowskip);

			residual.data(start_sample+windowskip-1:end_sample-2*windowskip) = ...
				sig(windowskip:end-2*windowskip) - estimate(windowskip:end-2*windowskip);
	
% 		figure(bat);
% 		hold on;
% 		n = [windowskip:3*windowskip-1];
% 		%plot([start_sample:end_sample],sig,'k');
% 		plot([start_sample+windowskip:start_sample+3*windowskip-1], ...
% 			estimate(windowskip:end-windowskip-1),'b');
% % 		plot(n+start_sample-1,residual.data(start_sample+windowskip:end_sample-windowskip),'r');
% 		
% 		if i==2
% 			stop
% 		end
		
% 		if i==10
% 			a = lpc(frame,order)
% 			figure;
% 			freqz(1,[1 -a(2:end)])		
% 			stop
% 		end
	end
	waitbar(i/num_windows,h);
end
close(h);


% --------------------------------------------------------------------
function audiodata = synthesize(analysis, excitation, windowsize, windowskip, windowshape)
audiodata.Fs = analysis.Fs;

num_channels = min(size(excitation.data));
num_windows = length(analysis.data);
audiodata.data = zeros(size(excitation.data,1),num_channels);
windowdata = get_windowdata(windowsize,windowshape);
order = length(analysis.data{1,1})-1;
orig_length=size(excitation.data,1);
windowskip = windowsize/4;

add_samples = num_windows*windowskip + windowsize - size(excitation.data,1);
if add_samples > 0
	excitation.data = [excitation.data; zeros(add_samples,1)];
	audiodata.data = [audiodata.data; zeros(add_samples,1)];
end
Zf = zeros(1,order);

h = waitbar(0, 'Synthesizing ...');
% for i=1:num_windows
% 	for k=1:num_channels
% 		start_sample = (i-1)*windowskip + 1;
% 		end_sample = start_sample + windowsize-1;
% 		a = analysis.data{i,k};
% 
% 		[sig,Zf] = filter(1,a,excitation.data(start_sample:end_sample,k),Zf);
% 		audiodata.data(start_sample+windowskip-1:end_sample-2*windowskip) = ...
% 			sig(windowskip:end-2*windowskip);
% 		
% % 		if i==10
% % 			figure;
% % 			plot(audiodata.data);
% % 			stop
% % 		end
% 		
% % 		audiodata.data(start_sample:end_sample,k) = ...
% % 			audiodata.data(start_sample:end_sample,k) + sig(1:windowsize);
% 	end
% 	waitbar(i/num_windows,h);
% end

% Try long ways
for i=1:num_windows
	for k=1:num_channels
		start_sample = i*windowskip;
		end_sample = start_sample + windowskip-1;
		a = analysis.data{i,k};

		for j=start_sample:end_sample
			audiodata.data(j,k) = -a(2:end)*audiodata.data(j-1:-1:j-order,k) + ...
				excitation.data(j,k);
		end
		
	end
	waitbar(i/num_windows,h);
end

close(h);

function cola_check(handles)
if get(handles.analcola,'Value')
	winsize = str2double(get(handles.analwindowsize,'String'));
	% if windowshape = ...
	set(handles.analwindowskip,'String',num2str(ceil(winsize/2)));
	set(handles.analwindowskip,'Enable','inactive');
else
	set(handles.analwindowskip,'Enable','on');
end