Rappresentare la convoluzione continua con Matlab

[gotthard]

Ahahahah!

A parte gli scherzi ( ), magari ora aspettiamo il continuo delle "risposte serie", quelle sul tuo vero problema.

[simo85]

quelle sul tuo vero problema.

Ma io non so se il problema è mio che non voglio capire, o di Matlab.

Come scrivevo, la funzione di Heaviside è rappresentata correttamente:

01.png (4.57 KiB) Osservato 3068 volte

La rappresentazione della convoluzione da problemi per qualsiasi intervallo di tempo usato:

02.png (8.73 KiB) Osservato 3068 volte

Suppongo per via di come Matlab tratta la funzione di Heaviside..

Codice: Seleziona tutto

fs = 1000;
t = -5:1/fs:5;
x = exp(2*t).*heaviside(-1*t);
h = heaviside(t - 3);
y = 1/fs * conv2(x,h, 'same');
plot(t, h);ylim([-1,2])
plot(t, y);ylim([0,0.7])

Concludo che per me la soluzione (se cosi la possiamo chiamare) al mio problema è quella usata nel messaggio 12.

Se invece la soluzione vera e propria non costa tanto e qualcuno ha voglia di condividerla facendomi notare dove sbaglio con Matlab, con un semplice esempio di codice, lo ringrazio in anticipo.

[gotthard]

Ma io non so se il problema è mio che non voglio capire, o di Matlab.

Secondo me il problema è di Matlab, e occorre trovare un modo per trattare la convoluzione tra segnali a tempo continuo quando almeno uno dei due termini della convoluzione è costituito da un gradino, aggirando il problema di Matlab dovuto al fatto che lo rappresenta (inevitabilmente) con un vettore finito di valori.

Ieri sera ho provato a ragionarci, ma non ti so dare la soluzione, visto che alcune prove che volevo fare le ho viste da te fatte nei post precedenti senza arrivare alla soluzione del problema..

Bisognerebbe ascoltare, quando hanno un po' di tempo, gli esperti DirtyDeeds e dimaios!

[simo85]

Secondo me il problema è di Matlab, e occorre trovare un modo per trattare la convoluzione tra segnali a tempo continuo

Io penso sia un MIO problema con Matlab, che non uso la testa e non lo uso bene.

Per esempio, ho preso un esempio dal libro Signals & Systems Using MATLAB by Luis F. Chaparro, giusto per vedere se io sto facendo qualcosa di malsano e di strano.

L'esempio è quello a pagina 274, io ho preso solo le prime righe di codice ed ho solo cambiato il timestep:

Codice: Seleziona tutto

Ts = 1/1000;
t = 0:Ts:2;
x = 10 * cos(20*pi*t + pi*(rand(1, 1) - 0.5));
h = 20 * exp(-10.^t).*cos(40 *pi*t);
y = Ts * conv(x, h, 'same');
plot(t, y);

Come si può osservare, si usa sempre la funzione conv... Allora il problema è nella funzione di heaviside.
Devo usare un array inizializzato a mano.

Ora ci provo.

[DirtyDeeds]

Ora ci provo.

No, fermati, e torniamo indietro.

Tu vuoi rappresentare il risultato della convoluzione in un certo intervallo di tempo, diciamo . Per poter calcolare la convoluzione in quell'intervallo di tempo hai bisogno di conoscere le due funzioni e in altri due intervalli di tempo, possibilmente infiniti. Poiché in Matlab non puoi memorizzare vettori infiniti, devi decidere dove fermarti: se memorizzi una porzione troppo stretta di e , il risultato della convoluzione non sarà corretto nell'intervallo e osserverai effetti di "bordo"; questo è il senso del mio messaggio [13].

Puoi anche ragionare in un altro modo: nel programma in [12] hai generato campioni di e tra -40 s e +40 s: in quale intervallo di tempo i risultati della convoluzione saranno validi?

[simo85]

No, fermati, e torniamo indietro.

Troppo tardi DirtyDeeds

Codice: Seleziona tutto

fs = 1000;
t = -5:1/fs:5;
hh = [t];
h = [t];

for i = 1:length(t)
   
   if t(i) < 0
      hh(i) = 1.0;
   else
      hh(i) = 0.0;
   end

   if t(i) > 3
      h(i) = 1.0;
   else
      h(i) = 0.0;
   end
end

x = exp(2*t).*hh;
x1 = exp(2*t).*heaviside(-1*t);

y = 1/fs * conv(x,h, 'same');
y1 = 1/fs * conv(x1,heaviside(t-3), 'same');

subplot(6,1,1);
plot(t,hh);ylim([-1,2]);grid on;xlabel('hh');
subplot(6,1,2);
plot(t,x);ylim([-1,2]);grid on;xlabel('exp(2t) * hh');
subplot(6,1,3);
plot(t,h);ylim([-1,2]);grid on;xlabel('h');
subplot(6,1,4);
plot(t,y);ylim([-1,2]);grid on;xlabel('1/fs * conv(x,h, same)');
subplot(6,1,5);
plot(t,x1);ylim([-1,2]);grid on;xlabel('exp(2t) * heaviside(- t)');
subplot(6,1,6);
plot(t,y1);ylim([-1,2]);grid on;xlabel('1/fs * conv(x1, heaviside(t-3), same)');

Avevo già fatto tutto in questi minuti.
Lo posto anche perché (a parte essere utile per la sintassi) poi lo perdo, so già come vanno queste cose.

Poiché in Matlab non puoi memorizzare vettori infiniti, devi decidere dove fermarti: se memorizzi una porzione troppo stretta di e , il risultato della convoluzione non sarà corretto nell'intervallo e osserverai effetti di "bordo"; questo è il senso del mio messaggio [13].

OK. Questo mi è chiaro.

Puoi anche ragionare in un altro modo: nel programma in [12] ... in quale intervallo di tempo i risultati della convoluzione saranno validi?

Penso per un intervallo , se non sbaglio.

[simo85]

Usando il codice del messaggio 12 per timestep differenti e con campioni da -3s a 3s, ottengo valori diversi sull'asse y.



È evidente che ho qualche concetto sballato nella mia testolina di rapa.

[DirtyDeeds]

Usando il codice del messaggio 12 per timestep differenti e con campioni da -3s a 3s, ottengo valori diversi sull'asse y.

Ricordati che stai campionando un segnale, c'è anche il problema dell'aliasing.

[simo85]

Mmm. Facciamo una cosa. Ci tornerò sopra più avanti e procedo per gradi, altrimenti finisce che non vado avanti con con il ripasso, e nemmeno con quello di Fisica.

Se continuo cosi rimango indietro, oltre a far si che il mio corpo cominci a formar parte della sedia e della scrivania.
Meglio andare a fare un po' di jogging sul lungo mare, magari mi si illumina la lampadina.

Poi sono anche riuscito a piantare il computer con Matlab e dover fare il reboot con REISUB.

[simo85]

Dunque, chissà che una passeggiata mi abbia schiarito le idee ?

Correndo ed osservando il mare, mi sono venute in mente varie cose..

Una è che stessi facendo confusione tra dominio continuo e discreto e che Matlab avesse qualche problema del rappresentare il risultato di una convoluzione come funzione a pezzi.
Era una ipotesi.. L'ho tenuta in tasca e quando sono arrivato a casa ho scritto questo script di prova:

Codice: Seleziona tutto

fs = 500;
t = 0:1/fs:2;
x = square(pi*t);
h = 3*square(4*pi*t);
y = 1/fs * conv(x,h,'same');
plot(t,y);


Ed il risultato è questo:

xxx.png (11.98 KiB) Osservato 3020 volte

È evidente che Matlab non ha questi problemi. I problemi ce li ho io...
Cosi, mi sono soffermato sul post 25 del grande DirtyDeeds (sottolineo in grassetto):

devi decidere dove fermarti: se memorizzi una porzione troppo stretta di e , il risultato della convoluzione non sarà corretto nell'intervallo e osserverai effetti di "bordo"; questo è il senso del mio messaggio [13].

Poi mi sono anche ricordato delle simulazioni sulla convoluzione discreta che ho fatto l'altro giorno a partire da questo thread: viewtopic.php?f=7&t=59863#p596740 e che in funzione della dimensione degli array (quindi del segnale, e nel caso continuo il concetto è lo stesso), cambia anche il risultato di uscita della convoluzione.

Di conseguenza, stavo trascurando completamente di simulare il tempo di (di lunghezza infinita) per un lunghezza finita abbastanza lunga da poter rappresentare il risultato della convoluzione con la lunghezza di :

Cosi ho scritto questo codice per applicare questo concetto:

Codice: Seleziona tutto

fs = 100;
t = -3:1/fs:6;
th = -12:1/fs:12;
x = exp(2*t).*heaviside(-1*t);
h = heaviside(th - 3) - heaviside(th - 12);
y = 1/(fs) * conv2(x, h, 'same');
plot(t,x);
plot(th,h);
plot(t,y);grid on; ylim([-0.5,1]);

Ed ecco il risultato:

yyy.png (4.82 KiB) Osservato 3020 volte

Ma anche quello che segue va bene, nel codice precedente si ha:

Codice: Seleziona tutto

fs = 100;
t = -3:1/fs:6;
th = -12:1/fs:12;
x = exp(2*t).*heaviside(-1*t);
h = heaviside(th - 3);
y = 1/(fs) * conv2(x, h, 'same');
plot(t,x);
plot(th,h);
plot(t,y);grid on; ylim([-0.5,1]);

Per il momento mi sento in dovere di ringraziare DirtyDeeds e dimaios, ed anche gotthard per i preziosi suggerimenti e correzioni.