compandexpofn.m

非常好的数字处理教程

function compandexpofn(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;
		movegui(f,'center');
		reset(handles.SQNRplot);
		reset(handles.pdfplot);
		handles.signal = getrandsig(handles);
		updatePDF(handles);
		updateSQNR(handles);
	case 'resize'
		movegui(f,'onscreen');
	case 'sqnrupdate'
		updateSQNR(handles);
		updatePDF(handles);
	case 'pdfupdate'
		handles.signal = getrandsig(handles);
		updatePDF(handles);
		updateSQNR(handles);
	case 'apply'
		handles.audiodata = load_audiodata;
		if ~isfield(handles.audiodata, 'filenamepath')
			return;
		end
		if (size(handles.audiodata.data,2) > 1)
			handles.audiodata.data = normalize(to_mono(handles.audiodata.data));
		end
		handles.qaudiodata = handles.audiodata;
		[handles.qaudiodata.data, D] = quantizedata(handles.audiodata.data, ...
			handles.quant.x, handles.quant.y);
		set(handles.play_orig,'Visible','on');
		set(handles.play_quant,'Visible','on');
	case 'play_orig'
		if isfield(handles,'audiodata')
			play_audiodata(handles.audiodata);
		end
	case 'play_quant'
		if isfield(handles,'audiodata')
			play_audiodata(handles.qaudiodata);
		end
	case 'print'
		print_figure(f);
	case 'close'
		close_figure(f,figname(1:end-4));
		return;
end
set(f,'UserData',handles);


function signal = getrandsig(handles)
	Q = 20000;

	signal.data = [];
	signal.type = '';
	
	contents = get(handles.pdfmenu,'String');
	switch lower(contents{get(handles.pdfmenu,'Value')})
		case 'uniform'
			data = 2*rand(Q,1)-1;
			signal.data = data./sqrt(var(data));
			signal.type = contents{get(handles.pdfmenu,'Value')};
		case 'gaussian'
			signal.data = randn(Q,1);
			signal.type = contents{get(handles.pdfmenu,'Value')};
		case 'laplacian'
			signal.data = exprnd(1,Q,1).*sign(randn(Q,1))./sqrt(2);
			signal.type = contents{get(handles.pdfmenu,'Value')};
		case 'audio'
			audiodata = load_audiodata;
			if ~isfield(audiodata, 'filenamepath')
				return;
			end
			if (size(audiodata.data,2) > 1)
				audiodata.data = normalize(to_mono(audiodata.data));
			end
			signal.data = audiodata.data./sqrt(var(audiodata.data));
			signal.type = audiodata.filenamepath;
		case 'audio difference'
			audiodata = load_audiodata;
			if ~isfield(audiodata, 'filenamepath')
				return;
			end
			if (size(audiodata.data,2) > 1)
				audiodata.data = normalize(to_mono(audiodata.data));
			end
			signal.data = diff(audiodata.data);
			signal.data = signal.data./sqrt(var(signal.data));
			signal.type = audiodata.filenamepath;
	end
	
	
function updatePDF(handles)

	xmax = str2num(get(handles.maxx,'String'));
	M = length(handles.signal.data)/100;
	[fx,x] = hist(handles.signal.data,M);

	normalize1 = sum(fx)*2*xmax/M;
	
	% Lowpass filter fx
	N = 3;
	b = ones(1,N)./N;
	fx = filtfilt(b,1,[ones(1,N-1).*fx(1) fx]);
	fx = fx(N:end);
	fx = [0 fx 0]./normalize1;
	x = [x(1) x x(end)];
	
	% Compress the signal and find the new pdf
	mu = str2num(get(handles.mu,'String'));
	
	lnmu = log(1+mu);
	lnAp1 = log(mu) + 1;
  VdA   = xmax/mu;
	
	contents = get(handles.compmenu,'String');
	switch lower(contents{get(handles.compmenu,'Value')})
		case 'mu-law'
			compsig = xmax*log(1+mu*abs(handles.signal.data)/xmax)/lnmu.*sign(handles.signal.data);
		case 'a-law'
			indx = find(abs(handles.signal.data) <= VdA);
			if ~isempty(indx)
					compsig(indx) = mu/lnAp1*abs(handles.signal.data(indx)).*sign(handles.signal.data(indx));
			end
			indx = find(abs(handles.signal.data) >  VdA);
			if ~isempty(indx)
					compsig(indx) = xmax/lnAp1* ...
						(1+log(abs(handles.signal.data(indx))/VdA)).*sign(handles.signal.data(indx));
			end
	end
	
	[fxc,xc] = hist(compsig,M);
	normalize2 = sum(fxc)*2*xmax/M;
	
	fxc = filtfilt(b,1,[ones(1,N-1).*fxc(1) fxc]);
	fxc = fxc(N:end);
 	fxc = [0 fxc 0]./normalize2;
	
	axes(handles.pdfplot);
	set(handles.pdfplot,'FontSize',12);
	
	plot(x,fx,'k',x,fxc,'b');
	axis([1.1*min(min(x),min(xc)) 1.1*max(max(x),max(xc)) -0.1 1.1*max(max(fx),max(fxc))]);
	xlabel('x');
	grid on;
	set(handles.pdfplot,'YTickLabel','');
	title([handles.signal.type ' PDF'],'Interpreter','none');


function updateSQNR(handles)
	Q = 1000;
	mu = str2num(get(handles.mu,'String'));
	xmax = str2num(get(handles.maxx,'String'));
	N = str2num(get(handles.levels,'String'));
	
	contents = get(handles.quantmenu,'String');
	switch lower(contents{get(handles.quantmenu,'Value')})
		case 'uniform'
			del = 2*xmax/N;
			if rem(N,2) == 0
				x = [-Inf [-N/2:N/2]*del Inf];
				y = x(2:end-2) + del/2;
				y = [y(1) y y(end)];
			else
				x = [[1:(N+1)/2]*del Inf]-del/2;
				x = [-fliplr(x) x];
				y = x(2:end-2) + del/2;
				y = [y(1) y y(end)];
			end
		case 'non-uniform'
			contents = get(handles.pdfmenu,'String');
			switch lower(contents{get(handles.pdfmenu,'Value')})
				case 'gaussian'
					[x,y,D,H] = qcgauss(N,1,10e-7,'optimal');
				case 'laplacian'
					[x,y,D,H] = qclaplace(N,1,10e-7,'optimal');
				otherwise
					del = 2*xmax/N;
					if rem(N,2) == 0
						x = [-Inf [-N/2:N/2]*del Inf];
						y = x(2:end-2) + del/2;
						y = [y(1) y y(end)];
					else
						x = [[1:(N+1)/2]*del Inf]-del/2;
						x = [-fliplr(x) x];
						y = x(2:end-2) + del/2;
						y = [y(1) y y(end)];
					end
			end
	end

	
	vari = 10.^([-60:1:20]/10); 
	csqnr = [];
	usqnr = [];
	randsig = handles.signal.data(1:Q,1);
	
	contents = get(handles.pdfmenu,'String');
	switch lower(contents{get(handles.pdfmenu,'Value')})
		case 'uniform'
			tit = 'Signal-to-Quantization Noise for Uniform PDF';
		case 'gaussian'
			tit = 'Signal-to-Quantization Noise Gaussian PDF';
		case 'laplacian'
			tit = 'Signal-to-Quantization Noise Laplacian PDF';
		otherwise
			tit = handles.signal.type;
	end
	
	lnmu = log(1+mu);
	contents = get(handles.compmenu,'String');
	lnAp1 = log(mu) + 1;
  VdA   = xmax/mu;
		
	for j=1:length(vari)
		signal = sqrt(vari(j))*randsig;
		switch lower(contents{get(handles.compmenu,'Value')})
			case 'mu-law'
				compsig = xmax*log(1+mu*abs(signal)/xmax)/lnmu.*sign(signal);
			case 'a-law'
				indx = find(abs(signal) <= VdA);
				if ~isempty(indx)
            compsig(indx) = mu/lnAp1*abs(signal(indx)).*sign(signal(indx));
        end
        indx = find(abs(signal) >  VdA);
        if ~isempty(indx)
            compsig(indx) = xmax/lnAp1* ...
							(1+log(abs(signal(indx))/VdA)).*sign(signal(indx));
        end
		end
		[qcsignal,D] = quantizedata(compsig,x,y);
		[qusignal,D] = quantizedata(signal,x,y);

		switch lower(contents{get(handles.compmenu,'Value')})
			case 'mu-law'
				creproduced = xmax/mu* ...
					(exp(abs(qcsignal)*lnmu/xmax)-1).*sign(qcsignal);
			case 'a-law'
				VdlnAp1 = xmax/lnAp1;
        indx = find(abs(qcsignal) <= VdlnAp1);
        if ~isempty(indx)
            creproduced(indx,1) = lnAp1/mu* ...
							abs(qcsignal(indx)).*sign(qcsignal(indx));
        end
        indx = find(abs(qcsignal) >  VdlnAp1);
        if ~isempty(indx)
            creproduced(indx,1) = VdA* ...
							exp(abs(qcsignal(indx))/VdlnAp1 - 1).*sign(qcsignal(indx));
				end
		end
		
		qcerror = creproduced - signal;
		querror = qusignal - signal;
		csqnr=[csqnr 10*log10(vari(j)/mean(qcerror.^2))]; %sample mean
		usqnr=[usqnr 10*log10(vari(j)/mean(querror.^2))]; %sample mean
	end
	
	axes(handles.SQNRplot);
	set(handles.SQNRplot,'FontSize',12);
	if get(handles.holdon,'Value')
		hold on;
	else
		hold off;
	end
	a = plot(10*log10(vari),csqnr,'b');
	hold on;
	b = plot(10*log10(vari),usqnr,'k');
	hold off;
	axis([-60 15 min(-30,1.1*max(min(csqnr,usqnr))) max(45,1.1*max(max(csqnr,usqnr)))]);
	xlabel('Signal Power (dB)');
	ylabel('SQNR (dB)');
	grid on;
	legend([b,a],'No Companding','Companding',4);
	
	title(tit,'Interpreter','none');
	
	
	
function [qdata, D] = quantizedata(data, x, y)

	qdata = zeros(length(data),1);
	N = length(y);	
	for i = 1:N
		ind = find(data >= x(i) & data < x(i+1));
		qdata(ind,1) = y(i);
	end

	D = 0;