Re: Fourier ed analisi dei segnali

From: Paolo Russo <paolrus_at_libero.it>
Date: Mon, 03 Nov 2003 16:50:48 GMT

[Drizzt:]
>A0 = 1/ND * E [y(I)]
>
>Bn = 2/ND * E [ y(I) * sin(n * I * 2 * PI / ND)
>Cn = 2/ND * E [ y(I) * cos(n * I * 2 * PI / ND)
>
>An = SQR( (Bn*Bn) + (Cn*Cn) )
>Fn = arctg(Cn/Bn)
>
>Ok, non funziona :-)
>No so' se i calcoli di cui sopra necessitano di qualche adattamento
>particolare. Io ho provato ad analizzare una sequenza di punti ottenuta
>campionando 360 punti di cos(x) con x che va da 0 a 6.28 Se sintetizzo
>l'armonica 1, ottengo il risultato di partenza sfasato di 180 gradi...

T'e` andata ancora bene... :-)
Dovresti aver ottenuto B1=0 (circa) e C1=1. Gia` nel
calcolare C1/B1 avresti dovuto ottenere una divisione per
zero, con risultato dipendente dall'hardware e software che
usi. Dato che, a causa degli errori di calcolo (e dell'aver
usato 6.28 anziche' 2*PI) probabilmente B1 non e` venuto
proprio esattamente zero, il rapporto C1/B1 puo` esserti
venuto un numero spropositato con segno a caso. Con segno
negativo, arctg(-oo)=-PI/2; con segno positivo,
arctg(oo)=PI/2. Come vedi, ci sono 180 gradi che vanno e
vengono come niente fosse...
Il problema e` intrinseco nel fatto di calcolare la fase in
quel modo. Gia` nel solo calcolo del rapporto Cn/Bn perdi un
bit di informazione sui segni; il rapporto e` uguale se i
segni sono, ad esempio entrambi positivi o entrambi negativi,
ma la fase e` ben diversa...
Teoricamente, la fase andrebbe calcolata nel seguente
abominevole modo (in pseudo C):

if (Bn>0) Fn=atan(Cn/Bn);
else if (Bn<0) Fn=atan(Cn/Bn)+PI;
else if (Cn<0) Fn=-PI/2;
else Fn=PI/2;

Ma questo in teoria; e` quello che potrebbe scrivere un
matematico puro, ma un programmatore dovrebbe preoccuparsi
degli overflow e, in generale, degli algoritmi instabili. Di
per se' la dipendenza della fase da Cn e Bn e` stabile,
tranne quando entrambi Cn e Bn tendono a zero, nel qual caso
la fase e` mal definita ma irrilevante. Altre instabilita`,
insite negli algoritmi ma non nella natura del problema, sono
gratuite e vanno evitate. Meglio allora scambiare, se
necessario, Cn e Bn per stabilizzare il rapporto:

double Phase(Cn,Bn)
    double Cn,Bn;
{
if (Bn<0) return PI-Phase(Cn,-Bn);
if (Cn<0) return -Phase(-Cn,Bn);
if (Cn>Bn) return PI/2-Phase(Bn,Cn);
if (Bn==0) return 0.0;
return atan(Cn/Bn);
}

Pero`, se fossi in te, manterrei la rappresentazione (Bn,Cn)
se dovessi rigenerare il segnale per qualunque scopo che non
fosse di semplice verifica.

Ciao
Paolo Russo
Received on Mon Nov 03 2003 - 17:50:48 CET