filtro numerico

Tanti auguri :-) È uno dei grandi problemi che l'umanità non è ancora riuscita a risolvere interamente.

Scherzi a parte, si chiama modellazione inversa o, nel caso specifico, deconvoluzione. Il modo più becero è passare nel dominio della frequenza (Fourier), dividere y per la trasformata del tuo filtro (in questo caso 1-a*exp(-j w)) e quindi tornare nel dominio del tempo. Ovviamente non otterrai esattamente il segnale iniziale, per via del rumore e(t) e delle condizioni iniziali su x (se non le conosci).

Ciao

Pasu

olvere

ividere

di

Mi genufletto.

Piccio.

"Pasu" ha scritto nel messaggio news:gc580r$p49$ snipped-for-privacy@aioe.org...

Caspita che sfiga :-))

Allora, vediamo se ho capito bene, ho il vettore Y = [y(1) y(2) ... y(n)], creo il vettore F = [ y(1)/(1-a*exp(-j w)) y(2)/(1-a*exp(-j w)) .... y(n)/(1-a*exp(-j w)) ] faccio la trasformata inversa di F e ottengo qualcosa di simile a x ...

corretto? hai un link utile che spieghi il metodo e la teoria che ci sta dietro?

Ciao e grazie

Non proprio Y=fft(y) dove y è il vettore delle uscite, nel tempo. H=fft(h) dove h è il filtro, nel tuo caso h=[1,a,0,0,0,0,0,...,0], riempito di zeri in modo da avere lunghezza pari a quella di y.

X=Y./H (divisione elemento per elemento)

xcap= ifft(X)

Link, a parte la solita wikipedia, non ne ho. Comunque qualsiasi libro di segnali o dsp lo spiega. È la proprietà più notevole della trasformata di Fourier: un prodotto di convoluzione (filtraggio nel tempo) si trasforma in un prodotto algebrico (moltiplicazione elemento per elemento).

Se hai Matlab puoi provare una cosa del genere:

%%%%%%%%%%%%%%%%%% N=256; x=sin(2*pi*(0:N-1)*10/N)'; %il tuo segnalino e=0.1*randn(N,1); %rumore h=[1; -0.5]; %filtro, a=0.5

y=filter(h,1,x)+e; Y=fft(y); H=fft(h,N); Xcap=Y./H; xcap=ifft(Xcap);

plot([x,xcap])

%%%%%%%%%%%%%%%%%%

Ovviamente più piccolo è il rumore, migliore sarà la stima di x. Se il rumore è nullo otterrai praticamente x, a parte un transitorio iniziale dovuto alle condizioni iniziali, tanto più lungo quanto a è vicina a 1 o -1.

Prego!

P.

Uh? Mavalà...

"Pasu" ha scritto nel messaggio news:gc7sig$ksa$ snipped-for-privacy@aioe.org...

Wow funziona, grazie ora però mi scontro con un problema logico a cui non riesco a dare una soluzione: ho un segnale y(t) che ha questa caratteristica, se chiamo z(t) = y(t) - a*y(t-1) allora so che z(t) = x(t) - a*x(t-1) + e(t) con e(t) rumore gaussiano bianco a media nulla, conosco y(t) ma ciò che voglio è x(t) che è un segnale incognito. Se questo algoritmo di deconvoluzione lo applico a z(t) ottengo y(t), sembra ovvio per come ho costruito z(t) ma sapendo che così costruito z(t) è legato anche a x(t) come detto in precedenza non ci sto capendo + nulla :-)) Hai qualche altra dritta? :-))) Ancora grazie ciao

Certo che funziona...L'ho provato prima di postarlo, che credi? :-)

z(t) = x(t) - a*x(t-1) + e(t)

Ma no, perchè?

x(t) deve passare due volte in quel filtro, quindi devi convolvere i due filtri (o moltiplicarli nel dominio della frequenza): h(t)**h(t) (** è il prodotto di convoluzione) Nel dominio della trasformata zeta: H(z)*H(z)= (1-az^-1)*(1-az^-1) = 1 -2az-1 + a^2 z^-2

quindi z(t)= x(t) - 2a x(t-1) +a^2 x(t-2) + e(t) - a e(t-1) (il rumore passa solo attraverso il secondo filtro)

(scusa il casino di simboli... Non è mai consigliabile usare z come nome di variabile!)

Pensaci bene, anche senza fare conti: se un segnale passa due volte attraverso un filtro, come fa ad essere ancora quello di partenza (o quello intermedio)? Il fatto che i filtri siano identici non importa, sono comunque due "scatolotti" distinti.

Ciao

P.

"Pasu" ha scritto nel messaggio news:gcb3ss$lne$ snipped-for-privacy@aioe.org...

CUTTONE

Caro Pasu, i casi sono due, o ti sei messo d'accordo con l'amico Jack per fare questa bellissima figura oppure hai veramente due maroni grossi come una casa, in ogni caso non ci ho capito nulla e ti faccio i miei piu' sinceri complimenti ;-))))

Da grande voglio essere come voi ^_^

Saluti da Drizzt.

Per ipotesi, il segnale y(t) è costruito in modo che y(t) - a*y(t-1) = x(t) - a*x(t-1) + e(t) i cui termini al secondo membro sono incogniti (eccetto a).

Perdonami, non ti seguo, capisco i passaggi matematici ma non colgo l'obiettivo (scusami sono ancora sbarbatello in materia :-))

Ciao e Ri-grazie :)

Eheh... Grazie Jack, ottimo lavoro! Ha funzionato! :-)

Beh, grazie, ma se non ci hai capito nulla vuol dire che i maroni non sono poi così grossi... Ti assicuro comunque che 'ste cose le trovi a pagina 4 di qualsiasi libro di segnali. La mia "spiegazione" non era molto matematica, era più orientata alla programmazione (Matlab)... Forse è per questo che non mi hai seguito.

Anch'io!

Ciao!

Pasu

Se conosci la trasformata Z è tutto molto più facile. In caso contrario, sostituisci le variabili nel tempo (y->x). È un po' più macchinoso, ma funziona lo stesso.

y(t)= x(t) - a x(t-1) y(t-1)=x(t-1) - a x(t-2) z(t)= y(t) - a y(t-1)= [x(t) - a x(t-1)] - a [x(t-1) - a x(t-2)] = = x(t) - 2a x(t-1) -a^2 x(t-2).

Prova a riscrivere tu questi passaggi.

Per semplificare, se il filtro fosse solo uno scalamento (un solo campione) la convoluzione sarebbe un normale prodotto, ma la sostituzione di variabili funziona allo stesso modo y= a x z= a y -> z=a (a x) = a^2 z.

Convinto ora?

P.

"Pasu" ha scritto nel messaggio news:gcbe3u$5jk$ snipped-for-privacy@aioe.org...

ecco questo non è vero, l'unica cosa che so del segnale y(t) è che se prendo la quantità y(t) - a*y(t-1) questa equivale a y(t) - a*y(t-1) = x(t) - a*x(t-1) + e(t) (dimentica il problema che avevo posto all'inizio del thread) allora ho pensato, toh che bello mi calcolo un nuovo segnale s(t) = y(t) - a*y(t-1) quindi sarà anche s(t) = x(t) - a*x(t-1) + e(t) quindi tu mi hai dato la soluzione di deconvoluzione, l'ho applicata ad s(t) ma mi torna y(t) e non x(t) come io volevo. Lo so c'è qualcosa che non quadra in maniera logica, ci sto ancora riflettendo.

nuovamente grazie

Ti giuro che ho letto più volte il tuo post ma non ci ho capito niente.

P.>> y(t)= x(t) - a x(t-1)

Non è vero perchè ho tralasciato il rumore e(t)? Volevo solo semplificare le cose, visto che togliendolo non cambia la sostanza del problema.

Perchè? Quello che hai scritto è una cosa parecchio diversa. Ma anche se fosse così, questo non implica y(t)=x(t), perchè non conosci le condizioni iniziali (il valore di y(t-1), x(t-1) e x(t) all'istante iniziale), che in generale sono diverse tra loro.

Quindi stai cambiando il problema?

E dai. Ma perchè?? A parte le condizioni iniziali che nessuno ti garantisce siano uguali per x ed y, qui ci hai anche aggiunto il rumore...

Riscrivila così: y(t)= x(t) - a x(t-1) + e(t) + a y(t-1)

In parole, per calcolare il valore di y(t) all'istante attuale, ti occorre sapere il valore di x all'istante attuale e a quello precedente, il valore di e all'istante attuale e il valore di y calcolato al passo precedente.

Quei x(t-1) e y(t-1), vogliono dire che ci occorre il valore di x ed y all'istante (al campione) precedente a quello considerato. Per averlo, occorre una memoria nella quale conservare il dato attuale (salvato una volta calcolato il valore di y(t) ) per l'istante successivo. Spero sia chiaro a questo punto che quella per x e quella per y sono due memorie distinte, che in generale contengono valori diversi e che vengono inizializzate indipendentemente.

Ti assicuro che per capire le cose a questo livello sono sufficienti strumenti matematici che si imparano in prima superiore. In più si tratta solo di avere chiaro il concetto di cella di memoria, moltiplicatore e sommatore. Son sicuro che sul libro o sulla dispensa su cui stai studiando c'è almeno un disegno dell'implementazione di un filtro FIR.

Ciao

Pasu

Ok

Perfetto, è proprio quello che vorrei avere, x e y diverse, nella equazione di cui sopra io conosco solo y(t) e y(t-1) (oltre ad a). Voglio ricavarmi x(t). (posso anche fare un'ipotesi sul valore iniziale delle x e y).

Sì i concetti base li ho + o - chiari, ma non riesco ad applicarli al mio caso.

grazie ancora per la pazienza

Prima di tutto l'ultima volta ti ho detto una mezza (ok, 3/4 di) vaccata. Il filtro di cui sopra ( y(t)= x(t) - a x(t-1) + e(t) + a y(t-1) ) lo si può implementare con una cella di memoria sola. Infatti basta sommare y(t) e x(t) (o in questo caso ay(t)-ax(t) ) PRIMA di memorizzarli.

Allora, appurato che ora il problema è deconvolvere

y(t)= x(t) + (-a) x(t-1) + e(t) - (-a) y(t-1)

o in forma più generale y(t)= x(t) + B x(t-1) + e(t) - A y(t-1)

(A e B sono le convenzioni che usa Matlab, le segno maiuscole per distinguerle dalla tua a).

Tu, come dici, conosci y(t) e h(t) (il filtro) e, al più le condizioni iniziali della famosa memoria di cui sopra.

Per i valori di A e B che hai scelto, hai che il filtro h, in z è:

dunque, forse è questo che intendevi, hai una cancellazione polo-zero, e quindi x(t) e y(t) sono effettivamente uguali, MA a meno del rumore e delle condizioni iniziali. Le condizioni iniziali, in particolare, fanno sì che x ed y siano diverse per un certo transitorio iniziale, ovvero sono uguali asintoticamente.

Nel caso generale

ovviamente questo non è vero.

Ci sono due modi per ricondursi alla soluzione che ti ho dato nel primo post. Il primo, più "tonto" è calcolare la risposta all'impulso di H(z) il quale, avendo ora un denominatore, ha durata infinita (è un IIR, Infinite Impulse Response). Ma a noi interessa cosa accade nei primi N istanti, quindi possiamo calcolarci solo quelli.

Quindi, posto B= - a e A = - a, si calcola la risposta all'impulso e la sua trasformata di Fourier

%%%%%%%% h_resp = impz(B,A,N); % N è la lunghezza di x e y H_resp = fft(h_resp); %%%%%%%

oppure, in un'istruzione sola

%%%%%%% H_resp = freqz(B,A,N,'whole'); %%%%%%%

Dopodichè procedi come prima: X= Y./H_resp .

Un altro modo è notare che, siccome conosci y(t), puoi conoscere anche

y(t) + A y(t-1)

(ancora una volta, se non conosci le condizioni iniziali devi porre y(-1)=0 )

quindi ora hai di nuovo un problema simile a quello iniziale

dove y'(t)= y(t) + A y(t-1)

ed H'(z)= 1 + B z^-1,

quindi ancora

y'(t) = x(t) + B x(t-1),

MA SOLO perchè conosci già i valori di y(t) e perchè sai che y' è una funzione di y(t) ED y(t-1) !!!

Quindi puoi scrivere y'(t) in quel modo, ma solo allo scopo di deconvolvere y(t), non per calcolare y(t) in funzione di x(t)

Era questo che intendevi?

Ciao

Pasu

"Pasu" ha scritto nel messaggio news:gcib9o$e5p$ snipped-for-privacy@aioe.org... [CUT]

ok quindi non riuscirò mai a separare i due segnali asintoticamente, è brutto saperlo, ma almeno ora lo so e non ci perdo più tempo.

Grazie di tutto

ciao

Beh dai, non la mettere giù così pesante... A seconda del valore della tua 'a' il transitorio iniziale potrebbe esaurirsi in pochi campioni (esaurirsi=diventare più piccolo del rumore o comunque della precisione numerica della macchina).

Ma cosa devi fare esattamente?

Ciao

Pasu