1、%-testViterbi.m-% Annotated by Quan Xu on Nov. 22, 2011clear all;clc;close;delete *.mat;EsNodB = 0:1:10;EsNo = 10.(EsNodB/10);N_samples = 100;N_bits = 512;for k = 1 : length(EsNodB) for m = 1 : N_samples % noise power WGNPower = 1/EsNo(k); % for uncoded system WGNPower1 = 2/EsNo(k); % for 1/2 convol

2、utional coded system % random message msg = randint(1,N_bits);% unencode data b = msg; % Encode data bits trel = poly2trellis(3,7 5); % Define trellis b1 = convenc(msg,trel); % Encode % modulate x = 1 - 2*b;x1 = 1 - 2*b1; % AWGN channel h = ones(size(x); h1 = ones(size(x1); % Gaussian noise generati

3、on w = sqrt(WGNPower/2)*(randn(size(x) + sqrt(-1)*randn(size(x);w1 = sqrt(WGNPower1/2)*(randn(size(x1) + sqrt(-1)*randn(size(x1); % recieved signals y = h.*x + w;y1 = h1.*x1 + w1;% demodulate uncoded data outUncod = uncode(y,h); % Soft-decision viterbi decoding outVit1 = simeudistVit(trel, y1); errr

4、ate(m) = mean(abs(outUncod-msg); errrate1(m) = mean(abs(outVit1-msg); end % for m ber(k) = mean(errrate) ber1(k) = mean(errrate1) end % for ksemilogy(EsNodB,ber,'k',EsNodB,ber1,'-r');grid onaxis(1 10 10-5 1)xlabel('10log(Eb/No)')ylabel('BER')title('AWGN channel vi

5、terbi decoder');h = legend('uncoded', '2,1,2convolutional encoding');set(h,'Interpreter','none')save conViterbi.mat;% EOF%- testViterbi.m -%-simEudist.m -% Calulating Euclidean distance% (x-y).2 -> xy ; since x = +1 or -1; y is the same for every branch % Annot

6、ated by Quan Xu on Nov. 22, 2011function dist = simEudist(x,y)leth = length(x);realx = real(x);realy = real(y);dist = 0;for i=1:leth dist = dist + realx(i)*realy(i); end% EOF%-simEudist.m -%-uncode.m -% Annotated by Quan Xu on Nov. 22, 2011% demodulating the uncoded datafunction unOut = uncode(y,h)

7、for l = 1 : length(y) llr(l) = abs(y(l)+h(l)2 - abs(y(l)-h(l)2; x_hat(l) = sign(llr(l); b_hat(l) = (1-x_hat(l)/2;end % for lunOut = b_hat;% EOF%-uncode.m -%-trellis2enc.m -% Annotated by Quan Xu on Nov. 22, 2011function enc = trellis2enc( trl ),% put the trellis structure into a more user friendly m

8、anner enc.k = log2( trl.numInputSymbols ); % number of inputs enc.n = log2( trl.numOutputSymbols ); % numbor of outputs enc.r = enc.k / enc.n; % code rate enc.ksym = trl.numInputSymbols; % number of possible input combinations enc.nsym = trl.numOutputSymbols; % number of possible output combinations

9、 enc.stat = trl.numStates; % number of encoder states % forward transitions: enc.next.states = trl.nextStates + 1; % NEXT states enc.next.output = trl.outputs; % NEXT outputs for i = 1:enc.ksym, % NEXT (binary) outputs enc.next.binout( :,:,i ) = 1 - 2*de2bi( oct2dec( trl.outputs(:,i) ), enc.n, '

10、left-msb' ); end % store possible binary outputs and inputs: enc.inp = de2bi( oct2dec( 0:enc.ksym-1 ), enc.k, 'left-msb' ); % all possible binary inputs enc.out = de2bi( oct2dec( 0:enc.nsym-1 ), enc.n, 'left-msb' ); % all possible binary outputs enc.bininp = 2*enc.inp-1;return;%-

11、trellis2enc.m -%-convenc.m -% This function could be found in C:Program FilesMATLABR2008btoolboxcommcommconvenc.m% Annotated by Quan Xu on Nov. 22, 2011function varargout = convenc(msg, trellis, varargin)% Typical error checking.error(nargchk(2,4,nargin,'struct');nvarargin = nargin - 2;% Set

12、 defaultspunctVec = ;initialstate = 0;switch (nvarargin)case 1 if isempty(varargin1) if isscalar(varargin1) initialstate = varargin1; else punctVec = varargin1; end endcase 2 punctVec, initialstate = deal(varargin:);endif nargout > 2 error('comm:convenc:TooManyOutputArg','Too many out

13、put arguments.');end% check trellisif istrellis(trellis), error('comm:convenc:InvalidTrellis','Trellis is not valid.');end% Get info out of trellis structurek = log2(trellis.numInputSymbols);n = log2(trellis.numOutputSymbols);outputs = oct2dec(trellis.outputs);% Check msgif isemp

14、ty(msg) msg_dim = size(msg); if ( isnumeric(msg) | islogical(msg) ) | . length(msg_dim)>2 | . min(msg_dim)>1 error('comm:convenc:InvalidMsg','The input message must be a logical or numeric vector.'); end outLog = islogical(msg); % for output data type msg = double(msg); % for p

15、roper numerical operation if max(max(msg < 0) | . max(max(isfinite(msg) | . isreal(msg) | . max(max(floor(msg) = msg) | . max(max(msg) > 1 error('comm:convenc:InputNotBinary','The input message must contain only binary values.'); end if mod(length(msg), k) =0 error('comm:co

16、nvenc:InvalidMsgLength','Length of the input message must be a multiple of the ' . 'number of bits in an input symbol.'); end % Get message orientation if msg_dim(1)>1 msg_flip = 1; msg=msg' else msg_flip = 0; endend% Check Puncture vectorif isempty(punctVec) % Validity ch

17、eck if ( isnumeric(punctVec) | islogical(punctVec) ) | . length(size(punctVec) > 2 | . ( isvector(punctVec) && isscalar(punctVec) ) | . max(max(isfinite(punctVec) | . isreal(punctVec) error('comm:convenc:InvalidPuncPat', 'The puncture pattern parameter',. ' must be a v

18、ector of real or logical values.'); end % Binary value check if any(punctVec=0 & punctVec=1) error('comm:convenc:PuncPatNotBinary', . 'The puncture pattern parameter must be a binary vector of 1''s ',. 'and 0''s only.'); end % Length checks if length(p

19、unctVec) < n error('comm:convenc:InvalidPuncPatLength',. 'The puncture pattern parameter length must be at least the ',. 'number of bits in an output symbol.'); end if mod(length(msg)/k)*n, length(punctVec) =0 error('comm:convenc:InvalidCodeLengthPunc', . 'The

20、input message length divided by the base code rate must be an',. 'ninteger multiple of the length of the puncture pattern parameter.'); endend% Check initial stateif isnumeric(initialstate) | . isscalar(initialstate) | . max(max(initialstate < 0) | . max(max(isfinite(initialstate) | .

21、 isreal(initialstate) | . max(max(floor(initialstate) = initialstate) | . max(max(initialstate) > trellis.numStates-1 error('comm:convenc:InvalidInitialState','The initial state must be an integer scalar between 0 and ' . '(TRELLIS.numStates-1). See POLY2TRELLIS.');end% Re

22、turn if input message is emptyif isempty(msg) varargout1 = ; varargout2 = initialstate; return;end% Actual call to core function 'convcore.c'code, fstate = . convcore(msg,k,n,trellis.numStates,outputs,trellis.nextStates,initialstate);if isempty(punctVec) % Expand punctVec if needed if length

23、(code) = length(punctVec) pVec = punctVec(:); bb = pVec(:, ones(1, length(code)/length(punctVec); % repeat punctVec = bb(:); % vector end % Puncture the encoded output code(logical(punctVec) = ;end% Change code back to same orientation as input MSGif msg_flip code=code'end% Set output data type

24、to logical if appropriateif outLog, code = logical(code); end;% Set outputsvarargout = code, fstate;% EOF%-convenc.m -%-simeudistVit.m -% Using the simplified Euclidean distance% trellis2enc transfer the Matlab trellis to be user-friendly form% input arguments% trel : matlab trellis structure% y : r

25、eceived symbols from the channel% output % decVit : decoding bits% Annotated by Quan Xu on Nov. 22, 2011function decVit = simeudistVit(trel, y);enc = trellis2enc( trel );% some parametersk = enc.k;n = enc.n;r = enc.r;ksym = enc.ksym;nsym = enc.nsym;numStat = enc.stat;states = enc.next.states;outputs

26、 = enc.next.binout;leth = length(y)/n;%Input table; for decoding with current state and previous statedecodBit = 2*ones(numStat); %Row numbers represent current states plus 1 while columns represent next states plus 1for i=1:numStat for j=1:2 nxt = states(i,j); decodBit(i,nxt) = j-1; end %jend %i%-c

27、aculating path metric array &surviving predecessor states array initial valuenumMem = log2(numStat); % Memories of the encoder: number of branches before constructing four statespm = zeros(1,numStat);pmtmp = pm;surpred = zeros(numStat,leth);inipred = 1 zeros(1,numStat-1);initmp = inipred;for i =

28、 1:numMem iy = zeros(1,n); for p = 1:n iy(p) = y(i-1)*n+p); end for j = 1:numStat if initmp(j) = 1 ininxt0 = states(j,1); pm(ininxt0) = pmtmp(j)+simEudist(outputs(j,:,1),iy); surpred(ininxt0,i) = j; inipred(ininxt0) = 1; ininxt1 = states(j,2); pm(ininxt1) = pmtmp(j)+simEudist(outputs(j,:,2),iy); surpred(ininxt1,i) = j; inipred(ininxt1) =1; end end %j initmp = inipred; pmtmp = pm;end %i%-caculating path metric array &surviving predecessor states.%.array initial value% indiceid = zeros(2*numStat, 2);for i = 1:2*numStat id(i,1) = mod(i-1,numStat)+1; id(i,2) = floor(i-1)/numStat)+1; en