www.gusucode.com > signal 工具箱matlab源码程序 > signal/+fdesign/@abstracttype/supernoiseshape.m
function nsres = supernoiseshape(this,b,linearphase,wl,cb,f,a,args) %SUPERNOISESHAPE <short description> % Copyright 2008-2011 The MathWorks, Inc. % Process the arguments to set up the respective results structs, get % the number of iterations to loop, and find which metric to use. nsres = processargs(b,linearphase,cb,f,a,args); % Quantize coeffs to WL bits warnstate = warning('off','fixed:fi:overflow'); bfi = fi(b,1,wl); fl = bfi.FractionLength; bq = double(bfi); bp = bq; k = nsres.trials; aux = fi(zeros(k,1),1,wl,fl); % Initialize the temporary results vector to zeros and the initial % pareto frontier matrix with the quantized filter results. bands = noiseshapeparetobands(this); nbands = length(bands); m = zeros(k,nbands+1); p = [zeros(1,nbands) length(b)]; for i=1:nbands, % Lowpass, highpass uses {'cb','nb'} while halfband uses {'cb'} for % example p(i) = measurefreqz(bq,nsres,bands{i}); end % This runs the noise shaping algorithm using a noise shaping % filter of length 3, then of length 4. These values have % experimentally shown to be the best lengths, but future research % could be done to see if higher values give better results. refcoeffs = nsres.refcoeffs; lsbq = double(lsb(bfi)); nb = length(b); n = nsres.n; w = nsres.w; metric = nsres.metric; pnorm = nsres.pthnorm; for nnsf = 3:4 % The first thing we do is compute a noise shaping filter of % the specified length and critical band. This is generally % just a linear phase remez filter, with the pass band set to % the critical band. nsf = getnoiseshapefilter(this,nnsf,cb); % This noise shaping algorithm requires that the noise shaping % filter be scaled to have the first coefficient be 1, and then % have the first coefficient removed from the filter structure. nsf = -1*nsf(2:end)'./nsf(1); lnsf = length(nsf); % Initialize the filter state to be random uniform noise % within plus or minus 1 LSB. r = (rand(k,lnsf)-0.5)*lsbq; % Initialize the new filter coefficients to zero for now. z = zeros(k,n); % This is the main noise shaper here; essentially, it is % just a simple FIR digital filter with a quantization % after each time step. It takes in the original filter % coefficient as the input and gives noise shaped filter % coefficients as the output. It uses the 'nsf' (noise % shaping filter) to move the quantization noise out of the % critical band. for j = 1:1:n x1 = refcoeffs(j)+r*nsf; aux(:) = x1; % This is the fastest way to quantize z(:,j) = double(aux); r(:,2:end) = r(:,1:(end-1)); r(:,1) = x1-z(:,j); end % Ok, now that the noise shaping is done we need to compute % the metric in both the critical band ('m(:,1)') and the % non-critical band(s) ('m(:,2)'). % % These MxN matrices simple compute an M point FFT on an N point % input. For example, w1*z computes the M point FFT on N point % vector z in the critical band (w2 is for the non-critical band). for i = 1:nbands, res = nsres.(['res_' bands{i}]); m(:,i) = metric(w{i}*z',res.des,pnorm)'; end m(:,end) = length(b); % If the filter was linear phase, mirror the coefficients. if linearphase z = [fliplr(z) z(:,(1+mod(nb,2)):end)]; end % Only keep the results that lie on the pareto frontier to % save memory and time. [p,ix] = paretofrontier([p;m]); % Keep the filter coefficients for points on the pareto % front. z2 = [bp;z]; bp = z2(ix,:); end nsres = nsresaddresults(nsres,bp,bands); % Restore warning warning(warnstate); end %-------------------------------------------------------------------------- function nsres = processargs(b,linearphase,cb,f,a,args) % Define frequency mask [res1,res2] = definefreqmask(cb,f,a); grid = [res1.fgrid; res2.fgrid]; des = [res1.des; res2.des]; [grid,ix] = sort(grid,'ascend'); nsres.res.fgrid = grid; nsres.res.des = des(ix); nsres.res_cb = res1; nsres.res_nb = res2; % Set the number of trials per batch. nsres.trials = 1000; % Define metric switch args.noiseShapeNorm, case inf nsres.metric = @maximumripple; case 2 nsres.metric = @leastsquares; otherwise nsres.metric = @leastpthnorm; end nsres.pthnorm = args.noiseShapeNorm; % Compute DFT constants j = sqrt(-1); n = size(b,2); w1 = sum(exp(-1*j*pi*bsxfun(@times,nsres.res_cb.fgrid,0:(n-1))),3); w2 = sum(exp(-1*j*pi*bsxfun(@times,nsres.res_nb.fgrid,0:(n-1))),3); if linearphase % If the input filter has a linear phase numerator, then the % algorithm will preserve linear phase. The advantage here is % that optimizations can be used to make noise shaping run % faster. odd = mod(n,2); n = ceil(n/2); refcoeffs = b((n+(~odd)):end); if odd w1 = w1(:,n:end)+[zeros(size(w1,1),1) fliplr(w1(:,1:(n-1)))]; w2 = w2(:,n:end)+[zeros(size(w2,1),1) fliplr(w2(:,1:(n-1)))]; else w1 = w1(:,(n+1):end)+fliplr(w1(:,1:n)); w2 = w2(:,(n+1):end)+fliplr(w2(:,1:n)); end else % Non-linear phase, so no optimizations. refcoeffs = b; end nsres.w = {w1;w2}; nsres.refcoeffs = refcoeffs; nsres.n = n; end function [res1,res2] = definefreqmask(cb,f,a) r = 250; res.fgrid = zeros(length(f)/2*r,1); res.des = res.fgrid; % Mask (assumes piecewise constant amplitude vector a) for i = 1:length(f)/2 res.fgrid((((i-1)*r)+1):(i*r)) = linspace(f(2*i-1),f(2*i),r); res.des ((((i-1)*r)+1):(i*r)) = linspace(a(2*i-1),a(2*i),r); end % We don't want to measure the frequency response at pi idx = res.fgrid~=1; res.fgrid = res.fgrid(idx); res.des = res.des(idx); cfgridnum = sum((res.fgrid >= cb(1)) & (res.fgrid <= cb(2))); res1.fgrid = zeros(cfgridnum,1); res2.fgrid = zeros(length(res.fgrid)-cfgridnum,1); fnames = {'fgrid','des'}; for f = 1:1:length(fnames) if length(res.(fnames{f})) ~= length(res.fgrid) continue; end res1.(fnames{f}) = res1.fgrid; res2.(fnames{f}) = res2.fgrid; j = 1; k = 1; for i = 1:1:length(res.fgrid) if (res.fgrid(i) >= cb(1)) && (res.fgrid(i) <= cb(2)) res1.(fnames{f})(j) = res.(fnames{f})(i); j = j+1; else res2.(fnames{f})(k) = res.(fnames{f})(i); k = k+1; end end end end function m = measurefreqz(b,nsres,band) res = nsres.(['res_' band]); m = nsres.metric(freqz(b,1,pi*res.fgrid),res.des,nsres.pthnorm); end function m = maximumripple(x,des,pnorm) %#ok<INUSD> % This function determines that maximum difference between array 'x' % and array 'des'. It automatically expands singleton dimensions and % returns an array of values. m = 20*log10(max(abs(bsxfun(@minus,abs(x),des))))'; end function m = leastsquares(x,des,pnorm) %#ok<INUSD> % This function determines the mean square error between array 'x' and % array 'des'. It automatically exands singleton dimensions and % returns an array of values. m = sqrt(sum(bsxfun(@minus,abs(x),des).^2)'/length(des)); end function m = leastpthnorm(x,des,pnorm) % This function determines the P-th norm between array 'x' and array % 'des'. L = size(x,2); m = zeros(L,1); y = bsxfun(@minus,abs(x),des); for i = 1:1:L m(i) = norm(y(:,i),pnorm); end end function [p,ix] = paretofrontier(d) % This function computes the pareto frontier of an MxN dataset 'd'. % It returns the sorted points on the frontier as 'p' and the indexes % of those points in 'd' as 'ix'. A point is considered to be on the % pareto frontier if no other point exists with all N dimensions lower % than it. % Compute the pareto frontier efficiently [d1,m] = unique(d,'rows'); p = true(size(d1,1),1); for i=1:size(d1,1) if p(i) p = p & ~(all(bsxfun(@ge,d1,d1(i,:)),2) & ((1:size(d1,1))'~=i)); end end ix = m(p); p = d1(p,:); [p1,pix] = sort(p(:,1),'ascend'); p = p(pix,:); ix = ix(pix); end function nsres = nsresaddresults(nsres,bp,bands) % bp contains filters on the pareto front % Keep filter with the best attenuation in the critical band (i.e. % stopband) nsres.filters.bns = bp(1,:); nbands = length(bands); for i = 1:nbands, nsres.filters.metrics.(['bns_' bands{i}]) = measurefreqz(bp(1,:),nsres,bands{i}); end nsres.filters = orderfields(nsres.filters); nsres = orderfields(nsres); end % [EOF]