Implementation and Verification of Waveform Coding Using MATLAB with Applications to Speech Compression
OBJECTIVES
2. To implement and experience the 
-law PCM coding
3. To implement and experience ADPCM coding
4. To implement and experience DCT and W-DCT coding
PROCEDURE:
Part A: Pulse Code Modulation (PCM Coding)
Given the data file “speech.dat” with 16 bits per sample and a sampling rate of 8 kHz, develop the PCM codec (see Chapter 11) and perform compression and decompression, and complete the following table.
Examine the given speech and determine
File size = _________________ (bytes)
The sound quality will be rated as the following: excellent, good, intelligent, unacceptable.
>>sound(speech/max(abs(speech)), 8000); % listen to the normalized original speech
>>pause
>>sound(quspeech/max(abs(quspeech)),8000); %listen to the normalized decoded speech
| PCM with Midtread Quantizer | 4 bits/sample | 6 bits/sample | 8 bits/sample |
| File size (bytes) | | | |
| Compression ratio | | | |
| Signal-to-noise ratio | | | |
| Sound quality | | | |
Plot the original speech, decompressed speech, and quantized error for each case using the MATLAB subplot.
Part B: -law PCM Coding
Given the data file “speech.dat” with 16 bits per sample and a sampling rate of 8 kHz, develop the -law compression codec (see Chapter 11) and perform compression and decompression, and complete the following table.
Examine the given speech and determine
File size = _________________ (bytes)
The sound quality will be rated as the following: excellent, good, intelligent, unacceptable.
|
| 4 bits/sample | 6 bits/sample | 8 bits/sample |
| File size (bytes) | | | |
| Compression ratio | | | |
| Signal-to-noise ratio | | | |
| Sound quality | | | |
Plot the original speech, decompressed speech, and quantized error for each case using the MATLAB subplots.
Part C: Adaptive Differential PCM (ADPCM Coding)
Given the data file “speech.dat” with 16 bits per sample, a sampling rate of 8 kHz, and ADCPM codec (see Chapter 11), perform compression and decompression, and complete the following table
Examine the given speech and determine
File size = _________________ (bytes)
The sound quality will be rated as the following: excellent, good, intelligent, unacceptable.
| ADPCM codec | 4 bits/sample |
| File size (bytes) | |
| Compression ratio | |
| Signal-to-noise ratio | |
| Sound quality | |
Plot the original speech, decompressed speech, and quantized error using the MATLAB subplot.
Part D: Discrete-cosine Transform (DCT) Coding and M-DCT (Modified DCT) Coding
Given the data file “speech.dat” with 16 bits per sample, a sampling rate of 8 kHz, and DCT and MDCT codec in Chapter 11, perform compression and decompression using the following specified parameters and complete the following table.
Speech data: 16 bit per sample
Sampling rate: 8000 Hz
Block size: 16 samples
Scale factor: see Table below
Coefficients: 3 or 4-bit linear quantizer (see Table below)
Compare the sound quality.
Examine the given speech and determine
File size = _________________ (bytes)
Plot and listen to the original speech, DCT decompressed speech, and M-DCT decompressed speech, respectively.
| Number of bits nbits = | Scale factor scalef= | DCT sound quality | M-DCT sound quality |
| 3 | scalef1bit | | |
| 3 | scalel2bits | | |
| 4 | scalef1bit | | |
| 4 | scalef2bits | | |
Comment on your results
Some keys for Lab 9
close all; clear all
disp('load speech data');
load speech.dat;
lg=length(speech);
t=[0:1:lg-1]/8000;
disp('loading finished');
disp('mulaw companding')
nspeech=speech/(2^15); % 15 bits
mu=input('input mu =>');
for x=1:lg
munspeech(x) =mulaw(nspeech(x),1,mu); %mulaw compression
end
disp('finished mu-law companding');
disp('start to quantization')
bits = input('input bits=>');
for x=1:lg
[pq uindx(x)]= midtread(bits,1,nspeech(x));
[pq muindx(x)]= midtread(bits,1,munspeech(x));
end
%%
% transmission
%
disp('expander');
for x=1:lg
qunspeech(x)= mtrdec(bits,1,uindx(x));
qmunspeech(x)=mtrdec(bits,1,muindx(x));
end
for x=1:lg
expnspeech(x)= muexpand(qmunspeech(x),1,mu);
end
quspeech=qunspeech.*2^15;
qspeech =expnspeech.*2^15;
disp('finished')
qerr = speech-qspeech;
subplot(2,1,1),plot(t, speech, 'w', t, qspeech, 'c',t,qspeech-speech,'r');grid
subplot(2,1,2),plot(t, speech, 'w', t, quspeech,'b',t,quspeech-speech,'r');grid
disp('speech:orginal data 15 bits');
disp('quspeech: PCM in quantized');
disp('qspeech: mulow deccoded');
disp('SNR speech and qspeech');
snr(speech,qspeech);
disp('SNR speech quspeech');
snr(speech,quspeech);
function qvalue = mulaw(vin, vmax, mu)
vin = vin/vmax;
qvalue = vmax*sign(vin)*log(1+mu*abs(vin))/log(1+mu);
function rvalue = muexpand(y,vmax, mu)
y=y/vmax;
rvalue=sign(y)*(vmax/mu)*((1+mu)^abs(y) -1);
function [ pq, indx ] = midtread(NoBits,Xmax,value)
% function [pq indx] = midtread(NoBits, Xmax, value)
% this routine is created for simulation of uniform quatizer.
%
% NoBits: number of bits used in quantization.
% Xmax: overload value.
% value: input to be quantized.
% pq: output of quantized value
% indx: codeword integer
% Note: the midtread method is used in this quantizer.
%
if NoBits == 0
pq = 0;
indx=0;
else
delta = 2*abs(Xmax)/(2^NoBits-1);
Xrmax=delta*(2^NoBits/2-1);
if abs(value) >= Xrmax
tmp = Xrmax;
indx=(2^NoBits/2-1);
else
tmp = abs(value);
end
indx=round(tmp/delta);
pq =round(tmp/delta)*delta;
if value <>
pq = -pq;
indx=-indx;
end
end
function pq = mtrdec(NoBits,Xmax,indx)
% function pq = mtrdec(NoBits, Xmax, indx)
% this routine is created for simulation of uniform quatizer.
%
% NoBits: number of bits used in quantization.
% Xmax: overload value
% pq: output of quantized value
% indx: codeword integer
% Note: the midtread method is used in this quantizer.
%
if NoBits == 0
pq = 0;
else
delta = 2*abs(Xmax)/(2^NoBits-1);
pq=indx*delta;
end
function snr = calcsnr(speech, qspeech)
% function snr = calcsnr(speech, qspeech)
% this routine is created for calculation of SNR
%
% speech: original speech waveform.
% qspeech: quantized speech.
% snr: output SNR in dB.
%
% Note: midrise method is used in this quantizer.
%
qerr = speech-qspeech;
snr = 10*log10(sum(speech.*speech)/sum(qerr.*qerr))
% Waveform coding using DCT and MDCT for a block size of 16 samples
% main program
close all; clear all
load speech.dat % provided by the instructor
% create scale factors
N=16; % block size
scalef4bits=sqrt(2*N)*[1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768];
scalef3bits=sqrt(2*N)*[256 512 1024 2048 4096 8192 16384 32768];
scalef2bits=sqrt(2*N)*[4096 8192 16384 32768];
scalef1bit=sqrt(2*N)*[16384 32768];
scalef=scalef2bits;
nbits =3;
% ensure the block size to be 16 samples.
x=[speech zeros(1,16-mod(length(speech),16))];
Nblock=length(x)/16;
DCT_code=[]; scale_code=[];
% encoder
for i=1:Nblock
xblock_DCT=dct(x((i-1)*16+1:i*16));
diff=(scalef-(max(abs(xblock_DCT))));
iscale(i)=min(find(diff=min(diff(find(diff>=0))))); %find a scale factor
xblock_DCT=xblock_DCT/scalef(iscale(i)); % scale the input vector
for j=1:16
[DCT_coeff(j) pp]=biquant(nbits,-1,1,xblock_DCT(j));
end
DCT_code=[DCT_code DCT_coeff ];
end
%decoder
Nblock=length(DCT_code)/16;
xx=[];
for i=1:Nblock
DCT_coefR=DCT_code((i-1)*16+1:i*16);
for j=1:16
xrblock_DCT(j)=biqtdec(nbits,-1,1,DCT_coefR(j));
end
xrblock=idct(xrblock_DCT.*scalef(iscale(i)));
xx=[xx xrblock];
end
% Transform coding using MDCT
xm=[zeros(1,8) speech zeros(1,8-mod(length(speech),8)), zeros(1,8)];
Nsubblock=length(x)/8;
MDCT_code=[];
% encoder
for i=1:Nsubblock
xsubblock_DCT=wmdct(xm((i-1)*8+1:(i+1)*8));
diff=(scalef-max(abs(xsubblock_DCT)));
iscale(i)= iscale(i)=min(find(diff=min(diff(find(diff>=0))))); %find a scale factor
xsubblock_DCT=xsubblock_DCT/scalef(iscale(i)); % scale the input vector
for j=1:8
[MDCT_coeff(j) pp]=biquant(nbits,-1,1,xsubblock_DCT(j));
end
MDCT_code=[MDCT_code MDCT_coeff];
end
%decoder
% recover thr first subblock
Nsubblock=length(MDCT_code)/8;
xxm=[];
MDCT_coeffR=MDCT_code(1:8);
for j=1:8
xmrblock_DCT(j)=biqtdec(nbits,-1,1,MDCT_coeffR(j));
end
xmrblock=wimdct(xmrblock_DCT*scalef(iscale(1)));
xxr_pre=xmrblock(9:16) % recovered first block for overlap and add
for i=2:Nsubblock
MDCT_coeffR=MDCT_code((i-1)*8+1:i*8);
for j=1:8
xmrblock_DCT(j)=biqtdec(nbits,-1,1,MDCT_coeffR(j));
end
xmrblock=wimdct(xmrblock_DCT*scalef(iscale(i)));
xxr_cur=xxr_pre+xmrblock(1:8); % overlap and add
xxm=[xxm xxr_cur];
xxr_pre=xmrblock(9:16); % set for the next overlap
end
subplot(3,1,1);plot(x,'k');grid; axis([0 length(x) -10000 10000])
ylabel('Original signal');
subplot(3,1,2);plot(xx,'k');grid;axis([0 length(xx) -10000 10000]);
ylabel('DCT coding')
subplot(3,1,3);plot(xxm,'k');grid;axis([0 length(xxm) -10000 10000]);
ylabel('W-MDCT coding');
xlabel('Sample number');
function [ tdac_coef ] = wmdct(ipsig)
%
% This function transforms the signal vector using the W-MDCT
% usage:
% ipsig: inpput signal block of N samples (N=even number)
% tdac_coe: W-MDCT coefficents (N/2 coefficients)
%
N = length(ipsig);
NN =N;
for i=1:NN
h(i) = sin((pi/NN)*(i-1+0.5));
end
for k=1:N/2
tdac_coef(k) = 0.0;
for n=1:N
tdac_coef(k) = tdac_coef(k) + ...
h(n)*ipsig(n)*cos((2*pi/N)*(k-1+0.5)*(n-1+0.5+N/4));
end
end
tdac_coef=2*tdac_coef;
function [ opsig ] = wimdct(tdac_coef)
%
% This function transform the W-MDCT coefficients back to the signal
% usage:
% tdac_coeff: N/2 W-MDCT coeffcients
% opsig: output signal black with N samples
%
N = length(tdac_coef);
tmp_coef = ((-1)^(N+1))*tdac_coef(N:-1:1);
tdac_coef = [ tdac_coef tmp_coef];
N = length(tdac_coef);
NN =N;
for i=1:NN
f(i) = sin((pi/NN)*(i-1+0.5));
end
for n=1:N
opsig(n) = 0.0;
for k=1:N
opsig(n) = opsig(n) + ...
tdac_coef(k)*cos((2*pi/N)*(k-1+0.5)*(n-1+0.5+N/4));
end
opsig(n) = opsig(n)*f(n)/N;
end
Posts
0 comments:
Post a Comment