MATLAB bruit blanc bizarre

[invitea2e53836]

Bonjour,

je cherche à générer un bruit blanc de fréquence avec matlab.

Voilà ce que je fais :

L = 10000;
sigma = 1;
V = sqrt(sigma);
S = V*randn(1,L);
Pxx = pwelch(S);
loglog(Pxx);

en faisant le loglog j'obtiens ce graphe :



je suis assez étonné que la psd ne soit pas centrée sur 1 (10^0). En effet un mean(Pxx) donne environ 0.32.
Le carré de mon V (1 donc) ne devrait t'il pas me donner le niveau autour duquel la psd va se centrer c'est à dire 1 et non 0.32 ?
Quelqu'un voit mon erreur (de code ou de raisonnement) ?

merci d'avance.
-----

[inviteb9f49292]

Je ne connais pas du tout pwelch, et je n'ai jeté qu'un coup d'oeil rapide, donc ne prends pas la suite pour argent comptant...

-la moyenne calculée de Pxx se retrouve bien sur le graphe, pourquoi voudrais tu que le graphe affiche autre chose que Pxx ?
-la densité spectrale de puissance est... une densité spectrale, il faut donc normaliser la la fréquence d'échantillonnage pour retrouver le niveau absolue.
-un autre moyen de vérification est "parseval", la variance de ton signal temporel doit "se retrouver" dans sa densité spectrale de puissance, pour du bruit blanc sa variance est identique a l'intégrale de son spectre (somme cumulée des canaux fréquentiels)
-pour la psd de longueur "n" j'utilise

s_hat = fft (s, n) / n;
psd = s_hat .* conj(s_hat);

Attention au repliement des fréquences négatives

Je pense que ces deux remarques devraient te mettre sur la voie...

[invitea2e53836]

-la moyenne calculée de Pxx se retrouve bien sur le graphe, pourquoi voudrais tu que le graphe affiche autre chose que Pxx ?

oui oui pas de soucis là dessus, c'était juste pour illustrer quantitativement ce que le graphe affichait.

la densité spectrale de puissance est... une densité spectrale, il faut donc normaliser la la fréquence d'échantillonnage pour retrouver le niveau absolue.

je n'ai pas trop compris ce que tu voulais dire par là, normaliser où et pourquoi ?

un autre moyen de vérification est "parseval", la variance de ton signal temporel doit "se retrouver" dans sa densité spectrale de puissance, pour du bruit blanc sa variance est identique a l'intégrale de son spectre (somme cumulée des canaux fréquentiels)

on est d'accord, c'est pour ça justement que je pensais de mémoire que le carré de la variance donne le niveau théorique de la psd pour le bruit blanc or ce n'est pas le cas vu que je n'ai pas 1.

Pour ta définition de la psd, je suppose que ça vient de la psd définie à partir de l'autocorrélation. Cependant Pwelch marche très bien selon moi, le soucis ne se situe pas au niveau du calcul de la psd je pense.

[invitea2e53836]

Up, j'ai tenté de tracer la psd avec :

L = 10000;
sigma = 1;
V = sqrt(sigma);
S = V*randn(1,L);

Noverlapp = 512;
Nfft = 1024;
Fs = 757;
[Pxx,F] = pwelch(S,kaiser(1024), Noverlapp, Nfft, Fs);
loglog(F/Fs,Pxx);

Même si j'ai pas compris à 100% le rôle de chaque paramètre (notamment le paramètre que prend la fenêtre kaiser et le Noverlapp).
En tout cas ma psd n'est toujours pas centrée sur 1, elle est juste décalée vers la gauche en basse fréquence :/

[A voir en vidéo sur Futura]

[inviteb9f49292]

je n'ai pas trop compris ce que tu voulais dire par là, normaliser où et pourquoi ?

C'est là que réside ton problème.

Pour ta définition de la psd, je suppose que ça vient de la psd définie à partir de l'autocorrélation. Cependant Pwelch marche très bien selon moi, le soucis ne se situe pas au niveau du calcul de la psd je pense.

Il n'y a qu'une seule définition de la densité spectrale de puissance, il y a par contre plusieurs estimateurs différents. Effectivement le soucis ne se situe pas dans le calcul de la psd mais dans l'interprétation de son résultat.

Même si j'ai pas compris à 100% le rôle de chaque paramètre (notamment le paramètre que prend la fenêtre kaiser et le Noverlapp).

vraisemblablement la PSD est faite par FFT de 1024 échantillons, avec un taux de recouvrement de 50%.

[invitea2e53836]

C'est là que réside ton problème.

J' peux avoir une explication ici ?

Sinon j'ai cherché un peu et pour la normalisation j'ai trouvé que l'intégrale du carré de mon bruit blanc doit être égale à N au carré. Le facteur de normalisation s'exprimant comme sqrt(2*ts/N).

Du coup j'ai fait ça :

%définition bruit blanc
L = 10000;
sigma = 1;
V = sqrt(sigma);
S = V*randn(1,L);

%calcul PSD
Noverlapp = 512;
Nfft = 1024;
Fs = 757;
Ts = 1/Fs;
[Pxx,F] = pwelch(S,kaiser(1024), Noverlapp, Nfft, Fs);

%normalisation et tracé
INT = max(cumtrapz(S));
N = sqrt(INT);
NORM = sqrt(2*Ts/N);
PxxN = Pxx/NORM;
loglog(F/Fs,PxxN);

Ceci me donne une PSD centrée sur 10^-1.
Du coup j'ai remarqué que si je la multiplie par 10 je tombe bien sur une psd centrée sur 1, sauf que cette dernière partie c'est du bidouillage et que je n'ai pas la justification pour ça

[invitea2e53836]

En fait j'ai trouvé mieux :

L = 10000;
sigma = 1;
V = sqrt(sigma);
S = V*randn(1,L);

%calcul PSD
Noverlapp = 512;
Nfft = 1024;
Fs = 757;
Ts = 1/Fs;
[Pxx,F] = pwelch(S,kaiser(1024), Noverlapp, Nfft, Fs);

%normalisation et tracé
INT = trapz(S.*S);
N = sqrt(INT);
NORM = sqrt(2*Ts/N);
PxxN = Pxx/NORM;
loglog(F/Fs,PxxN);
MOY = mean(PxxN)

du coup j'ai une PSD constamment centrée sur 0.5
Pourquoi un facteur 1/2 par rapport à ce que je m'attends à avoir ?

[inviteb9f49292]

J' peux avoir une explication ici ?

Si tu fais ta FFT sur N points échantillonnés à Fsr Hz, la "largeur" de ton échantillon fréquenciel est de Fsr/N Hz, et sa valeur est la puissance transportée pour cette fréquence, sur une largeur de Fsr/N Hz, que tu veux ramener à une valeur de puissance transportée sur une largeur de 1Hz (densité spectrale -> unité / Hz). Je pense que ton passage par l'intégrale et tout ça est une autre manière de dire la même chose...

Pourquoi un facteur 1/2 par rapport à ce que je m'attends à avoir ?

Car il y a de la puissance transportée dans les fréquences négatives...

[invitea2e53836]

Bon j'ai tenté un autre truc :

N = 10000;
sigma = 280;
f = rand(1,N);
f = sqrt(sigma)*(f-mean(f))/sqrt(var(f,1));
F = fft(f)/N;
DSP = F.*conj(F);

Il se trouve qu'avec ça en faisait un trapz(DSP) je trouve pile poil mon sigma=280. Sauf que du coup je suis perdu et y a un truc que j'ai complètement loupé d'ou la question suivante :

De quelle manière le coefficient sigma est répercuté sur l'allure de la PSD ?
Si je suis ce que je pensais jusqu'à maintenant, alors dans ce cas si les choses sont bien faites alors ça donne que la PSD est centrée sur la valeur du sigma, dans le cas que j'ai pris on serait sur .
Or, avec ce que je viens de faire, il se trouve que c'est l'intégrale de la DSP qui donne 280.

Autement dit : si je veux une PSD centrée sur alors que dois-je mettre comme valeur sur mon sigma ?

[invitea2e53836]

up ?

[invitea2e53836]

personne ?

[invitea2e53836]

euh y a toujours quelqu'un ?

[invitea2e53836]

help please

[Antoane]

Bonjour,

déplacé sur demande de l'auteur.

[invitea2e53836]

Ici non plus personne ?