Podstawy szerokopasmowej transmisji przewodowej
Autorzy: Jarosław Bułat
Pliki: rozdzial_23.zip
przykład 23.1
W przykładzie 23.1 (plik przyklad_23_1.m) przedstawiono ilustrację transmisji według schematu z rysunku 23.4. Jest to prosty symulator transmisji DMT, zestrojony dla usługi ADSL. Uwzględnia on modulacje i demodulację QAM, kanał transmisyjny w postaci filtra FIR, korektor czasowy TEQ, korektor częstotliwościowy FEQ, blok funkcyjny OFDM oraz zakłócenia AWGN. Program oblicza średni SNR w każdym podpaśmie częstotliwościowym dla przesłanych Ns ramek sygnału. Do poprawnego działania należy dostarczyć odpowiedź impulsową kanału komunikacyjnego – c, odpowiedź impulsową filtra TEQ – w oraz opóźnienie Dx, skojarzone ze skróconą odpowiedzią impulsową. Symulacja jest wykonywana poprzez zakodowanie sygnału pseudolosowego, przesłanie go w dziedzinie czasu, zdekodowanie i następnie porównanie z wysłanym sygnałem. Założono, że odbiornik jest zsynchronizowany z nadajnikiem.% Przykład 23.1 - ilustracja transmisji ADSL zgodnie ze schematem przedstawionym
% na rysunku 23.4
clear;
fs = 2.208e6;
codingGain = 4.2;
margin = 6;
% inicjalizacja parametrów transmisji
N = 512;
Ld = 32;
Le = 16;
Nd = N + Ld;
Ns = 200;
inputPower = 23;
inputImp = 100;
lineNoise = -120;
load c_w_Dx;
% korektor TEQ, IR kanału transmisyjnego
c = [ c zeros( 1, N-length(c)) ];
sir = conv( w, c );
% skrócona odpowiedz impulsowa kanału SIR
sir = [ sir zeros(1,N-length(sir)) ];
SIR = fft( sir(1:N ) ).'; k = (0:N-1)';
FEQ = SIR .* (1/(2*sqrt(2))) .* exp( j*2*pi/N *k* Dx );
% wagi korektora FEQ
rb = zeros(2*Nd,1);
tf_sig = zeros(N,1);
tf_err = zeros(N,1); % inicjalizacja modemu
U = zeros(N,1);
U0 = zeros(N,1);
c_res = zeros( 1, length( c ) -1 );
w_res = zeros( 1, length( w ) -1 );
noiseAWGN = randn(Ns*Nd,1) * sqrt(inputImp*0.001*fs/2*10^(lineNoise/10)); % szumy
signalPower = inputImp * 0.001 * 10^( inputPower/10 ); % wzmocnienie sygnału
scalingFactor = sqrt( 10^( codingGain/10 ) ) * sqrt( signalPower / 1 );
for iter=1:Ns
U0 = U;
Ur=sqrt(3)*(2*rand(N/2-1, 1)-1);
Ui=sqrt(3)*(2*rand(N/2-1, 1)-1);
U = Ur+j*Ui;
U = [ 0; U; 0; conj( flipud(U) ) ];
x = ifft( U );
x = ( real( x )*sqrt(N) ) / sqrt(2);
x = [(x(N-Ld+1:N)); x];
s = x * scalingFactor;
[ y c_res ] = filter( c, 1, s, c_res ); % kanał
y = y + noiseAWGN( (iter-1)*Nd+1:iter*Nd ); % dodanie szumu addytywnego
y = y / scalingFactor; % normowanie mocy sygnału
[ ye w_res ] = filter( w, 1, y, w_res ); % korektor TEQ
rb(1:Nd) = rb(Nd+1:2*Nd);
rb(Nd+1:2*Nd) = ye;
y = rb(1+Dx+Ld : N+Dx+Ld); % usunięcie prefiksu
Y= fft(y); % blok funkcyjny OFDM
Y = Y/(2*sqrt(N)); % normowanie sygnału
U_r = Y ./ FEQ; % korektor FEQ
if( iter > 2 ) % wyznaczanie SNR
e = U0 - U_r; % błąd rekonstrukcji
tf_sig = tf_sig + U0.*conj(U0); % moc sygnału
ce = e .* conj(e);
tf_err = tf_err + ce; % moc zakłóceń
end
end
idx = find( tf_err == 0 ); tf_err( idx ) = eps;
snr = tf_sig(1:N/2) ./ tf_err(1:N/2); % stosunek mocy sygnału do szumu
Przykład 23.2
Przykład 23.2 (plik adsl.m) prezentuje implementację w języku MATLAB kodu źródłowego głównej pętli zbudowanego symulatora ADSL. Sposób użycia tej funkcji oraz inicjalizacja niezbędnych parametrów zostaną podane w przykładzie 23.3.% Przykład 23.2 - plik adsl.m
function y = adsl( iter, trType, bitMask, teqIR, Dx, feq, noiseType )
% iter: liczba iteracji (ramek)
% trType: typ transmisji:
% 0 = sygnał treningowy {1,2,3}, bez prefiksu, bez TEQ, bez FEQ, bitMask ignorowany
% 1 = estymacja SNR: jak dla 0 plus prefiks, TEQ, FEQ, bitMask ignorowany
% 2 = transmisja danych: wszystkie elementy modemu działają
% bitMask: (N/2,1) liczba bitów transmitowany w podkanale, 0 - brak transmisji
Lp = 32; % długość prefiksu
N = 512; % długość ramki bez prefiksu
fp = 2.208e6; % częstotliwość próbkowania
sigPower = 23; % moc sygnału [dBm]
% noiseType = 'AWGN'; % %typ zakłóceń addytywnych (szumu): AWGN, NBI, impulsowy, etc...
noisePower = -100; % moc szumu AWGN [dBm/Hz]
noisePowerNBI = -20; % moc szumu NBI [dBm]
noiseChannelNBI = 100; % numer podkanału w którym wystąpiło zakłócenie NBI
channelNr = 2; % numer kanału z bazy csaloop: do wyboru {1...8, 10, 11}
H = zeros( 1, N/2 );
sig = zeros( 1, N/2 );
noise = zeros( 1, N/2 );
ber = zeros( 1, N/2 );
% główna pętla programu - jedna iteracja = transmisja jednej ramki danych
for i = 1:iter
X1 = genBinData( trType, bitMask );
X2 = QAM( X1, trType, bitMask );
x3 = Mod( X2, N );
x4 = Prefix( x3, trType, Lp );
x5 = NormEnergy( x4, sigPower );
y1 = channel( x5, channelNr );
y2 = AddNoise( y1, noiseType, noisePower, noisePowerNBI, noiseChannelNBI, fp, iter );
y3 = DeNormEnergy( y2, sigPower );
y4 = TEQ( y3, teqIR );
y5 = DePrefix( y4, trType, N, Dx, Lp );
Y6 = DeMod( y5, N );
Y7 = FEQ( Y6, feq );
Y8 = DeQAM( Y7, trType, bitMask );
H = HEstim( i, trType, X2, Y6, H );
[sig, noise] = SNR( i, trType, X2, Y7, sig, noise );
ber = BER( i, trType, bitMask, X1, Y8, ber );
% err = cmp( X1, Y8, X2, Y7 );
end
if( trType == 0 ) % estymacja knanału
H( N/2+2:N ) = conj( fliplr( H(2:N/2) ) );
H = H / (iter-1);
h = ifft( H );
h = real(h);
y = h;
elseif( trType == 1 ) % obliczanie SNR
snr = sig./noise;
y = snr;
elseif( trType == 2 ) % transmisja bitów
y = ber;
end
%===================== genBinData ====================
function y = genBinData( trType, mask )
% trType: typ transmisji:
% 0 = sygnał treningowy {1,2,3,4}, bez: prefiksu, TEQ, FEQ; bitMask ignorowany
% 1 = esytmacja SNR: wygnał jak w 0, uwzględniono: prefiks, TEQ, FEQ, bitMask ignorowany
% 2 = transmisja danych
% mask(N/2,1): ilość bitół w każdy podpaśmie
% y(N/2,1): wektor wyjściowy zawiarający losowe dane
persistent state;
if isempty( state ) % zapewnia brak korelacji pomiędzy szumem danymi (losowymi)
rand( 'state', 111 );
disp( 'initialize persistent random data generator' );
else
rand( 'state', state );
end
if( trType == 0 || trType == 1 ) % ReverbSequence (ADSL norma C-REVERB1 s.92)
N = 512;
b = zeros( N, 1 );
b( 1:9 ) = 1;
for i =10:length(b)
b(i) = xor( b(i-4), b(i-9) );
end
for i = 1:N/2
y( i ) = ( 2 * b( 2*(i-1)+1 ) + b( 2*(i-1)+2 ) ) +1;
end
elseif( trType == 2 )
for s=1:length(mask)
if( mask(s) == 0 )
y(s) = 0;
else
y(s) = randint( 1, 1, 2^mask(s) );
end
end
else
y = zeros( size(mask) );
end
state = rand( 'state' );
%===================== QAM ====================
function y = QAM( x, trType, mask )
% x(N/2,1): dane binarne {0, 1, 2, 3, 4, 5, ... 2^n-1}
% mask(N/2,1): liczba bitów w każdym podkanale
% y(N/2,1): QAM
y=zeros( size(mask) );
if( trType == 0 || trType == 1 ) % estymacja SNR i kanału, używane wyłącznie QAM4
QAM4 = [ 1+j, 1-j; -1+j, -1-j ];
y = zeros( size( mask ) );
for i = 1:length( mask )
if( mask(i) ~= 2 )
y(i) = 0;
else
y(i) = QAM4( x(i) );
end
end
elseif( trType == 2 ) % QAM4, QAM8, ..., QAM32768
qamMax = [ 0 1 3 3 5 7 11 15 23 31 47 63 95 127 191 ];
for s=1:length(mask)
M = 2^mask(s);
if( mask( s ) == 0 )
y(s) = 0;
else
y(s) = qammod( x(s), M ) / qamMax( mask(s) );
end
end
else
y = zeros( size(mask) );
end
%===================== Mod ====================
function y = Mod( x, N )
% x(N/2,1): QAM
% N: długość ramki (512 dla adsl)
% y(N,1): dane zmodulowane (ifft dla adsl)
% uwaga: zmienna y nie może zawierać części urojonej
% sygnał wyjściowy musi być unormowany do energii jednostkowej
x(1) = sqrt( 2 );
x( N/2+1 ) = sqrt( 2 );
% x(1) = 0;
% x( N/2+1 ) = 0;
x( N/2+2: N ) = conj( fliplr( x(2:N/2 ) ) );
y = ifft( x );
y = y.*sqrt(N/2);
%===================== prefix ====================
function y = Prefix( x, trType, Lp )
% x(N,1): dane zmodulowane (w dziedzinie czasu)
% trType: typ transmisji:
% 0: bez: prefiksu, TEQ, FEQ; bitMask ignorowany
% 1,2: transmisja, wymagany przefiks
if( trType == 0 )
y = x;
else
y = [ x(end-Lp+1:end), x ];
end
%===================== NormEnergy ===============
function y = NormEnergy( x, power )
% x(N+Lp,1): jedna ramka danych w dziedzinie czasu
% power: moc syngnału w dBm
inputImp = 100;
signalPower = inputImp * 0.001 * 10^( power/10 );
y = x .* sqrt( signalPower );
%===================== channel ==================
function y = channel( x, nr )
% x(N+Lp,1): jedna ramka sygnału przed kanałem
% nr: numer symulowanego kanału {1...8, 10, 11}
% y(N+Lp,1): jedna ramka sygnału po przejściu przez kanał
persistent ir_residue;
persistent ir;
if isempty( ir_residue )
load ir; % ładowanie odpowiedzi impulsowych różnych kanałów
eval( sprintf( 'ir=ir%d;', nr ) );
ir_residue = zeros( length( ir ) -1, 1 );
disp( 'initialize persistent ir_residue array' );
end
[ y ir_residue ] = filter( ir, 1, x, ir_residue );
%===================== noise ===================
function y = AddNoise( x, type, power, powerNBI, channelNBI, fs, iter )
% x(N+Lp,1)
% type: typ szumu: {AWGN, NBI, Impulse}
% power: moc szumu w dBm/Hz
% fp: częstotliwość pró”kowania
persistent state;
persistent NBI;
if isempty( state ) % zapewnia brak korelacji pomiędzy szumem danymi (losowymi)
rand( 'state', 112 );
disp( 'initialize persistent random noise' );
NBI = narrowNoiseGen( fs, iter*544, channelNBI/256/2*fs, 'AM-FIR', fs/512, 12345 )';
else
rand( 'state', state );
end
y = x;
if( strcmp(type,'AWGN') || strcmp( type, 'AWGN+NBI' ) )
inputImp = 100;
npower = inputImp * 0.001 * fs/2 * 10^( power/10 );
xn = randn( size(x) );
y = y + xn .* sqrt( npower );
end
if( strcmp( type, 'NBI' ) || strcmp ( type, 'AWGN+NBI' ) )
inputImp = 100;
noiseNotchGain = inputImp * 0.001 * 10^( powerNBI/10 );
y = y + NBI( 1+(iter-1)*length(x) : iter*length(x) ) * sqrt( noiseNotchGain );
end
if( strcmp( type, 'Impulse' ) )
error( 'AddNoise - Impulse noise: missing code' );
end
state = rand( 'state' );
%===================== NormEnergy ===============
function y = DeNormEnergy( x, power )
% x(N+Lp,1)
% y(N+Lp,1)
inputImp = 100;
signalPower = inputImp * 0.001 * 10^( power/10 );
y = x ./ sqrt( signalPower );
%===================== TEQ ====================
function y = TEQ( x, ir )
% x(N+Lp,1)
% ir(Le,1): odpowiedź impulsowa korektora TEQ
% y(N+Lp,1): ramka poddana korekcji TEQ w dziedzinie czasu
persistent eq_residue;
if isempty( eq_residue )
eq_residue = zeros( length( ir ) -1, 1 );
disp( 'initialize persistent eq_residue array' );
end
[ y eq_residue ] = filter( ir, 1, x, eq_residue );
%===================== DePrefix =================
function y = DePrefix( x, trType, N, D, Lp )
% x(N+Lp,1)
% type: typ transmisji
% 0: sekwencja treningowa - bez prefiksu
% 1: transmisjia - prefiks wymagany
% N: liczba próbek w ramce (długość ramki w dziedzinie czasu)
% D: opóźnienie transmisji powodowane przez kanał i TEQ
% Lp: długość prefiksu
% y(N,1): jedna ramka danych w dziedzinie czasu bez prefiksu
persistent lastFrame;
if isempty( lastFrame )
lastFrame = zeros( size( x ) );
disp( 'initialize last frame delay' );
end
cx = [ lastFrame, x ];
lastFrame = x;
if( trType == 0 )
y = x;
else
y = cx( D+Lp+1 : N+D+Lp );
end
%===================== DeMod ==================
function y = DeMod( x, length )
% x(N,1)
% length: długość ramki w dziedzinie czasu
% y(N/2,1): zespolone dane w dziedzinie częstotliwości (zdemodulowane)
y = fft( x );
y = y( 1:length/2 );
y = y ./ sqrt( length/2 );
%===================== FEQ ====================
function y = FEQ( x, feq )
% x(N/2,1): dane zespolone
% feq(N/2,1): korektor częstotliwościowy - FEQ
y = x .* feq;
%===================== DeQAM ==================
function y = DeQAM( x, trType, mask )
% x(N/2,1): dane po korekcji FEQ - dziedzina częstotliwości
% mask(N/2,1): liczba bitów w podkanałach
% y(N/2,1): dane binarne
if( trType == 0 || trType == 1 )
%QAM4 = [ 1+j, 1-j; -1+j, -1-j ];
x1 = round( x );
for i = 1:length( x )
rx = real( x1(i) );
ix = imag( x1(i) );
if( rx==1 && ix==1 )
y(i) = 1;
elseif( rx==1 && ix==-1 )
y(i) = 2;
elseif( rx==-1 && ix==1 )
y(i) = 3;
elseif( rx==-1 && ix==-1 )
y(i) = 4;
else
y(i) = 0;
end
end
elseif( trType == 2 )
qamMax = [ 0 1 3 3 5 7 11 15 23 31 47 63 95 127 191 ];
for s=1:length(mask)
M = 2^mask(s);
if( mask( s ) == 0 )
y(s) = 0;
else
y(s) = qamdemod( x(s)*qamMax( mask(s) ), M );
end
end
else
y=zeros( size(mask) );
end
%===================== H estimation ===============
function H = HEstim( i, trType, X, Y, Hp )
% estymacja odpowiedzi impulsowej, pierwsza ramka jest błędna, należy ją zignorować
% iter: numer iteracji (przesłanej ramki)
% type: typ transmisji
% 0: sekwencja treningowa = identyfikacja kanału
% 1: estymacja SNR
% 2: transmisjia bitów
% X(N/2,1): transmitowane dane w dziedzinie częstotliwości
% Y(N/2,1): odebrane dane w dziedzinie częstotliwości
% Hp(N/2,1): odpowiedź częstotliwościowa kanału z poprzedniej ramki
% H(N/2,1): odpowiedź częstotliwościowa kanału - uśrednianie
if( trType ~= 0 )
H = Hp;
return;
else
if( i == 1 )
H = zeros( size(Hp) );
else
H = Hp + Y./X;
end
end
function [sig, noise] = SNR( i, trType, X, Y, sigPrev, noisePrev )
% estymacja SNR (Signal-to-Noise-Ratio), pierwsza ramka jest błędna, należy ją zignorować
% iter: iteration number
% iter: numer iteracji (przesłanej ramki)
% type: typ transmisji
% 0: sekwencja treningowa = identyfikacja kanału
% 1: estymacja SNR
% 2: transmisjia bitów
% X(N/2,1): transmitowane dane w dziedzinie częstotliwości
% Y(N/2,1): odebrane dane w dziedzinie częstotliwości
% sigPrev: moc sygnału z poprzedniej ramki
% noisePrev: moc szumół z poprzedniej ramki
% sig: moc sygnału - uśrednianie
% noise: moc szumu - uśrednianie
if( trType ~=1 )
sig = zeros( size( sigPrev ) );
noise = zeros( size( noisePrev ) );
return;
else
if( i == 1) % pierwsza ramka
sig = zeros( size( sigPrev ) );
noise = zeros( size( noisePrev ) );
else
e = Y - X;
e = e.*conj(e);
sig = sigPrev + (X.*conj(X));
noise = noisePrev + e;
end
end
function ber = BER( i, trType, mask, X, Y, berp )
% estymacja BER, pierwsza ramka jest błędna, należy ją zignorować
% type: typ transmisji
% 0: sekwencja treningowa = identyfikacja kanału
% 1: estymacja SNR
% 2: transmisjia bitów
% mask: liczba bitów w podkanałach
% X(N/2,1): tansmitowane dane (binarnie)
% Y(N/2,1): odebrane dane (binarnie)
persistent lastX;
if isempty( lastX )
lastX = zeros( size( X ) );
disp( 'initialize last frame delay' );
end
if( trType ~=2 )
ber = zeros( size(X) );
return;
else
%error( 'BER: missing code' );
if( i>2 )
ber = zeros( size(X) );
else
ber = berp + ( (lastX-Y) ~= 0 );
end
end
lastX = X;
%===================== cmp ====================
function err = cmp( x1, y1, x2, y2 )
% porównywanie, wykresy, etc...
persistent lastX2;
if isempty( lastX2 )
lastX2 = zeros( size( x2 ) );
disp( 'initialize last frame delay' );
err = zeros( size(x2) );
return;
end
figure(1);
t = 1:length(y2);
plot( t, real( lastX2(t)), 'b.', t, real(y2(t)), 'bo', t, imag( lastX2(t)), 'rx', t, imag(y2(t)), 'ro' );
err = zeros( size(x1) );
lastX2 = x2;
Przykład 23.3
Przykład 23.3 (plik przyklad_23_3.m) prezentuje implementację w języku MATLAB sposobu użycia symulatora.% przyklad 23.3 clear all; close all; N = 512; % długość ramki Le = 16; % długość filtru TEQ Lp = 32; % długość prefiksu % etap 1, identyfikacja kanału: trType = 0; % estymacja kanału mask = 2*ones( 1, N/2 ); % maska (estymowane wszystkie podkanały) teqIR = [ 1, 0, 0, 0, 0, 0, 0, 0 ]; % ,,pusty'' TEQ - delta Kroneckera D = 0; % opóźnienie wprowadzane przez kanał FEQ = ones( 1, N/2 ); % ,,pusty'' FEQ h = adsl( 500, trType, mask, teqIR, D, FEQ, 'AWGN' ); % estymacja IR kanału % etap 2, obliczanie SNR clear adsl % kasowanie zmiennych persistant trType = 1; % estymacja SNR [ teqIR, D ] = AdslTEQ( h, Le, Lp, 0 ); % wyznaczanie korektora TEQ sir = conv( teqIR, h ); % obliczanie skróconej odpowiedzi impulsowej FEQ = fft( sir(1:N) ); % wyznaczanie korekora FEQ FEQ = 1./ ( FEQ .* exp( j*2*pi/N * (0:N-1)*(D) ) ); FEQ = FEQ( 1:N/2 ); snr = adsl( 512, trType, mask, teqIR, D, FEQ, 'AWGN' ); % estymacja wartości SNR % wizualizacja wyników: SNR, przydział bitów figure( 1 ); plot( 10*log10( snr ), 'k' ); axis( [0 255, 0, 50] ); title( 'SNR' ); xlabel( 'nr podkanału' ); ylabel( 'snr [dB]' ); bit = bitLoad( snr, 0 ); % przydział bitów - scenariusz ,,maksimum bitrate'' figure( 2 ), plot( 1:256, bit, 'k'); axis( [0, 255, 0, 15] ); title( 'przydział bitów' ); xlabel( 'numer podkanału' ); ylabel( 'ilość bitów' ); % etap 3, transmisja danych clear adsl trType = 2; % estymacja BER [ teqIR, D ] = AdslTEQ( h, Le, Lp, 0 ); sir = conv( teqIR, h ); FEQ = fft( sir(1:N) ); FEQ = 1./ ( FEQ .* exp( j*2*pi/N * (0:N-1)*(D) ) ); FEQ = FEQ( 1:N/2 ); mask = bit; ber = adsl( 16, trType, mask, teqIR, D, FEQ, 'AWGN' ); figure( 3 ); plot( ber ); title( 'liczba błędów w podkanałach' ); xlabel( 'nr podkanału' ); ylabel( 'liczba błędnie przesłanych symboli' );
Przykład 23.4
Przykład 23.4 (plik przyklad_23_4.m) prezentuje sposób estymacji odpowiedzi impulsowej dla modemu ADSL. Funkcja HEstim(...) jest częścią składową pliku adsl.m.% Przykład 23.4 % funkcja obliczająca i uśredniająca iloraz widma sygnału wyjściowego i wejściowego function H = HEstim( i, trType, X, Y, Hp ) % (...) H = Hp + Y./X; % obliczanie estymaty odpowiedzi impulsowej z uśrednionych danych H( N/2+2:N ) = conj( fliplr( H(2:N/2) ) ); H = H / (iter-1); h = ifft( H ); h = real( h );
Przykład 23.5
Program z przykładu 23.5 (przyklad_23_5.m) oblicza wagi korektora TEQ maksymalnie skracającego odpowiedź impulsową kanału.% Przykład 23.5
function [ w ] = EigMinISI( c, Ld, Lw, D )
% c - odpowiedź impulsowa kanału
% Ld - długość prefiksu
% Lw - długość korektora TEQ
% D - opóźnienie skróconej odpowiedzi impulsowej
% w - obliczona odpowiedź impulsowa filtra TEQ
Lc = length( c );
Lz = Lc+Lw-1-Ld-D;
C = zeros( Lw, Lc+Lw-1 );
for i = 1:Lw
C( i, i:i+Lc-1 ) = c;
end
Wd = zeros( Lc+Lw-1, Lc+Lw-1 ); Wr = Wd;
b = [ zeros( 1, D ), ones( 1, Ld ), zeros( 1, Lz ) ];
Wd = diag ( b ); Wr = diag ( 1-b );
Xd = C*Wd*C';
G = chol ( Xd );
Xr = C*Wr*C';
T = inv( G )'*Xr*inv( G );
[ V, D ] = eig( T );
LambdaMin = min( diag( D ) );
LambdaMinIdx = find( diag( D ) == LambdaMin );
w = V( :, LambdaMinIdx )'*inv( G )';