
% Illustration et comparaison des egaliseurs ZF et MMSE (structure filtre
% transverse RIF)
%
% Marc Castella, 28/03/2006
% INT (Institut National des Telecommunications), Departement CITI
% Evry, France


clear; close all;

choix_canal = 3;        
M = 20;               %longueur de l'egaliseur
SNR = 10;            %rapport signal/bruit en dB

Naffich = 60;       %nombre de symboles affiches

% canaux de Proakis p.631
h1 = [0.04 -0.05 0.07 -0.21 -0.5 0.72 0.36 0 0.21 0.03 0.07].';
h2 = [0.407 0.815 0.407].';
h3 = [0.227 0.460 0.688 0.460 0.227].';
switch choix_canal
    case 1,
        c = h1;
    case 2,
        c = h2;
    case 3,
        c = h3;
end;

[C,W] = freqz(c,1,128,1);

L = length(c);                  %longueur du canal
Ns = 10000;                      %nombre echantillons
a = 2*(unidrnd(2,Ns,1)-1.5);    %symboles
x = filter(c,1,a);              %observation en sortie canal
figure(5); subplot(2,1,1); plot(x(1:Naffich)); title('Sortie canal non bruite');
x = x + 10^(-SNR/20)*randn(Ns,1);
figure(5); subplot(2,1,2); plot(x(1:Naffich)); title('Sortie canal bruite');

% calcul du filtre ZF
decal = floor((L+M-1)/2);
h = zeros(L+M-1,1); h(decal) = 1;    %RI globale souhaitee
e_zf = convmtx(c,M)\h;          %egaliseur zero-forcing
y = filter(e_zf,1,x);

disp(strcat('EQM du ZF : ',num2str(norm(y(decal:Ns)-a(1:Ns-decal+1)).^2/length(y))))

figure(1);
subplot(2,1,1); 
plot((1:Naffich),y(decal:decal+Naffich-1),'b',(1:Naffich),a(1:1+Naffich-1),'r-.');
title('Egaliseur ZF');
subplot(2,1,2);
stem(conv(e_zf,c));
figure(2);
[E_zf,W] = freqz(e_zf,1,128,1); semilogy(W,abs(C),'b',W,abs(E_zf),'r');
title('Reponse en frequence de l''egalizeur ZF');
legend('Canal','Egaliseur','Location','SouthWest');

% calcul du filtre MMSE
X = convmtx(x,M); X_short = X(M:end-M+1,:);

%estimateurs des correlations
Rxx = X_short'*X_short/(Ns-M+1);
Rxa = X(1+(decal-1):end-M,:)'*a(1:Ns-(decal-1)-1)/(Ns-decal); %attention decalage = decal - 1

e_mmse = Rxx\Rxa;       %egaliseur mmse (Wiener)
y = filter(e_mmse,1,x);

disp(strcat('EQM du MMSE : ',num2str(norm(y(decal:Ns)-a(1:Ns-decal+1)).^2/length(y))))

figure(3);
subplot(2,1,1);
plot((1:Naffich),y(decal:decal+Naffich-1),'b-',(1:Naffich),a(1:1+Naffich-1),'r-.');
title('Egaliseur MMSE');
subplot(2,1,2);
stem(conv(e_mmse,c));
figure(4);
[E_mmse,W] = freqz(e_mmse,1,128,1); semilogy(W,abs(C),'b',W,abs(E_mmse),'r');
title('Reponse en frequence de l''egaliseur MMSE (Wiener)');
legend('Canal','Egaliseur','Location','SouthWest');

% autres canaux possibles (pour illustration)
% ----------------------------------
% c = fir1(50,[0.6 0.62],'stop').'; %bof, filtre moyennenement convaincant
% % ------
% z1 = 0.99*exp(i*pi/2);
% num_h = conv([1 -z1],[1 -conj(z1)]); num_h = conv(num_h,num_h);num_h = conv(num_h,num_h);
% den_h = 1;
% c = impz(num_h,den_h);
