6.5  DMT / OFDM in Matrix Notation

AK-TODO: DMT / OFDM in Matrix Notation

Botar em algum lugar: o canal eh linear, queremos usar FFT, dribla-se a convolucao circular da FFT com CP...

Extended Listing E.48 to OFDM is Listing 6.4.

Listing 6.4: MatlabOctaveCodeSnippets/snip_multicarrier_DMTInMatrixNotation.m
1%% This code demonstrates DMT using convolution in matrix notation. 
2debugMe=0; %enable for debugging 
3useSimplifiedConvolution=1; %if enabled, convolution is calculated 
4%only for valid samples. Otherwise, the complete conv matrix is used 
5% 
6%The samples contained in the vectors below are ordered from the 
7%least recent to the most recent, the samples x[0], x[1], x[2] 
8%become the column vector [x[0];x[1];x[2]]. 
9%For simplicity, all tones transmit the same QAM constellation. 
10M=32; %modulation order, number of symbols in alphabet 
11Nfft=16; %must be even 
12numDMTSymbols=50; %number of DMT symbols to be transmitted 
13h=[8.4 -2 -2.8 2 -2.8 2 -2.8 -2]; %longer example 
14%h=[1 0.9 0.8]; %impulse response using typical h=[h[0] h[1] h[2]] 
15if debugMe==1 
16    M=4; %modulation order, number of symbols in alphabet 
17    Nfft=6; %must be even, 6 will lead to only 2 tones carrying info. 
18    numDMTSymbols=5; %number of DMT symbols to be transmitted 
19    h=[1 0.9 0.8]; %impulse response using typical h=[h[0] h[1] h[2]] 
20end 
21%% Pre-processing / initialization 
22constellation=ak_qamSquareConstellation(M); %QAM constellation 
23numOfTones=(Nfft-2)/2; %chose not to use tones at DC and Nyquist 
24D=length(h)-1; %channel dispersion 
25CP=D; %cyclic prefix 
26if Nfft <= CP 
27    error('Nfft must be increased'); 
28end 
29hmatrix = convmtx(transpose(h),Nfft+CP); 
30symbolIndicesTx=1+floor(M*rand(numOfTones,numDMTSymbols)); %randomize 
31if debugMe==1 %create indices 1,2,3...M,1,2,3...,M up to the end 
32    tempInteger = floor(numOfTones*numDMTSymbols/M); 
33    tempRemainder = numOfTones*numDMTSymbols-tempInteger*M; 
34    symbolIndicesTx = [repmat(1:M,1,tempInteger) 1:tempRemainder]; 
35    %symbolIndicesTx=reshape(tempSequence,numOfTones,numDMTSymbols); 
36end 
37symbolsTx=constellation(symbolIndicesTx); %random QAM symbols 
38X=reshape(symbolsTx,numOfTones,numDMTSymbols); %QAM symbols as matrix 
39Xk=[zeros(1,numDMTSymbols) %DC level is assumed zero 
40    X 
41    zeros(1,numDMTSymbols) %Nyquist level is assumed zero 
42    conj(flipud(X))]; %Xk has Hermitian symmetry 
43x=ifft(Xk)*sqrt(Nfft); %adopt orthonormal FFT 
44P=[zeros(CP,Nfft-CP) eye(CP,CP) %permutation matrix 
45    eye(Nfft,Nfft)]; 
46xp=P*x; %add cyclic prefix using permutation matrix 
47if useSimplifiedConvolution==1 
48    %take out CP and convolution tail. The convolution tail of 
49    %current block will affect only the CP samples of the next block. 
50    %Because the CP samples are discarded, both CP samples and tails 
51    %can be ignored unless the complete convolution is desired. 
52    hmatrixValid = hmatrix(CP+1:CP+Nfft,:); 
53    yvalid=hmatrixValid*xp; %convolve with the channel 
54else 
55    y=hmatrix*xp; %convolve with the channel 
56    yvalid = y(CP+1:CP+Nfft,:); %extract only valid samples 
57end 
58if debugMe==1 
59    A=flipud(xp); 
60    xp2=A(:); 
61    y2=conv(xp2,h); 
62    y2=y2(1:(Nfft+CP)*numDMTSymbols); 
63    y2valid=reshape(y2,Nfft+CP,numDMTSymbols); 
64    y2valid(1:CP,:)=[]; %take out CP 
65    discrepancy=yvalid - flipud(y2valid); 
66    plot(discrepancy(:)); 
67end 
68%% Receiver 
69H = fft(h,Nfft); %channel frequency response 
70H = H(2:Nfft/2); %only positive frequencies, no DC nor Nyquist 
71H = H(:); %make it a column vector 
72G = 1./H; %design the frequency equalizer (FEQ) 
73Y=fft(yvalid)/sqrt(Nfft); %an orthonormal FFT has been adopted 
74Y=Y(2:Nfft/2,:); %only positive frequencies, without DC and Nyquist 
75Yeq = repmat(G,1,numDMTSymbols).*Y; %perform equalization 
76[symbolIndicesRx, symbolsRx]=ak_qamdemod(Yeq,M); %demodulate 
77symbolIndicesRx=symbolIndicesRx+1; %make first index to be 1, not 0 
78if debugMe==1 %plot constellations 
79    clf, plot(real(symbolsTx),imag(symbolsTx),'xb'); hold on 
80    plot(real(Yeq),imag(Yeq),'or'); 
81end 
82disp(['SER=' num2str(mean(symbolIndicesTx(:)~=symbolIndicesRx(:)))])