log[1??y[n]?Xmaxx[n]Xmax?sign(x[n])
log(1??)]给出的mulaw1()函数是最大数值假设为Xmax=1的信号向量上实现该μ律压扩器。
（1）为实现上图的系统，必须为逆μ律压扩器编写一个M文件。这个M文件应有下面的调用顺序和参数
Function x=mulawinv(y,mu)
（2）实现一个6bit的μ律量化器，即压缩的样本用6个bit表示，考察相应的量化误差，并用wavwrite()函数写到WAV文件，倾听量化后的声音。
>> y=mulaw1(x,255);
>> yh=fxquant(y,6,'round','sat'); 4、信噪比
比较量化器的一个方便方法是计算信号功率与量化噪声功率的比值。对于matlab的实验，SNR的一个方便定义是
SNR?10log(L?1n?0)
?[n]?x[n])2?(xn?0?(x[n])L?12 （1）使用提供的snr()函数各个bit下均匀量化的SNR，计算出的结果是否
与预计的数量不同？
>> [snratio E]=snr(x1,x);
（2）均匀量化和μ律量化的比较
给出的qplot()函数对一个固定量化器的均匀量化与μ律量化进行比较，而该固定量化器以减小的振幅（以因子2减小）作为输入。
绘出的图应表明，μ律量化器在约64 :1的输入幅度范围上保持一个常数信噪比。均匀量化器需要多少个比特数以在相同的区间上至少维持与6bit的μ律量化器相同的信噪比？
>> qplot(x,6,255,7);
附录： （1）
function [P, c] = welch(x, N, M, wintype, nplot) %WELCH Power Spectrum Estimation by Welch's method %-----
% Usage: [P,c] = welch(x,N,M,wintype,nplot) %
% P : power spectrum by Welch's method
% c : correlation function = inverse of Welch power spectrum % x : input signal vector % N : FFT length % M : window length
% wintype = 'r' for rectangular
% otherwise, Hamming window is used
% nplot = omitting this argument suppresses plotting
%--------------------------------------------------------------- % copyright 1994, by C.S. Burrus, J.H. McClellan, A.V. Oppenheim, % T.W. Parks, R.W. Schafer, & H.W. Schussler. For use with the book % \% (Prentice-Hall, 1994).
%---------------------------------------------------------------
x = x(:);
if(wintype=='r') w = ones(M,1); else
w = hamming(M); end
K = fix((length(x)-M/2)/(M/2)); KLU = K*sum(w.^2); P = zeros(N,1); for i=1:K
ibeg = (M/2)*(i-1) + 1;
X = fft(x(ibeg:ibeg+M-1).*w, N); P = P + abs(X).^2; end
P = P/KLU; c = ifft(P);
c = real(c(1:M/2)); P = P(1:N/2+1); if( nargin==5 )
plot((0:N/2)/(N/2),P)
axis([0 1 0 max(P)]) %<-- version 4.0
title(['Welch Power Spectrum: Window length=',... num2str(M), ' Number sections=', num2str(K)]) xlabel('FREQUENCY (omega/pi)') pause %---
plot((0:N/2)/(N/2),10*log10(P))
title(['Welch Log Power Spectrum: Window length=',... num2str(M), ' Number sections=', num2str(K)]) xlabel('FREQUENCY (omega/pi)') pause %---
stem(0:M/2-1, c) %<-- version 4.0
title(['Welch Correlation Function: Window length=',... num2str(M), ' Number sections=', num2str(K)]) xlabel('LAG INDEX') end （2）
function X = fxquant( s, bit, rmode, lmode ) %FXQUANT simulated fixed-point arithmetic %-------
% Usage: X = fxquant( S, BIT, RMODE, LMODE ) %
% returns the input signal S reduced to a word-length
% of BIT bits and limited to the range [-1,1). The type of % word-length reduction and limitation may be chosen with % RMODE: 'round' rounding to nearest level % 'trunc' 2's complement truncation % 'magn' magnitude truncation % LMODE: 'sat' saturation limiter
% 'overfl' 2's complement overflow % 'triang' triangle limiter % 'none' no limiter
%--------------------------------------------------------------- % copyright 1994, by C.S. Burrus, J.H. McClellan, A.V. Oppenheim, % T.W. Parks, R.W. Schafer, & H.W. Schussler. For use with the book % \% (Prentice-Hall, 1994).
%---------------------------------------------------------------
if nargin ~= 4;
error('usage: fxquant( S, BIT, RMODE, LMODE ).'); end;
if bit <= 0 | abs(rem(bit,1)) > eps;
error('wordlength must be positive integer.'); end;
Plus1 = 2^(bit-1);
X = s * Plus1;
if strcmp(rmode, 'round'); X = round(X); elseif strcmp(rmode, 'trunc'); X = floor(X); elseif strcmp(rmode, 'ceil'); X = ceil(X); elseif strcmp(rmode, 'magn'); X = fix(X);
else error('unknown wordlength reduction spec.'); end;
if strcmp(lmode, 'sat'); X = min(Plus1 - 1,X); X = max(-Plus1,X);
elseif strcmp(lmode, 'overfl');
X = X + Plus1 * ( 1 - 2*floor((min(min(X),0))/2/Plus1) ); X = rem(X,2*Plus1) - Plus1; elseif strcmp(lmode, 'triang');
X = X + Plus1 * ( 1 - 2*floor((min(min(X),0))/2/Plus1) ); X = rem(X,4*Plus1) - Plus1; f = find(X > Plus1); X(f) = 2*Plus1 - X(f); f = find(X == Plus1); X(f) = X(f) - 1;
elseif strcmp(lmode, 'none'); % limiter switched off else error('unknown limiter spec.'); end;
X = X / Plus1; （3）
function y = mulaw1(x, mu)
%MULAW mu-law compression for signals with %----- assumed maximum value of 32767 %
% Usage: y = mulaw(x, mu); %
% x : input signal, column vector with max value 32767 % mu : compression parameter (mu=255 used for telephony) %
% see also MULAWINV
%--------------------------------------------------------------- % copyright 1994, by C.S. Burrus, J.H. McClellan, A.V. Oppenheim, % T.W. Parks, R.W. Schafer, & H.W. Schussler. For use with the book % \% (Prentice-Hall, 1994).
%---------------------------------------------------------------
Xmax = 1;%wang:32767;
y = (Xmax/log(1+mu))*log(1+mu*abs(x)/Xmax).*sign(x); （4）
function [s_n_r,e] = snr(xh, x) e=xh-x;
sumX=sum(x.*x);
sumE=sum(e.*e);
s_n_r=10*log10(sumX/sumE); （5）
function qplot(s, nbits, mu, ncases)
%QPLOT for plotting dependence of signal-to-noise ratio %----- on decreasing signal level %
% Usage: qplot(s, nbits, mu, ncases) %
% s : input test signal
% nbits : number of bits in quantizer % mu : mu-law compression parameter % ncases : number of cases to plot %
% NOTE: assumes ROUNDING for quantizer
% and requires user-written MULINV and SNR %
% see also MULAW and FXQUANT
%--------------------------------------------------------------- % copyright 1994, by C.S. Burrus, J.H. McClellan, A.V. Oppenheim, % T.W. Parks, R.W. Schafer, & H.W. Schussler. For use with the book % \% (Prentice-Hall, 1994).
%---------------------------------------------------------------
P = zeros(ncases,2); x = s;
for i=1:ncases
sh = fxquant(x, nbits, 'round', 'sat'); P(i) = snr(sh,x); y = mulaw1(x, mu);
yh = fxquant(y, nbits, 'round', 'sat'); xh =mulawinv(yh, mu); Q(i) = snr(xh,x); x = x/2; end
plot(P); hold on; plot(Q,'r'); hold off;
title(['SNR for ', num2str(nbits), '-bit Uniform and ',... num2str(mu), '-Law Quantizers'])
xlabel('power of 2 divisor'); ylabel('SNR in dB')