Tesina del corso di Comunincazioni Numeriche

Sistema di trasmissione da simulare:

La frequenza di campionamento nella simulazione fs=1/Delt va scelta in modo tale che:

Procedere come segue:

Facoltativo:

spostare l'equalizzatore prima del campionatore e misurare il diagramma ad occhio dopo l'equalizzatore.

Nota 1: il ritardo introdotto da ciascun elemento dello shift register all'interno dell'equalizzatore e' sempre pari a Tb; se l'equalizzatore e' posto dopo il campionatore (che agisce negli istanti t_0+kTb) questo ritado corrsiponde ad un campione (di durata Tb), ma se l'equalizzatore e' posto prima del campionatore (che agisce negli istanti t_0+kTb) questo ritardo corrisponde a Ns campioni (ciascuno di durata Delt).

Nota 2: Questa misura serve per vedere se l'equalizzatore "apre" l'occhio in corrispondenza dell'istante di campionamento. Nel sistema reale l'equalizzatore andra' comunque posizionato dopo il campionatore (che agisce negli istanti t_0+kTb) e questa misura non sara' possibilile. In fase di progetto la misura puo' essere utile per capire se si sta operando correttamente o no.

Indice

close all
clear all
clc

Configurazioni

1. Parametri della simulazione

% Numero di bit da trasmettere (ignorato se si usa sequenza di m)
Nbit = 200;

% Frequenza di campionamento
Fs = 80;

% Durata delta di Kronecker
Nx = 100;

% Signal Noise Ratio
SNR_db = 10;
SNR_db_v = 0:5:25;

2. Parametri del canale di trasmissione e Bit rate (gruppo 5)

B0 = 1; % zero
z0 = -2*pi*B0;
B1 = 2; phi1 = pi/3; % due poli complessi coniugati
p1 = -2*pi*B1*(cos(phi1)+j*sin(phi1)); p2= conj(p1);
B2 = 3; phi2 = pi/6; % due poli complessi coniugati
p3 = -2*pi*B2*(cos(phi2)+j*sin(phi2)); p4 = conj(p3);

% Bit rate (bit/s)
Rbit = 4;

3. Flags

% Funzione di trasferimento del canale
fdB = 0; % show fdt in dB
fHc_a = 1; % show fdt "analogica" del canale
fHc_bt = 1; % show fdt del canale sintetizzata con trasformazione bilineare
fhc_iir = 1; % show risposta all'impulso della cascata di filtri iir

% Segnali del sistema di trasmissione e ricezione
fxm = 1; % show segnale in uscita dal mapper (DEBUG_SMCHTX = 1)
fx = 1; % show segnale trasmesso
fr = 1; % show segnale ricevuto (dopo il canale, prima del filtro hrx)
fw = 1; % show segnale in uscita dal filtro hrx

fhtot = 1; % show risposta all'impulso complessiva del sistema

% Diagramma ad occhio senza equalizzazione
fed = 1; % show eye diagram
fdp = 1; % show andamento distorsione di picco

fhtotsc =1; % show risposta all'impulso complessiva sottocampionata

% Equalizzatore zero-forcing
fdp_zf_Ne = 1; % show andamento dp(N_taps)
fPe_zf_Ne = 1; % show andamento Pe_zf(N_taps)
fg_zf_Ne = 0; % show g(t) in funzione del numero di prese
fg_zf_Ne_verbose = 1; % " in modalita' verbose
fed_zf = 1; % show eye diagram eq zero-forcing a monte del campionatore

% Equalizzatore di Wiener
fdp_mse_Ne =10; % show andamento dp(N_taps)
fPe_mse_Ne = 1; % show andamento Pe_mse(N_taps)
fg_mse_Ne = 0; % show g(t) in funzione del numero di prese
fg_mse_Ne_verbose = 1; % " in modalita' verbose
fPe_mse_SNR_Ne = 1; % show andamento Pe_mse(N_taps, SNR)
fed_mse = 1; % show eye diagram eq Wiener a monte del campionatore

4. Impostazioni di Debug

DEBUG_IIR = 1; % 1: 0 IIR poli semplici;  1: 3 IIR poly complex conj
DEBUG_GRUPPO = 5; % 5 e' il nostro gruppo! (1-14)
DEBUG_SMCHTX = 1; % 1: Sorgente di bit + Mapper + Convoluzione Htx
                  % 0: Sorgente di bit + ( Mapper & Htx in ROM )
DEBUG_LFSR = 1; % sorgente di bit costruita sulla base di sequenza di m
                % disponibile solo con DEBUG_SMCHTX = 1
DEBUG_LFSR_verbose = 0; % show le caratteristiche della sorgente di bit

% Debug altri gruppi
for G=1:1
switch DEBUG_GRUPPO
    case 1
        phi1=pi/8;phi2=pi/4;B0=3;B1=2;B2=1;Rbit=2.;
    case 2
        phi1=pi/8;phi2=pi/4;B0=1;B1=2;B2=3;Rbit=4.;
    case 3
        phi1=pi/8;phi2=pi/4;B0=3;B1=2;B2=1;Rbit=2.;
    case 4
        phi1=pi/3;phi2=pi/6;B0=1;B1=1;B2=1;Rbit=2.;
    case 5
        phi1=pi/3;phi2=pi/6;B0=1;B1=2;B2=3;Rbit=4.;
    case 6
        phi1=pi/3;phi2=pi/6;B0=1;B1=2;B2=1;Rbit=4.;
    case 7
        phi1=pi/6;phi2=pi/5;B0=0.5;B1=1.5;B2=1.5;Rbit=3.;
        Fs = 90;
    case 8
        phi1=pi/6;phi2=pi/5;B0=1;B1=1.5;B2=3;Rbit=3.;
        Fs = 90;
    case 9
        phi1=pi/4;phi2=pi/5;B0=3;B1=1.7;B2=2.5;Rbit=3;
        Fs = 90;
    case 10
        phi1=pi/4;phi2=pi/5;B0=2;B1=2;B2=1.5;Rbit=3;
        Fs = 90;
    case 11
        phi1=pi/8;phi2=pi/6;B0=1;B1=3;B2=2;Rbit=4;
    case 12
        phi1=pi/8;phi2=pi/6;B0=1;B1=2;B2=4;Rbit=4;
    case 13
        phi1=3*pi/10;phi2=2*pi/9;B0=1;B1=3;B2=2;Rbit=4;
    case 14
        phi1=3*pi/10;phi2=2*pi/9;B0=1;B1=2;B2=4;Rbit=4;
end
z0=-2*pi*B0; %zero
p1=-2*pi*B1*(cos(phi1)+j*sin(phi1));p2=conj(p1); % 2 poli complex coniugati
p3=-2*pi*B2*(cos(phi2)+j*sin(phi2));p4=conj(p3); % 2 poli complex coniugati
end

Init

Tbit = 1 / Rbit; % Durata del simbolo
Delt = 1 / Fs;   % Intervallo di campionamento

% Handle figures
fig_Hc_a = 1; % funzione di trasferimento del canale
fig_Hc_Vs_Hc_a = 2; % confronto H del canale con quella della trasf. bilin.
fig_hc_Vs_hc_a = 3; % confronto H del canale quella sintetizzata con l'IIR
fig_xm = 4; % segnale in uscita dal mapper
fig_x = 5; % segnale trasmesso
fig_r = 6; % segnale ricevuto
fig_w = 7; % segnale in uscita dal ricevitore
fig_htot = 8; % funzione di trasferimento complessiva del sistema
fig_e = 9; % diagramma ad occhio
fig_dp = 10; % andamento distorsione di picco
fig_htotsc = 11; % fdt complessiva del sistema sottocampionanta
fig_dp_zf_Ne = 21; % andamento dp(N_taps_zero-forcing)
fig_Pe_zf_Ne = 22; % andamento Pe_zf(N_taps_zero-forcing)
fig_g_zf_Ne = 23; % risposta all'impulso g dello zero-forcing
fig_ed_zf = 24; % eye diagram dopo zero-forcing a monte del campionatore
fig_dp_mse_Ne = 31; % andamento dp(N_taps_Wiener)
fig_Pe_mse_Ne = 32; % andamento Pe_mse(N_taps_Wiener)
fig_g_mse_Ne = 33; % risposta all'impulso g del Wiener
fig_Pe_mse_SNR_Ne = 100; % andamento Pe_mse(N_taps_Wiener, SNR)
fig_ed_mse = 34; % eye diagram dopo Wiener a monte del campionatore

disp(['Fs = ' num2str(Fs) ' Hz']);
disp(['Delt = ' num2str(Delt) ' s']);

SNR = 10.^(SNR_db./10);
SNR_v = 10.^(SNR_db_v./10);

% Controlli
if (Tbit*Fs - floor(Tbit*Fs)) ~= 0
    error('Tbit x Fs non e'' un numero intero di campioni!');
end

Fs = 80 Hz
Delt = 0.0125 s

Funzione di trasferimento del canale "analogica" Hc_a(f)

f = linspace(0, Fs/2, 1e4); % frequenza max osservabile = 1 / 2T
s = j*2*pi*f;

Hc_a = (s - z0) ./ ( (s - p1) .* (s - p2) .* (s - p3) .* (s - p4) );

% normalizzazione H(0) = 1
beta = Hc_a(1);
Hc_a = Hc_a ./ beta;

if fHc_a
    figure(fig_Hc_a)
    if fdB
        % dB
        semilogx(f, 20*log10(abs(Hc_a)))
        ylabel('|H_{c\_a}(f)|^2 [dB]');
        ylim([-40 10]);
    else
        % lineare
        semilogx(f, (abs(Hc_a)).^2)
        ylabel('|H_{c\_a}(f)|^2');
    end
    grid on
    title('Funzione di trasferimento del canale "analogica" H_{c\_a}(f)');
    xlabel('f [Hz]');
end;

clear s

Sintesi fdt del canale tramite trasformazione bilineare

z = exp(i*2*pi*f/Fs);

% Trasformazione bilineare (frequenza max osservabile = 1 / 2T)
s = 2/Delt .* (1 - z.^-1) ./ (1 + z.^-1);
Hc = (s - z0) ./ ( (s - p1) .* (s - p2) .* (s - p3) .* (s - p4) );
Hc = Hc ./ beta;

if fHc_bt
    figure(fig_Hc_Vs_Hc_a)
    if fdB
        % dB
        semilogx(f, 20*log10(abs(Hc)), 'b')
        hold on
        semilogx(f, 20*log10(abs(Hc_a)), 'r')
        ylabel('|H(f)|^2 [dB]');
        ylim([-40 10]);
    else
        % lineare
        semilogx(f, (abs(Hc)).^2, 'b')
        hold on
        semilogx(f, (abs(Hc_a)).^2, 'r')
        ylabel('|H(f)|^2');
    end
    grid on
    title('Confronto H_{c\_a}(f) con quella sintetizzata H_c(f)')
    xlabel('f [Hz]')
    legend('|H_c(f)|^2', '|H_{c\_a}(f)|^2', 'Location', 'BestOutside')
end

clear z s

Implementazione filtro IIR (cascata di IIR)

t = 0 : Delt : (Nx -1)*Delt;
x = [Fs zeros(1, Nx-1)]; % delta di Kronacker normalizzata a energia unitaria


% 1° IIR: zero
% y[k] = a1 x[k] + a2 x[k-1] - a3 y[k-1]
a1 = (2 - z0 * Delt) / Delt;
a2 = -(2 + z0 * Delt) / Delt;
a3 = -1;
y = zeros(size(x));
for k=1:length(x)
    if k > 1
        y(k) = a1 * x(k) + a2 * x(k-1) + a3 * y(k-1);
    else
        y(1) = a1 * x(1);
    end
end

if DEBUG_IIR

Coppia di poli complessi e coniugati

    % 2° IIR: prima coppia di poli complessi e coniugati
    % y[k] = a1 x[k] + a2 x[k-1] + a3 x[k-2] + a4 y[k-1] + a5 y[k-2]
    denominatore = (abs(p1) * Delt)^2 -4 * Delt * real(p1) +4;
    a1 = Delt^2 / denominatore;
    a2 = 2 * Delt^2 / denominatore;
    a3 = Delt^2 / denominatore;
    a4 = -2 * ((abs(p1) * Delt)^2 -4) / denominatore;
    a5 = -((abs(p1) * Delt)^2 +4 * Delt * real(p1) +4) / denominatore;

    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k>2
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * x(k-2) + a4 * y(k-1) +...
                   a5 * y(k-2);
        else
            if k > 1
                y(k) = a1 * x(k) + a2 * x(k-1) + a4 * y(k-1);
            else
                y(1) = a1 * x(1);
            end
        end
    end

    % 3° IIR: seconda coppia di poli complessi e coniugati
    % y[k] = a1 x[k] + a2 x[k-1] + a3 x[k-2] + a4 y[k-1] + a5 y[k-2]
    denominatore = (abs(p3) * Delt)^2 -4 * Delt * real(p3) +4;
    a1 = Delt^2 / denominatore;
    a2 = 2 * Delt^2 / denominatore;
    a3 = Delt^2 / denominatore;
    a4 = -2 * ((abs(p3) * Delt)^2 -4) / denominatore;
    a5 = -((abs(p3) * Delt)^2 +4 * Delt * real(p3) +4) / denominatore;

    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k>2
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * x(k-2) + a4 * y(k-1) +...
                   a5 * y(k-2);
        else
            if k > 1
                y(k) = a1 * x(k) + a2 * x(k-1) + a4 * y(k-1);
            else
                y(1) = a1 * x(1);
            end
        end
    end

else

Poli semplici

    % y[k] = a1 x[k] + a2 x[k-1] + a3 x[k-3]

    % coefficienti primo polo
    a1 = Delt / (2 - p1 * Delt);
    a2 = Delt / (2 - p1 * Delt);
    a3 = (2 + p1 * Delt) / (2 - p1 * Delt);
    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k>1
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * y(k-1);
        else
            y(1) = a1 * x(1);
        end
    end

    % coefficienti secondo polo
    a1 = Delt / (2 - p2 * Delt);
    a2 = Delt / (2 - p2 * Delt);
    a3 = (2 + p2 * Delt) / (2 - p2 * Delt);
    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k>1
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * y(k-1);
        else
            y(1) = a1 * x(1);
        end
    end

    % coefficienti terzo polo
    a1 = Delt / (2 - p3 * Delt);
    a2 = Delt / (2 - p3 * Delt);
    a3 = (2 + p3 * Delt) / (2 - p3 * Delt);
    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k==1
            y(1) = a1 * x(1);
        else
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * y(k-1);
        end
    end

    % coefficienti quarto polo
    a1 = Delt / (2 - p4 * Delt);
    a2 = Delt / (2 - p4 * Delt);
    a3 = (2 + p4 * Delt) / (2 - p4 * Delt);
    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k==1
            y(1) = a1 * x(1);
        else
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * y(k-1);
        end
    end

end

% Antitrasformata Hc_a(f) --> risposta all'impulso del canale hc_a(t)
A = (p4 - z0) / (p4 - p1) / (p4 - p2) / (p4 - p3);
B = (p3 - z0) / (p3 - p1) / (p3 - p2) / (p3 - p4);
C = (p2 - z0) / (p2 - p1) / (p2 - p3) / (p2 - p4);
D = (p1 - z0) / (p1 - p2) / (p1 - p3) / (p1 - p4);

hc_a = real( A*exp(p4*t) + B*exp(p3*t) + C*exp(p2*t) + D*exp(p1*t) );

% Normalizzazione a energia unitaria
sig2 = hc_a * hc_a' * Delt;
hc_a = hc_a ./ sqrt(sig2);
hc = real(y) ./ sqrt(sig2);

% plot grafici
if fhc_iir
    figure(fig_hc_Vs_hc_a)
    plot(t, hc_a)
    hold all
    plot(t, hc)
    title(['Confronto tra la risposta all''impulso teorica del filtro ' ...
           'analogico e quella all''uscita dell''IIR sintetizzato'])
    legend('h_{c\_a}(t)', 'h_c[z]', 'Location', 'BestOutside')
    xlabel('t [s]')
    grid on
end

clear a1 a2 a3 a4 a5 denominatore x y A B C D k

Trasmettitore

% Filtro di trasmissione: porta di durata Tbit a energia unitaria
htx = 1 / sqrt(Tbit) * ones(1, Tbit * Fs);

if DEBUG_SMCHTX

Sorgente di bit + Mapper + Convoluzione htx a blocchetti elementari

Sorgente di bit

    if DEBUG_LFSR
        % Sorgente di bit basata su sequenza di m (periodo 255)
        Nbit = 255;
        N_FF = 8;
        warning(['Il numero di bit da trasmettere e'' stato reimpostato',...
                 ' a ' num2str(Nbit)]);
        FF = [round(rand(1,N_FF-1)) 1]; % randomize
        b_t = zeros(1, Nbit); % vettore contenten i bit da trasmettere
        for K = 1:Nbit
            b_t(K) = FF(N_FF);
            % shift
            FF =[ xor(FF(4), xor(FF(5), xor(FF(6),FF(8)))) FF(1:7)];
        end

        if DEBUG_LFSR_verbose
            % Caratterizzazione della sorgente di bit
            disp('=============================================');
            disp('Caratterizzazione sorgente di bit');
            P1 = sum(b_t); P0 = Nbit - P1;
            disp(['P(1) / P(0) : ' num2str(P1) ' / ' num2str(P0) ' = ',...
                 num2str(P1/P0)]);
            % test sequenze 00 01 10 11
            i00 = 0; i01 = 0; i10 = 0; i11 = 0;
            for i=1:length(b_t)-1
                if b_t(1*i:1*i+1) == [1 0], i10=i10+1; end
                if b_t(1*i:1*i+1) == [0 1], i01=i01+1; end
                if b_t(1*i:1*i+1) == [1 1], i11=i11+1; end
                if b_t(1*i:1*i+1) == [0 0], i00=i00+1; end
            end
            disp('Test sequenze:');
            disp(['00 : ' num2str(i00)]); disp(['01 : ' num2str(i01)]);
            disp(['10 : ' num2str(i10)]); disp(['11 : ' num2str(i11)]);
            disp('=============================================');
        end
        clear N_FF FF i i00 i01 i10 i11 P1 P0 K
    else
        % Sorgente di bit statisticamente indipendenti {0 1}
        b_t = round( rand(1, Nbit) );
    end

Warning: Il numero di bit da trasmettere e' stato reimpostato a 255

Mapper 2 PAM

    s = 2*b_t -ones(size(b_t)); % {-1 1}
    m = [Fs zeros(1, Tbit * Fs -1)];
    x_temp = m' * s;
    xm_t = reshape(x_temp, [], Nbit * Tbit * Fs);

    % Plot segnale in uscita dal mapper
    t = 0 : Delt : (Nbit * Tbit - Delt);
    if fxm
        figure(fig_xm)
        stem(t, xm_t)
        title('Segnale xm(kT) in uscita dal mapper')
        xlabel('t [s]')
        grid on
    end

Filtro di trasmissione

    % Convoluzione con filtro di trasmissione htx
    x = conv(xm_t, htx) * Delt;
    x_t = x(1 : Nbit * Tbit * Fs); % tolgo Tbit * Fs -1  zeri in fondo

else

Versione ottimizzata

    % WARNING: il suo sviluppo e' morto presto: TODO check!

    % Sorgente di simboli + Mapper & convoluzione con htx

    % Sorgente di simboli 2 PAM statisticamente indipendenti +1, -1
    a = fix( rand(1, Nbit) * 2) * 2 - 1; % *Fs ma va via poi con conv

    % Mapper 2 PAM con in ROM direttamente la convoluzione con la
    % porta di durata Tbit normalizzata a energia unitaria (htx)
    x_temp = htx' * a;
    x_t = reshape(x_temp, [], Nbit * Tbit * Fs);

end

% Plot segnale trasmesso
if fx
    figure(fig_x)
    plot(t, x_t)
    title('Segnale trasmesso x(t)')
    xlabel('t [s]')
    grid on
    ylim([min(x_t) max(x_t)]*1.1) % stilistica!
end

clear m x_temp a s x

Canale di trasmissione

% 1° IIR: zero
% y[k] = a1 x[k] + a2 x[k-1] - a3 y[k-1]
a1 = (2 - z0 * Delt) / Delt;
a2 = -(2 + z0 * Delt) / Delt;
a3 = -1;
x = x_t; % uscita trasmettitore come ingresso
for k=1:length(x)
    if k > 1
        y(k) = a1 * x(k) + a2 * x(k-1) + a3 * y(k-1);
    else
        y(1) = a1 * x(1);
    end
end

if DEBUG_IIR

Coppia di poli complessi e coniugati

    % 2° IIR: prima coppia di poli complessi e coniugati
    % y[k] = a1 x[k] + a2 x[k-1] + a3 x[k-2] + a4 y[k-1] + a5 y[k-2]
    denominatore = (abs(p1) * Delt)^2 -4 * Delt * real(p1) +4;
    a1 = Delt^2 / denominatore;
    a2 = 2 * Delt^2 / denominatore;
    a3 = Delt^2 / denominatore;
    a4 = -2 * ((abs(p1) * Delt)^2 -4) / denominatore;
    a5 = -((abs(p1) * Delt)^2 +4 * Delt * real(p1) +4) / denominatore;

    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k>2
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * x(k-2) + a4 * y(k-1) +...
                   a5 * y(k-2);
        else
            if k > 1
                y(k) = a1 * x(k) + a2 * x(k-1) + a4 * y(k-1);
            else
                y(1) = a1 * x(1);
            end
        end
    end

    % 3° IIR: seconda coppia di poli complessi e coniugati
    % y[k] = a1 x[k] + a2 x[k-1] + a3 x[k-2] + a4 y[k-1] + a5 y[k-2]
    denominatore = (abs(p3) * Delt)^2 -4 * Delt * real(p3) +4;
    a1 = Delt^2 / denominatore;
    a2 = 2 * Delt^2 / denominatore;
    a3 = Delt^2 / denominatore;
    a4 = -2 * ((abs(p3) * Delt)^2 -4) / denominatore;
    a5 = -((abs(p3) * Delt)^2 +4 * Delt * real(p3) +4) / denominatore;

    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k>2
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * x(k-2) + a4 * y(k-1) +...
                   a5 * y(k-2);
        else
            if k > 1
                y(k) = a1 * x(k) + a2 * x(k-1) + a4 * y(k-1);
            else
                y(1) = a1 * x(1);
            end
        end
    end

else

Poli semplici

    % coefficienti primo polo
    a1 = Delt / (2 - p1 * Delt);
    a2 = Delt / (2 - p1 * Delt);
    a3 = (2 + p1 * Delt) / (2 - p1 * Delt);
    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k>1
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * y(k-1);
        else
            y(1) = a1 * x(1);
        end
    end

    % coefficienti secondo polo
    a1 = Delt / (2 - p2 * Delt);
    a2 = Delt / (2 - p2 * Delt);
    a3 = (2 + p2 * Delt) / (2 - p2 * Delt);
    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k>1
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * y(k-1);
        else
            y(1) = a1 * x(1);
        end
    end

    % coefficienti terzo polo
    a1 = Delt / (2 - p3 * Delt);
    a2 = Delt / (2 - p3 * Delt);
    a3 = (2 + p3 * Delt) / (2 - p3 * Delt);
    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k==1
            y(1) = a1 * x(1);
        else
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * y(k-1);
        end
    end

    % coefficienti quarto polo
    a1 = Delt / (2 - p4 * Delt);
    a2 = Delt / (2 - p4 * Delt);
    a3 = (2 + p4 * Delt) / (2 - p4 * Delt);
    x = y; % uscita iir precedente come ingresso
    for k=1:length(x)
        if k==1
            y(1) = a1 * x(1);
        else
            y(k) = a1 * x(k) + a2 * x(k-1) + a3 * y(k-1);
        end
    end

end

% normalizzazione filtro IIR
r_t = real(y) ./ sqrt(sig2);

if fr
    figure(fig_r)
    plot(t, r_t)
    title('Segnale ricevuto r(t)')
    xlabel('t [s]')
    grid on
end

clear a1 a2 a3 a4 a5 denominatore x y A B C D k

Filtro di ricezione

% Convoluzione con il filtro di ricezione
% porta di durata Tbit a energia unitaria
hrx = 1/sqrt(Tbit) * ones(1, Tbit * Fs);
w_t = conv(r_t, hrx) * Delt;
w_t = w_t(1 : Nbit * Tbit * Fs); % tolgo Tbit * Fs -1  zeri in fondo

if fw
    figure(fig_w)
    plot(t, w_t)
    title('Segnale in uscita dal filtro di ricezione y(t)');
    xlabel('t [s]')
    grid on
end

Risposta all'impulso complessiva del sistema

h_tot = conv(htx, conv(hc, hrx) * Delt) * Delt;

% Determinazione istante ottimo di campionamento (teorico sulla h_tot)
[htot_to, t0_teo] = max( abs(h_tot) );
t0_teo_sec = (t0_teo -1) * Delt;
disp(['t0_teo = ' num2str( t0_teo_sec ) ' s']);

if fhtot
    figure(fig_htot)
    plot(h_tot)
    title('Risposta all''impulso complessiva del sistema h_{tot}(t)')
    xlabel('t [z]')
    grid on
end

t0_teo = 0.35 s

Diagramma ad occhio

% scarto prime e ultima righe rispettivamente init filtri, zeri finali
ed_in = w_t(4*Tbit*Fs+1:floor(length(w_t)/(2*Tbit*Fs))*(2*Tbit*Fs));
ed = reshape(ed_in, 2 * Tbit * Fs, []);

ted = 0 : Delt : 2 * Tbit - Delt;
if fed
    figure(fig_e)
    plot(ted, ed)
    title('Diagramma ad occhio')
    xlabel('t mod Tbit [s]')
    grid on
end

% Determinazione istante ottimo di campionamento (dal diagramma ad occhio)
% come punto di minimo della distorsione di picco
ed_temp = ed';
eM = max( abs(ed_temp) );
em = min( abs(ed_temp) );

% Plot eM em
if fed
    hold on
    plot(ted, eM, 'bd');
    plot(ted, em, 'rd');
end

% Andamento distorsione di picco
dp = (eM(1:Tbit*Fs)-em(1:Tbit*Fs))./(eM(1:Tbit*Fs)+em(1:Tbit*Fs));

if fdp
    figure(fig_dp)
    % Per una migliore visualizzazione grafica faccio la spline e
    % visualizzo in pallini rossi i campioni per cui viene calcolata la
    % distorisione di picco.
    tdp = ted(1:Tbit*Fs);
    tdp_c = linspace(min(tdp), max(tdp), 1e3);
    dp_c = spline(tdp, dp, tdp_c);
    plot(tdp, dp, 'or')
    hold on
    plot(tdp_c, dp_c, '-b')
    title('Andamento distorsione di picco')
    xlabel('t [s]');
    grid on
    clear tdp tdp_c dp_c
end

% Ricerco l'istante ottimo di campionamento nel primo occhio, da 0 : Tbit s
[dp_noeq t0] = min(dp);

t0 = t0 + Tbit*Fs;
h_t0 = h_tot(t0);
t0_sec = (t0 -1) * Delt;
disp(['t0_ed = ' num2str(t0_sec) ' s (campione numero: ' num2str(t0) ')']);

% Distorsione di picco
disp(['Distorsione di picco dp = ' num2str(dp_noeq)]);

clear eM em ed_in

t0_ed = 0.3 s (campione numero: 25)
Distorsione di picco dp = 0.19399

Probabilita' di errore senza equalizzatore

% Probabilita' di errore del sistema 2 PAM ideale
Pe_2PAM = 0.5 * erfc( sqrt(SNR) );

% Probabilita' di errore del sistema 2 PAM in esame
Pe_noeq = 0.5 * erfc( h_t0 * (1-dp_noeq) * sqrt(SNR) );

display(['Pe(Ne) = ' num2str(Pe_noeq,'%.2e') ' (2 PAM ideale = ',...
         num2str(Pe_2PAM,'%.2e') ')']);

Pe(Ne) = 1.41e-01 (2 PAM ideale = 3.87e-06)

Sottocampionamento della risposta complessiva del sistema

% successori
t_sc_succ = t0 - floor(t0/(Tbit*Fs))*(Tbit*Fs) : (Tbit*Fs) : t0 - 1;
h_tot_sc_succ = h_tot(t_sc_succ);
disp(['Numero di successori = ' num2str(length(t_sc_succ))]);
N_succ = length(t_sc_succ);

% precursori
t_sc_prec = t0 + (Tbit*Fs) : (Tbit*Fs) : length(h_tot);
h_tot_sc_prec = h_tot(t_sc_prec);
disp(['Numero di precursori = ' num2str(length(t_sc_prec))]);
N_prec = length(t_sc_prec);

% sottocampionamento funzione di trasferimento complessiva del sistema
t_sc = [t_sc_succ t0 t_sc_prec];
h_tot_sc = h_tot(t_sc);

if fhtotsc
    figure(fig_htotsc)
    %t = 0:Delt:Delt*(length(h_tot)-1);
    plot(h_tot)
    hold on
    plot(t_sc_succ, h_tot_sc_succ, 'ro')
    plot(t0, h_tot(t0), 'go')
    plot(t_sc_prec, h_tot_sc_prec, 'md')
    grid on
    title(['Risposta all''impulso complessiva del sistema h_{tot}(t) ',...
           'campionata'])
end

Numero di successori = 1
Numero di precursori = 5

Equalizzatore Zero-Forcing

Ne = 1;
dp_old = 1;
dp_new = dp_noeq;
dp_zf_Ne(1) = dp_noeq;
sigma_nu_zf_Ne(1) = NaN;

%while (dp_old - dp_new)/dp_old > 0.20
while dp_old - dp_new > 0.01

Ciclo per la determinazione del numero di prese

    % incremento del numero di prese del filtro FIR di equalizzazione
    Ne = Ne +1;
    dp_old = dp_new;

    % costruzione della matrice H (N_prec + N_suc + Ne) x Ne
    H = zeros(N_prec + N_succ + Ne, Ne);
    H(1:N_prec+N_succ+1,1) = h_tot_sc';
    for I = 2:Ne
        H(:,I) = circshift( H(:,I-1), 1);
    end

    Rm = eye(Ne);

Ciclo per la determinazione della combinazione di righe di H

    % in realta' ho trovato che va a prendere sempre dalla seconda riga

    for I = 1:N_succ+N_prec

        % estrazione sottomatrice quadrata
        H1 = H(I:I+Ne-1,:);

Selezione del miglior C_ott tra le colonne di inv H1 disponibili

        g = eye(Ne);

        C = H1 \ g; % e' quindi una matrice dei possibili vettori di taps

        g_real = H * C; % compaiono gli elementi che non sono stati forzati
                        % a zero

        % Distorsione di picco: definizione
        % sum con l!=0 { | h(t0 -l*T )| } / h(t0)
        dp_new_zf = ( sum( abs(g_real) ) -1 ) ./ max( abs(g_real) );

        % ricerca delle distorsione di picco minima data dai possibili taps
        [ dp_new_zf, ind] = min ( dp_new_zf );

        if dp_new > dp_new_zf
            % Abbiamo ottenuto una distorsione di picco minore!
            dp_new = dp_new_zf;
            c_ott_zf = C(:, ind);
            dp_zf_Ne(Ne) = dp_new;
            sigma_nu_zf_Ne(Ne) = c_ott_zf' * Rm * c_ott_zf;
            g_saver_zf(1:N_succ+N_prec+Ne, Ne) = g_real(:,ind);
        end

    end

end

dp_zf = dp_new;

display('==============================================================');
display('Equalizzatore zero forcing');
display(['Numero di prese :', num2str(Ne)]);
display(['Distorsione di picco ottenuta : ', num2str(dp_zf), ...
        ' (' num2str((dp_noeq-dp_zf)/dp_noeq*100) '%)' ]);
display(['C_ott : ', num2str(c_ott_zf')]);

if fg_zf_Ne
    % mostra la risposta all'impulso g(t)
    % risposta all'impulso del sistema di trasmissione equalizzata
    figure(fig_g_zf_Ne)
    if fg_zf_Ne_verbose
        for I=2:min(Ne,10)
            plot(g_saver_zf(:,I))
            legend(num2str(I))
            title('Risposta all''impulso g(t) con eq. zero-forcing')
            grid on
            pause
        end
    end
    plot(g_saver_zf(:,2:Ne))
    title('Risposta all''impulso g(t) con eq. zero-forcing')
    grid on
end

if fdp_zf_Ne
    % mostra l'andamento della distorsione di picco all'aumentare del
    % numero di taps dell'equalizzatore zero-forcing
    figure(fig_dp_zf_Ne)
    plot(1, dp_zf_Ne(1), 'rd')
    hold on
    plot(2:Ne, dp_zf_Ne(2:Ne), '-bd')
    title(['Andamento della dp al variare del numero di taps dell',...
           '''eq. zero-forcing'])
    xlabel('Numero di prese Ne')
    ylabel('dp_z_f')
    grid on
end

clear I dp_old dp_new dp_new_zf g C g_real ind H1

==============================================================
Equalizzatore zero forcing
Numero di prese :7
Distorsione di picco ottenuta : 0.0048124 (97.5192%)
C_ott : -0.042898
         3.3495
         0.12344
         0.34232
        -0.082805
         0.05366
        -0.018226

Probabilita' di errore con equalizzatore zero-forcing

% Pe_noeq = 0.5 * erfc( h_t0 * (1-dp_noeq) * sqrt(SNR) ); Senza eq.
% h_t0 qui e' forzato a 1 dall'equalizzatore

Pe_zf_Ne = 0.5 * erfc( (1 - dp_zf_Ne) .* sqrt(SNR./sigma_nu_zf_Ne) );

display(['Pe(Ne) : ' num2str(Pe_zf_Ne,'%.2e\t')]);

if fPe_zf_Ne
    figure(fig_Pe_zf_Ne)
    plot(Pe_zf_Ne)
    hold on
    plot(Pe_noeq*ones(length(Pe_zf_Ne)), 'r')
    title('Probabilita'' di errore per eq. zero-forcing')
    xlabel('Numero di prese Ne')
    ylabel('Pe_{zf}(Ne)')
    grid on
end

Pe(Ne) : NaN	1.28e-01	1.06e-01	1.01e-01	
               9.74e-02	9.45e-02	9.34e-02

Eye diagram dopo filtro eq zero-forcing (a monte del campionatore)

% Nota 1: il ritardo introdotto da ciascun elemento dello shift register
% all'interno dell'equalizzatore e' sempre pari a Tb; se l'equalizzatore e'
% posto dopo il campionatore (che agisce negli istanti t_0+kTb) questo
% ritado corrsiponde ad un campione (di durata Tb), ma se l'equalizzatore
% e' posto prima del campionatore (che agisce negli istanti t_0+kTb) questo
% ritardo corrisponde a Ns campioni (ciascuno di durata Delt).

taps = zeros( 1, Tbit * Fs * length(c_ott_zf) );
for k = 1:length(c_ott_zf)
    taps(1 + (k-1)*Tbit*Fs) = c_ott_zf(k);
end

ed_zf_in = conv(w_t, taps); % convoluzione numerica!

% scarto un po' di porkeria dal diagramma ad occhio
ed_zf_in = ed_zf_in(2*Tbit*Fs+1 : floor(length(ed_zf_in)/(2*Tbit*Fs))...
                                  *(2*Tbit*Fs) -ceil(Ne/2)*2*Tbit*Fs );
ed_zf = reshape(ed_zf_in, 2 * Tbit * Fs, []);

if fed_zf
    figure(fig_ed_zf)
    ted_zf = 0 : Delt : 2*Tbit - Delt;
    plot(ted_zf, ed_zf)
    title('Diagramma ad occhio dopo zero forcing a monte del campionatore')
    xlabel('t mod Tbit [s]')
    grid on
end

clear taps k ted_zf ed_zf_in

Equalizzatore di Wiener

Ne = 1;
dp_old = 1;
dp_new = dp_noeq;
dp_mse_Ne(1) = dp_noeq;
sigma_nu_mse_Ne(1) = NaN;
h_t0_mse_Ne(1) = NaN;

% i primi passi sono identici allo zero-forcing

while (dp_old - dp_new)/dp_old > 0.01

Ciclo per la determinazione del numero di prese

    % incremento del numero di prese del filtro FIR di equalizzazione
    Ne = Ne +1;
    dp_old = dp_new;

    % costruzione della matrice H (N_prec + N_suc + Ne) x Ne
    H = zeros(N_prec + N_succ + Ne, Ne);
    H(1:N_prec+N_succ+1,1) = h_tot_sc';
    for I = 2:Ne
        H(:,I) = circshift( H(:,I-1), 1);
    end

    Rm = eye(Ne);

    % Se hrx(t) e' tale che la sua funzione di autocorrelazione soddisfa il
    % I criterio di Nyquist ( Rhrx(kT) = 0 per ogni k != 0 ) allora
    % Rmu = N_0 / 2 * I
    % c_ott = ( H' * H + N_0 / 2*E * I )^-1 * H' u
    Rmu = 0.5 * SNR^-1 * eye(Ne);

    u = eye(Ne + N_succ + N_prec);

    C = ( H'*H + Rmu) \ H' * u;

    g_real = H * C;

    % Distorsione di picco: definizione
    % sum con l!=0 { | h(t0 -l*T | } / h(t0)
    dp_new_mse = (sum(abs(g_real)) -max(abs(g_real)))./max(abs(g_real));

    % ricerca delle distorsione di picco minima data dai possibili taps
    [ dp_new_mse, ind] = min ( dp_new_mse );

    % Abbiamo ottenuto una distorsione di picco minore!
    dp_new = dp_new_mse;
    c_ott_mse = C(:, ind);
    dp_mse_Ne(Ne) = dp_new;
    sigma_nu_mse_Ne(Ne) = c_ott_mse' * Rm * c_ott_mse;
    h_t0_mse_Ne(Ne) = max(abs(g_real(:,ind)));
    g_saver_mse(1:N_succ+N_prec+Ne, Ne) = g_real(:,ind);

end

dp_mse = dp_new;

display('==============================================================');
display('Equalizzatore Wiener');
display(['Numero di prese :', num2str(Ne)]);
display(['Distorsione di picco ottenuta : ', num2str(dp_mse), ...
        ' (' num2str((dp_noeq-dp_mse)/dp_noeq*100) '%)' ]);
display(['C_ott : ', num2str(c_ott_mse')]);

if fg_mse_Ne
    % mostra la risposta all'impulso g(t)
    % risposta all'impulso del sistema di trasmissione equalizzata
    figure(fig_g_mse_Ne)
    if fg_mse_Ne_verbose
        for I=2:min(Ne,10)
            plot(g_saver_mse(:,I))
            legend(num2str(I))
            title('Risposta all''impulso g(t) con eq. Wiener')
            grid on
            pause
        end
    end
    plot(g_saver_mse(:,2:Ne))
    title('Risposta all''impulso g(t) con eq. Wiener')
    grid on
end

if fdp_mse_Ne
    % mostra l'andamento della distorsione di picco all'aumentare del
    % numero di taps dell'equalizzatore Wiener
    figure(fig_dp_mse_Ne)
    plot(1, dp_mse_Ne(1), 'rd')
    hold on
    plot(2:Ne, dp_mse_Ne(2:Ne), '-bd')
    title(['Andamento della dp al variare del numero di taps dell',...
           '''eq. di Wiener'])
    xlabel('Numero di prese Ne')
    ylabel('dp_{MSE}')
    grid on
end

clear I dp_old dp_new dp_new_mse g C g_real ind

==============================================================
Equalizzatore Wiener
Numero di prese :9
Distorsione di picco ottenuta : 0.072459 (62.6472%)
C_ott : -0.047571
         2.1389
         0.059364
         0.1397
        -0.036436
         0.017454
        -0.0052259
         0.0023727
        -0.00085894

Probabilita' di errore con equalizzatore di Wiener

Pe_noeq = 0.5 * erfc( h_t0 * (1-dp_noeq) * sqrt(SNR) ); %Senza eq.

Pe_mse_Ne = 0.5 * erfc( h_t0_mse_Ne .* (1 - dp_mse_Ne) .* ...
                        sqrt(SNR./sigma_nu_mse_Ne) );

display(['Pe(Ne) = ' num2str(Pe_mse_Ne,'%.2e\t')]);

if fPe_mse_Ne
    figure(fig_Pe_mse_Ne)
    plot(Pe_mse_Ne)
    hold on
    plot(Pe_noeq*ones(length(Pe_mse_Ne)), 'r')
    title('Probabilita'' di errore per eq. di Wiener')
    xlabel('Numero di prese Ne')
    ylabel('Pe_{mse}(Ne)')
    grid on
end

Pe(Ne) = NaN	1.31e-01	1.16e-01	1.12e-01	1.10e-01
               1.09e-01	1.09e-01	1.08e-01	1.08e-01

Eye diagram dopo filtro eq MSE (a monte del campionatore)

% Nota 1: il ritardo introdotto da ciascun elemento dello shift register
% all'interno dell'equalizzatore e' sempre pari a Tb; se l'equalizzatore e'
% posto dopo il campionatore (che agisce negli istanti t_0+kTb) questo
% ritado corrsiponde ad un campione (di durata Tb), ma se l'equalizzatore
% e' posto prima del campionatore (che agisce negli istanti t_0+kTb) questo
% ritardo corrisponde a Ns campioni (ciascuno di durata Delt).

taps = zeros( 1, Tbit * Fs * length(c_ott_mse) );
for k = 1:length(c_ott_mse)
    taps(1 + (k-1)*int16(Tbit*Fs)) = c_ott_mse(k);
end

ed_mse_in = conv(w_t, taps); % convoluzione numerica!

% scarto prima e ultima riga rispettivamente init filtri, zeri finali
ed_mse_in = ed_mse_in(2*Tbit*Fs+1 : floor(length(ed_mse_in)/(2*Tbit*Fs))...
                                    *(2*Tbit*Fs) -ceil(Ne/2)*2*Tbit*Fs );
ed_mse = reshape(ed_mse_in, 2 * Tbit * Fs, []);

if fed_mse
    figure(fig_ed_mse)
    ted_mse = 0 : Delt : 2*Tbit - Delt;
    plot(ted_mse, ed_mse)
    title('Diagramma ad occhio dopo Wiener a monte del campionatore')
    xlabel('t mod Tbit [s]')
    grid on
end

clear taps k ted_mse ed_mse_in

Campionatore

Blocchetto inutile perche' non simulo il rumore!

z_t = w_t( t0 : (Tbit*Fs) : length(w_t) );

Decisore a soglia zero

Blocchetto inutile perche' non simulo il rumore!

b_r = z_t > 0;