% MULTI-STAGE DIFFERENTIAL LPCM CODING (order 1 to k) (lpcmdif.m) % 3-stage transparent differential coding for lossless compression % 2 clipping options including error correction to the 3rd integral (5 pulse substitution) % set all parameters in program header prior to execution % LPCM reference spectrum quantized to qin bit with an assumed sampling rate fs Hz % plots are saved in Temp directory on C: drive if sav = 1 in program header % April 2001 % MATLAB 5.2 compliantclose; clear; home; clc; colordef white; fprintf('MULTI-STAGE DIFFERENTIAL LPCM CODING (order 1 to k)\n') fprintf('3-stage transparent differential coding for lossless compression\n') fprintf('Error dispersion/correction operates under overload\n') fprintf('*******************************************************************\n')
% control and input parameters order=6; %noise shaper order range 1 to k (6) sec=3; %differential compressor order 1, 2 or 3 stages (2) adjc=0; %adjacent sample smoothing after compression (1 on) clip=2; %clipping: clip = 0 off, clip = 1 three-pulse substitution, clip = 2 five-pulse substitution pass=5; %maximum number of passes in error correction loop (5) A1=10000; %input amplitude A1 A2=10000; %input amplitude A2 fin1=1000; %input frequency of A1, Hz fin2=20000; %input frequency of A2, Hz fs=96000; %digital audio sampling frequency (96 kHz) m=4; %oversampling ratio (relative to fs) (4) v=16; %vector length, bits (16) qout=11; %output word length of final differentiator (11); determines clipping level if clip > 0 qin=24; %reference signal resolution, input range +1 to -1, bits (24), sampling rate set at fs Hz flo=20000; %flo pre-emphasis lower corner frequency fhi=8*fs; %fhi pre-emphasis upper corner frequency power=5; %order of equalizer {power < 1 de-activates pre-emhasis/de-emphasis} (5) rng=500; %time-domain plot range (500) zero=0; %plot zero-order spectrum (zero=1 to plot) win=1; %window output (win = 1 raised cosine, win = 0 rectangular) sav=0; %select figure save option: sav = 1 Qin=0; %quantize input encoder input to qin bits if Qin = 1, else 0 tit=0; %put tit = 1 to plot captions on Figures, else 0
% vector length n n=2^(v+1); nm=1:n; n2=2:n/2;
% fundamental frequency in spectrum, f0 f0=m*fs/(2^v);
% triangular PDF random noise vector for dither rnd=rand(1,n)+rand(1,n)-1;
% pre-emphasised input signal to coder q1=(m^.5)*2^(qin-1); har1=ceil(fin1/f0); har2=ceil(fin2/f0); wt0=2*pi/2^v; wt1=har1*wt0; wt2=har2*wt0; if power<1 A1p=A1; A2p=A2; A=ones(size(1:2^v)); else % 3-dB break frequencies for pre-emphasis/de-emphasis fl=flo*(2^(1/power)-1)^-.5; fh=fhi*(2^(1/power)-1)^+.5; A1p=A1*((1+(f0*har1/fl)^2)/(1+(f0*har1/fh)^2))^(power*.5); A2p=A2*((1+(f0*har2/fl)^2)/(1+(f0*har2/fh)^2))^(power*.5); A=((1+(f0*(1:2^v)/fl).^2)./(1+(f0*(1:2^v)/fh).^2)).^(power*.5); end in(nm)=A1p*sin(wt1*nm-1)+A2p*sin(wt2*nm-1); if Qin==1 mx=max(abs(in)); in=mx*round(q1*in/mx+rnd(1:n))/q1; end
% LPCM reference input (amplitude range +1 to -1) ina quantized to qin bit re-sampling rate fs ina=round(q1*sin(wt0*nm-1)+rnd(1:n))/q1;
% display system parameters clipl=2^(qout-1); fprintf('\nPCM channel sampling rate, fs (MHz)') LPCM_sample_rate=fs/1e1^6 fprintf('\nPCM channel bit rate relative to m*fs (Mbit/s)') bit_rate=qin*fs*m/1e1^6 fprintf('\nNoise shaper order') order fprintf('\nOversampling ratio') m fprintf('\nChannel sampling rate (Mbit/s)') sampling_rate=fs*m/1e1^6 fprintf('\nChannel word length (bit)') word_length=qout fprintf('\nInput reference word length (bit)') word_length=qin fprintf('\nChannel bit rate (Mbit/s)') bit_rate=qout*fs*m/1e1^6
% noise-shaping routine (zero order) out0=round(in(nm)+rnd(nm));
% noise-shaping routine (selectable, maximum of kth order) outk=zeros(size(nm)); I(1:(order+1),nm)=0; for z=n/2-20:n I(1,z)=I(1,z-1)+in(z-1)-outk(z-1); for x=2:order I(x,z)=I(x-1,z)+I(x,z-1); end outk(z)=round((sum(I(1:order,z)')+in(z))+rnd(z)); end outk(n)=outk(n)-sum(outk(1+n/2:n)); outkx=outk; clear I outks(1:n)=outkx(1:n); outks(1)=0;
% calculate compressed residual using difference coding: stage 1 f1=zeros(size(nm)); home; clc fprintf('COMPRESSING NOISE SHAPER OUTPUT SIGNAL: STAGE 1\n') fout1=zeros(size(1:n)); fout1(2:n)=outks(2:n)-outks(1:n-1); foutk=fout1; f1=fout1;
% calculate compressed residual using difference coding: stage 2 f2=zeros(size(nm)); if sec>1 fprintf('COMPRESSING NOISE SHAPER OUTPUT SIGNAL: STAGE 2\n') fout2=zeros(size(1:n)); fout2(2:n)=fout1(2:n)-fout1(1:n-1); foutk=fout2; f2=fout2; end
% calculate compressed residual using difference coding: stage 3 if sec>2 fprintf('COMPRESSING NOISE SHAPER OUTPUT SIGNAL: STAGE 3\n') fout3=zeros(size(1:n)); fout3(2:n)=fout2(2:n)-fout2(1:n-1); foutk=fout3; end
% sum adjacent samples of final differentiator stage output if adjc==1 foutk(1)=0; foutk(2:n)=foutk(1:n-1)+foutk(2:n); outks(1)=0; outks(2:n)=outks(1:n-1)+outks(2:n); end
% Channel clipping: no clipping if clip==0 fprintf('\nNo clipping consequently no error correction\n\n') ERR=zeros(size(1:n)); ERRCK=ERR; ERR0=ERR; DET(2^v:n)=(sign(abs(foutk(2^v:n))-clipl)+1)/2; ERR0(10:n-2)=(foutk(10:n-2)-clipl*sign(foutk(10:n-2))).*DET(10:n-2); fk=foutk; end
% clipping with error correction: TYPE 1 [.5*e clipl-e .5*e] pulse substitution if clip==1 err=1; sx=1; ERR0=zeros(size(1:n)); while err==1 ERR=zeros(size(1:n)); ERRCK=ERR; DET=ERR; DET(n/2+1:n-1)=(sign(abs(foutk(n/2+1:n-1))-clipl)+1)/2; ERR(10:n-2)=(foutk(10:n-2)-clipl*sign(foutk(10:n-2))).*DET(10:n-2); if sx==1 ERR0(10:n-2)=(foutk(10:n-2)-clipl*sign(foutk(10:n-2))).*DET(10:n-2); end ERRO=ERR-round(ERR/2)*2; foutk=foutk+ERRO; ERR=ERR-ERRO; foutk=(1-DET).*foutk+DET.*(ERRO+sign(foutk)*clipl); foutk(2:n-1)=foutk(2:n-1)+ERR(3:n)/2+ERR(1:n-2)/2; if max(abs(DET))>0 fprintf('\nActivating clipping error correction loop sx: \n\n') sx end fk=foutk;
% clipping check sx=sx+1; err=0; for p=2^v:n if abs(foutk(p))>clipl ERRCK(p)=-clipl*sign(foutk(p))+foutk(p); err=1; end; end; if sx>abs(pass) err=0; end; end; end
% clipping with error correction: TYPE 2 [-(1/6)*e (2/3)*e clipl-e (2/3)*e -(1/6)*e] pulse substitution % note e=6 always if clip==2 err=1; sx=1; ERR0=zeros(size(1:n)); ERRQ=ones(size(1:n)); while err==1 DET=zeros(size(1:n)); ERRCK=zeros(size(1:n)); signf=zeros(size(1:n)); DET(n/2+1:n-1)=(sign(abs(foutk(n/2+1:n-1))-clipl)+1)/2; signf=sign(foutk).*DET; if sx==1 ERR0(10:n-2)=(foutk(10:n-2)-clipl*sign(foutk(10:n-2))).*DET(10:n-2); end ERRQ(10:n-2)=1+abs(round(((foutk(10:n-2)-clipl*sign(foutk(10:n-2))).*DET(10:n-2))/6)); foutk(10:n-2)=foutk(10:n-2)-6*ERRQ(10:n-2).*signf(10:n-2); foutk(9:n-3)=foutk(9:n-3)+4*ERRQ(10:n-2).*signf(10:n-2); foutk(11:n-1)=foutk(11:n-1)+4*ERRQ(10:n-2).*signf(10:n-2); foutk(8:n-4)=foutk(8:n-4)-ERRQ(10:n-2).*signf(10:n-2); foutk(12:n)=foutk(12:n)-ERRQ(10:n-2).*signf(10:n-2); if max(abs(DET))>0 fprintf('\nActivating clipping error correction loop number sx: \n\n') sx end fk=foutk; % clipping check sx=sx+1; err=0; for p=n/2+1:n-3 if abs(foutk(p))>clipl ERRCK(p)=-clipl*sign(foutk(p))+foutk(p); err=1; end; end; if sx>abs(pass) err=0; end; end; end; clear ERRQ
% reconstruct compressed signal fout1 using digital integrator: stage 1 fprintf('DE-COMPRESSING SIGNAL: STAGE 1\n') fout1=zeros(size(1:n)); for p=2:n fout1(p)=foutk(p)+fout1(p-1); end; outkr=fout1;
% reconstruct compressed signal outkr using digital integrator: stage 2 if sec>1 fprintf('DE-COMPRESSING SIGNAL: STAGE 2\n') fout2=zeros(size(1:n)); for p=2:n fout2(p)=fout1(p)+fout2(p-1); end; outkr=fout2; end;
% reconstruct compressed signal outkr using digital integrator: stage 3 if sec>2 fprintf('DE-COMPRESSING SIGNAL: STAGE 3\n') fout3=zeros(size(1:n)); for p=2:n fout3(p)=fout2(p)+fout3(p-1); end; outkr=fout3; end;
% re-quantization test to check correctness of reconstructed data outkrt=round(outkr); if abs(sum(outkrt(n/2+1:n-3)-outkr(n/2+1:n-3)))>0 fprintf('\nFailed channel re-quantization test\n') else fprintf('\nPassed channel re-quantization test\n') end outkr=outkrt; clear outkrt
% de-emphasis performed only in frequency domain, time domain left unfiltered out=outkr;
% compensate input for summing adjacent samples (multiply ina by 2) if adjc==1 ina=ina*2; end
% delete preamble sequence and set average value of sequences to zero off=10^(-300/20); nk=n/2; p(1:2^v)=out0(nk+1:n); p=p-mean(p)+off; out0=p; p(1:2^v)=outkx(nk+1:n); p=p-mean(p)+off; outkx=p; p(1:2^v)=out(nk+1:n); p=p-mean(p)+off; out=p; p(1:2^v)=outks(nk+1:n); p=p-mean(p)+off; outks=p; p(1:2^v)=outkr(nk+1:n); p=p-mean(p)+off; outkr=p; p(1:2^v)=ina(nk+1:n); p=p-mean(p)+off; ina=p; p(1:2^v)=f1(nk+1:n); p=p-mean(p)+off; f1=p; p(1:2^v)=f2(nk+1:n); p=p-mean(p)+off; f2=p; p(1:2^v)=fk(nk+1:n); p=p-mean(p)+off; fk=p; p(1:2^v)=ERR0(nk+1:n); ERR0=p; p(1:2^v)=ERRCK(nk+1:n); ERRCK=p; clear p off
% redefine vector lengths to eliminate preamble and set spectral plot span n=2^v; nm=1:n; n2=2:n/2; fr=f0*n2; nb=2:round(fs/f0); if m==1 nb=n2; end fb=nb*f0;
% window time domain waveforms wind=ones(size(nm)); wx=-8+.5*n; if win==1 wind=[.5*(1-cos(pi*(0:wx-1)/(wx+1))) wind(wx:n-wx-1) .5*(1-cos(pi*(wx-1:-1:0)/(wx+1)))]; end
% fourier transforms and de-emphasis zinfs=abs(fft(ina.*wind))/(n/2); zoutf0=abs(fft(out0.*wind))/(n/2); zoutf=abs(fft(out.*wind))/(n/2); infs(n2)=20*log10(zinfs(n2)); outf0(n2)=20*log10(zoutf0(n2)./A(n2)); outf(n2)=20*log10(zoutf(n2)./A(n2)); clear zinfs zoutf0 zoutf
% normalise peak infs to 0 dB infs(n2)=infs(n2)-max(infs(n2));
% normalise peak outf to 0 dB outf(n2)=outf(n2)-max(outf(n2));
% full spectral plots 0 to mfs/2 maxin=round(20*max(infs(n2)))/20; minin=round(20*min(infs(n2)))/20; if maxin<0 maxin=0; end axis([fr(2) fr(n/2-1) minin maxin]) hold on; grid if zero==1 plot(fr,outf0(n2),'g') end plot(fr,infs(n2),'cyan') plot(fr,outf(n2),'k') plot(fr,-20*log10(A(n2)/A(1)),'b') if tit==1 if zero==1 title('Spectrum: cyan input, green zero-order, black nth-order, blue de-emphasis') else title('Spectrum: cyan input, black nth-order, blue de-emphasis') end; end; pause if sav==1 print -dbitmap c:\Temp\'Figure 1' end close; home; clc
% expanded spectral plots 0 to fs if m>1 axis([fb(2) max(fb) minin maxin]) hold on; grid if zero==1 plot(fb,outf0(nb),'g') end plot(fb,infs(nb),'cyan') plot(fb,outf(nb),'k') plot(fb,-20*log10(A(nb)/A(1)),'b') if tit==1 if zero==1 title('EXPANDED: Spectrum: cyan input, green zero-order, black nth-order, blue de-emphasis') else title('EXPANDED: Spectrum: cyan input, black nth-order, blue de-emphasis') end; end; pause if sav==1 print -dbitmap c:\Temp\'Figure 2' end close; home; clc end
% time domain plot of noise shaper plot(outkx(1:rng),'k') if tit==1 title('Time-domain output of noise-shaped quantizer') end; pause if sav==1 print -dbitmap c:\Temp\'Figure 3' end close; home; clc
% time domain plot of first differentiator output if sec==1 plot(fk(1:rng),'k') hold plot(sign(ERRCK(1:rng)),'r') plot(clipl*ones(size(1:rng)),'k') plot(-clipl*ones(size(1:rng)),'k') plot(zeros(size(1:rng)),'k') if tit==1 if adjc==1 title('Compressed time-domain output (2 sample average): stage 1') else title('Compressed time-domain output (no average): stage 1') end; end; pause if sav==1 print -dbitmap c:\Temp\'Figure 4(sec=1)' end close; home; clc end
% time domain plot of second differentiator output if sec==2 plot(f1(1:rng),'k') hold plot(sign(ERRCK(1:rng)),'r') plot(fk(1:rng),'b') plot(clipl*ones(size(1:rng)),'k') plot(-clipl*ones(size(1:rng)),'k') plot(zeros(size(1:rng)),'k') if tit==1 if adjc==1 title('Compressed time-domain output (2 sample average): black stage 1, blue stage 2, red overload') else title('Compressed time-domain output (no average): black stage 1, blue stage 2, red overload') end; end; pause if sav==1 print -dbitmap c:\Temp\'Figure 4(sec=2)' end close; home; clc end
% time domain plot of third differentiator output if sec==3 plot(f1(1:rng),'k') hold plot(f2(1:rng),'b') plot(sign(ERRCK(1:rng)),'r') plot(fk(1:rng),'cyan') plot(clipl*ones(size(1:rng)),'k') plot(-clipl*ones(size(1:rng)),'k') plot(zeros(size(1:rng)),'k') if tit==1 if adjc==1 title('Compressed time-domain output (2 sample average): black stage 1, blue stage 2, cyan stage 3, red overload') else title('Compressed time-domain output (no average): black stage 1, blue stage 2, cyan stage 3, red overload') end; end; pause if sav==1 print -dbitmap c:\Temp\'Figure 4(sec=3)' end close; home; clc end
% time domain output of reconstructed signals plot(outkr(1:rng),'k') hold on plot(outks(1:rng)-outkr(1:rng),'r') if tit==1 title('De-compressed time-domain output (black) before de-emphasis, reconstruction error (red)') end; pause if sav==1 print -dbitmap c:\Temp\'Figure 5' end close; home; clc plot(in(1:rng),'k') hold on plot(in(1:rng)-out(1:rng),'r') if tit==1 title('Input signal (black), input-output (red)): pre-emphasised') end;pause if sav==1 print -dbitmap c:\Temp\'Figure 6' end close; home; clc
% time domain error check plot(ERR0(1:rng),'k') hold on plot(ERRCK(1:rng),'r') plot(clipl*ones(size(1:rng)),'k') plot(-clipl*ones(size(1:rng)),'k') if tit==1 title('Clipping error (black) at output of differentiator before correction, independent clipping check (red)') end; pause if sav==1 print -dbitmap c:\Temp\'Figure 7' end close; home; clc
% close program return
%***********************************************************************