Décomposition en somme de sinusoides

[invite36c8b572]

Bonjour à tous.

Afin d'analyser des données issus d'un radar"step frequency" (je ne sais pas si ça parle à beaucoup de monde mais c'est pas si grave pour la suite), je l’intéresse à la décomposition de signaux en somme de sinusoïdes.

Voici mon problème:
Un radar envoie des ondes dans le sol. Chaque reflections sur des interfaces vas être analysée par le récepteur et va produire une exponentielle de fréquence f, avec f proportionnelle à la distance de l'interface, sur le signal de sortie.
On obtient en sorti du radar un signal discrétisé de la forme :
s= sum(ai*exp(2*j*pi*fi*t)) + bruit

De ce signal on cherche à extraire les paramètres des exponentielles, afin de reconstituer la composition du sol. Une fréquence fi est proportionnelle à la distance de l'interface i et l'amplitude ai à sa nature.

L'opération standard est d'effectuer une fft. Cette opération permet de remonter aux paramètres des exponentielles mais est assez imprécise. Les pics s'entrecroisent entre les différentes fréquences, et la technique du zero-padding introduit des rebonds qui peuvent être mal interprétés.

C'est pourquoi je cherche des méthodes de décompositions plus précises pour remonter aux paramètres ai et fi.
Il semblerais que la méthode de Prony réponde à la question, pourtant j’obtiens des résultats étonnamment mauvais en utilisant cet algorithme sur matlab.

Voila, si vous avez des idées, des explications ou des méthodes de décomposition en exponentielles je vous attends!

Merci
-----

[obi76]

Salut,

juste une idée comme ça (ce n'est pas mon domaine) : mais faire comme les radars automatiques, I.E. faire une autocorrélation du signal émit et à recevoir, ça permet de savoir exactement quand tu le reçois. Il y a sans doute moyen de tirer des avantages de cette piste...

[invite36c8b572]

Il s'agit ici de radar un petit peu particulier (step frequency radar), le signal de sortie du radar ne correspond pas au signal reçu par l'antenne de réception. Le radar émet des ondes sinusoïdales en continue en faisant varier progressivement leur fréquence, puis il traite le signal reçu afin de nous transmettre un signal de la forme s= sum(ai*exp(2*j*pi*fi*t)) + bruit comme j'explique dans mon post précédent (t prend des valeurs discrètes). Mais le fonctionnement du radar n'a pas besoin d'être compris, il faut juste comprendre la forme du signal s.

[obi76]

Effectivement, dans ce cas je ne peux vous être d'aucune aide, désolé...

[A voir en vidéo sur Futura]

[phuphus]

Bonjour EChevallier,

serait-il possible que tu partages avec nous quelques données ? Un signal typique, ainsi que le code Matlab associé seraient les bienvenus, sauf problème de confidentialité. Sinon, un signal synthétisé par tes soins et représentatif de ta mesure, faisant apparaître ce problème de pics qui s'entrecroisent, serait OK.

Au pire, donne-nous des résultats connus en termes de ai et fi, et je recomposerai le signal dans mon coin pour voir ce que l'on peut en faire.

En attendant, quelques pistes :
- l'échantillonnage est-il bien fait avec un filtre anti-repliement ? Le filtre est-il correct ? La fréquence d'échantillonnage correctement choisie ?
- quelle appodisation utilises-tu ?
- je suppose que si des pics s'entrecroisent, alors tout ce qui est intéressant se passe dans une bande de fréquences réduite par rapport à ta bande passante totale. Dans ce cas, un zoom FFT pourrait être une solution. As-tu déjà regardé de ce côté ?

[phuphus]

... petits détails supplémentaires, mais qui à priori ne doivent pas interférer avec tes résultats :

- je suppose qu'en sortie du radar le signal est réel, donc s = real(sum(ai*exp(2*j*pi*fi*t))) + bruit
- une FFT est une décomposition en cosinus et non en sinus, mais encore une fois cela n'a que peu d'importance pour ton problème, puisqu'apparemment la phase ne t'intéresse pas

[invite36c8b572]

Je ne sais pas si c'est moi qui ne comprend pas mais je n'arrive pas à insérer d'images depuis mon ordinateur et non depuis un lien internet, ni à mettre de pièces jointes... ça va pas être facile de vous montrer! Si quelqu'un a une solution...

Le code matlab de base consiste juste à appliquer la fft à nos données.

La fft utilise des fréquences choisies arbitrairement et régulièrement espacées. Il y a apparitions de pics autour des fréquences effectivement présentes dans le signal. Dans l'idéal j'aimerais bien trouver une décomposition où les fréquences puissent prendre n'importe quelles valeurs réelles, et ainsi transformer les pics en raies. (C'est ce qu'est sensé faire Prony je crois)

[invite36c8b572]

La sortie du radar est un signal complexe . (combinaison de deux signaux)

[phuphus]

Re-bonjour EChevallier,

pour les pièces jointes, normalement il y a juste à se laisser guider par la fenêtre qui apparaît lorsque l'on cloque sur "Gérer les pièces jointes". Joindre un fichier sur ton dur se fait dans la partie la plus haute de cette fenpetre, avec 5 fichiers simultanés possibles.

Pour le code Matlab, pourrais-tu en faire une copie ici ? Pour plus de clarté, tu peux par exemple le publier dans ta réponse à l'aide de la balise CODE :

% Exemple de programme Matlab vide
clear all ;
clc ;
% Fin du programme Matlab vide

Les fréquences de la FFT ne sont pas vraiment arbitraires, puisqu'elles dépendent de la longueur de ton signal analysé. Le DeltaF est donc égal à Fech/Nfft. Donc, si jamais tu fais une analyse sur un nombre exact de périodes pour une fréquence considérée, tu trouveras une belle raie à cette fréquence et des 0 partout autour. Dans la pratique, cette situation n'est jamais rencontrée, on utilise donc des artifices pour améliorer un peu le résultat :
- Zero padding, que tu as déjà expérimenté (mais peut-être pas correctement...)
- Apodisation, pour limiter les effets de bord
- Zoom FFT, je pense la technique qui te serait la plus utile. Regarde par exemple ici p. 27, et dis-moi ce que tu en penses :

http://www.altracustica.org/docs/fr_analyse_sig_sys.pdf

Par contre, j'ai vraiment besoin d'un signal typique et de ton code Matlab pour pouvoir aller plus loin. L'application d'une FFT n'est jamais triviale, et la qualité des résultats dépend de beaucoup de facteurs annexes.

Es-tu allé voir du côté des méthodes avec modèle ? Comme tu connais à priori la forme de ton signal, alors le modéliser par une équation (comme celle que tu as écrite dans ta première intervention) et effecteur une régression pour trouver le signal modélisé approchant le plus ton signal réel doit pouvoir marcher (c'est d'ailleurs l'esprit de la méthode de Prony). Voir par exemple ici :

http://perso.telecom-paristech.fr/~r...Atiam-Chap.pdf

La FFT est surtout intéressante sur l'analyse de signaux dits arbitraires. Ici, tu connais déjà le type de signal que tu as.

[phuphus]

La sortie du radar est un signal complexe . (combinaison de deux signaux)

Y-a-t-il une justification physique à cela ??

[phuphus]

La fft utilise des fréquences choisies arbitrairement et régulièrement espacées. Il y a apparitions de pics autour des fréquences effectivement présentes dans le signal. Dans l'idéal j'aimerais bien trouver une décomposition où les fréquences puissent prendre n'importe quelles valeurs réelles, et ainsi transformer les pics en raies. (C'est ce qu'est sensé faire Prony je crois)

Pour une décomposition arbitraire, tu peux reprendre l'algorithme de la DFT et te le palucher en introduisant tes propres fréquences. Ici, les effets de bord seront encore plus visibles que dans la FFT, puisque tu verras directement apparaître les annulations du spectre de ton apodisation, mais pourquoi pas.

Sinon, tu peux aussi utiliser une astuce qui vaut ce qu'elle vaut. L'idée est de faire un filtre pour la seule fréquence qui t'intéresse, donc de faire un produit de convolution de ton signal par la sinusoïde considérée, et de prendre la valeur RMS du résultat. En faisant cela sur un ensemble de fréquences autour de tes pics, tu devrais pouvoir mettre en exergue un maximum à une fréquence beaucoup plus précise qu'avec du Fourier.

Je viens de tester les deux méthodes en Matlab, pour détecter une fréquence de 103.7 Hz sur un signal échantillonné à 10 kHz. En prenant du Fourier avec DeltaF de 10Hz (Nfft = 1000), le maximum est comme attendu à 100Hz. Avec la DFT personnalisée (raffinage entre 90 et 110Hz avec DeltaF = 0.1Hz), le pic est à 103.8Hz. Avec la convolution par des sinusoïdes choisies, le max se trouve à 103.6Hz.

Je te mets le code Matlab correspondant, désolé pour le manque de clarté et l'absence de commentaires mais je l'ai fait à la va-vite :

% Essai basique de détection précise de fréquence sur signal connu

clear all ;
clc ;

% Constantes
  Fs = 10000 ;
  Nyquist = round(Fs/2) ;
  j = sqrt(-1) ;

% Génération du signal
  F = 103.7 ;
  A = 1 ;
  Duree = 0.1 ;
  t = [0:1/Fs:Duree-1/Fs] ;
  s = real(A*exp(2*j*pi*F*t)) ;
  Taille = size(s,2) ;

% Fourier
  S = fft(s) ;
  Freq = [0:1/Duree:Nyquist] ;

% Filtrage sinus
  Fmin = 90 ;
  Fmax = 110 ;
  Resolution = 0.1 ;
  Resultats_SinFilt = zeros(2,((Fmax-Fmin)/Resolution)+1) ;
  indice = 0 ;
  for ii = Fmin:Resolution:Fmax
      indice = indice + 1 ;
      Resultats_SinFilt(1,indice) = ii ;
      temp = cos(2*pi*ii*t) ;
      temp = filter(temp,1,s) ;
      Resultats_SinFilt(2,indice) = sqrt((1/Taille)*sum(temp.^2));
  end

% DFT personnalisée
  Resultats_DFT = zeros(2,((Fmax-Fmin)/Resolution)+1) ;
  indice = 0 ;
  for ii = Fmin:Resolution:Fmax
      indice = indice + 1 ;
      Resultats_DFT(1,indice) = ii ;
      temp = exp(-j*2*pi*ii*t) ;
      Resultats_DFT(2,indice) = sum(s.*temp) ;
  end

% Tracés
  figure(1) ;
  plot(Freq,abs(S(1:501))/max(abs(S(1:501))),'b') ;
  hold on ;
  plot(Resultats_SinFilt(1,:),Resultats_SinFilt(2,:)/max(Resultats_SinFilt(2,:)),'r') ;
  plot(Resultats_DFT(1,:),abs(Resultats_DFT(2,:))/max(abs(Resultats_DFT(2,:))),'g') ;
  hold off ;

[invite36c8b572]

Tout d'abord merci de ton aide!

Les signaux issus du radar sont complexes, je crois que c'est qu'en fait le radar envoie 2 sinusoïdes en quadrature dans le sol. Cependant, je m'y connait très peu en électronique, et n'ai en fait jamais vu le radar. Je travaille sur des données qu'on me fournit.

Pour le zoom fft, je n'ai pas encore tout à fait bien compris mais est ce que ça ne ressemble pas au zero-padding?

Le lien vers la page de telecom-paris me parait très intéressant ( c'est assez ardu , je vais y passer du temps).

Pour les pieces jointes:
Le fichier données est un tableau complexe à 2 dimensions qui contient la sortie du radar. Chaque colonne contient une sortie. Les différentes colonnes correspondent au déplacement du radar sur le sol

Le fichier image contient, plot du module des données, plot de la fft, plot de la fft avec zero padding, et 2 images représentant le résultat final qui est en fait constitué d'un acollement de plusieurs tableaux de données.

J'obtient le zero padding "manuellement" en ajoutant des zeros aux données:

data_cmx=[data_cmx ; zeros(size_ifft-size(data_cmx,1),d)];
data_cmx_time=ifft(data_cmx,size(data_cmx,1), 1);

Oui, en fait il s'agit d'une ifft, mais cela revient presque au même (je n'ai pas encore compris pourquoi je dois appliquer une ifft et non une fft au vu de la forme du signal s).

Voila les opérations classiques pour lire les données

[invite36c8b572]

J'ai en fait essayer de programmer mon algorithme de décomposition:

Le but est de trouver les fréquences optimales de décomposition d'un signal
<f,g> est le produit scalaire de f et g
x un signal discret
On trouve la fréquence f0 tel que <x(.),exp(-2iPif0.)> soit maximum.
On calcule x1=x-<x(.),exp(-2iPif0.)>exp(-2iPif0.)
On a maintenant la première composante du signal (f0,<x(.),exp(-2iPif0.)>)
On recommence avec le signal x1...
On continue jusqu’à ce que le signale xn soit égale à zéro ou presque.
A la fin on obtient une décomposition du signale sur les fréquences f(n).
A chaque étapes on a choisit la fréquence optimale.

Les résultats ne sont pas mauvais, plutôt intéressant même, le problème c'est que le signal xn n'atteint pas zéros, et que les exponentielles trouvées ne sont pas forcement indépendantes. C'est pourquoi je cherche une méthode un peu plus propre .

"On trouve la fréquence f0 tel que <x(.),exp(2iPif0.)> soit maximum. " En fait pour trouver ce maximum j'utilise le zero padding pour augmenter la résolution et je prend le max du zero padding .

J'ai joint quelque images de cet algo et le code de la décomposition.
Il prend en etré un signal x , un nombre N_padd qui est la uniquement pour la précision de recherche du maximum et un nombre a d’itération de l'algorithme.

Je te remercie Phuphus pour ton intérêt et ton aide.

[phuphus]

Bonjour EChevallier,

merci pour tes données, je regarde cela dès que possible.

Concernant le produit scalaire, cela revient au même que ce que j'ai appelé précédemment "filtrage sinus"
Tu trouveras quelques précisions ici :

http://www.cours.polymtl.ca/ele3700/...u_syn_chII.pdf

[invite36c8b572]

Petite erreur dans la première pièce jointe, il manquait les images.

[phuphus]

Bonjour EChevallier,

impossible d'ouvrir tes fichiers .fig, je pense que c'est dû à ma version de Matlab qui est trop ancienne (6.1). Le mieux serait vraiment que je puisse avoir les données brutes et le programme Matlab les exploitant, sauf bien sûr si cela pose problème (mais comme tu as pu diffuser les .fig, je suppose que ça sera OK pour les données sources).

J'obtient le zero padding "manuellement" en ajoutant des zeros aux données:
Code:

data_cmx=[data_cmx ; zeros(size_ifft-size(data_cmx,1),d)];
data_cmx_time=ifft(data_cmx,size(data_cmx,1), 1);

Oui, en fait il s'agit d'une ifft, mais cela revient presque au même (je n'ai pas encore compris pourquoi je dois appliquer une ifft et non une fft au vu de la forme du signal s).

En voyant le code, j'aurais tendance à dire que :
- data_cmx c'est du fréquentiel
- data_cms_time, c'est du temporel
- le zero padding est réalisé sur le fréquentiel, c'est donc pour faire du suréchantillonnage

Me trompe-je ?

Du coup, cela expliquerait peut-être pourquoi le radar renvoie deux colonnes (amplitude / phase ou partie réelle / partie imaginaire). C'est d'ailleurs facile pour lui, si l'excitation n'est pas un chirp mais un balayage "lent", alors obtenir une analyse fréquentielle directe sans faire appel à une décomposition quelconque est simple : analyse du max du signal de retour pour avoir l'amplitude à la fréquence d'excitation courante, et analyse du déphasage entre max de l'excitation et max local du signal reçu pour le déphasage. Pour que cela fonctionne, il faut que le temps de retour du signal soit petit devant le temps caractéristique d'évolution des fréquences du balayage. Est-ce le cas ?

Pour expliquer le truc autrement : tu traces le signal temporel de retour, l'enveloppe de ce signal représente de manière directe une courbe amplitude / fréquence puisque ce signal est issu d'un balayage en fréquence. Mais cela ne fonctionne que si le balayage est "lent" voire discret (d'où le nom de step frequency, non??), sinon plusieurs fréquences peuvent revenir en même temps.

[phuphus]

... du coup, je viens de faire une recherche vite fait. La partie gauche du diagramme page 131 du document suivant correspond à ton matos ?

http://www.jpier.org/PIER/pier28/06....1.W.Chew.M.pdf

Si oui, alors les résultats issus du radar sont bien des mesures fréquentielles et non un signal temporel imaginaire.

Mais le fonctionnement du radar n'a pas besoin d'être compris, il faut juste comprendre la forme du signal s

Les deux sont liés, mon capitaine !

[invite36c8b572]

Personnellement je ne connais pas du tout le radar en question, et je ne m'y connais pas vraiment dans ce genre d'appareils. Je traite juste des données qui en sont issues (oui oui je suis un pauvre stagiaire qui ne reçoit pas beaucoup d'aide ).

Je pense que ce que tu dis sur le fonctionnement du radar est juste. Les fréquences sont incrémentées de manières discrètes.
Cependant la sortie "finale" du radar est un peu plus complexe. Le signal que tu décrit subit un traitement avant de "sortir" du radar.

Simplifions d'abord et considérons le cas réel.
Supposons que l'on ai une réflexion sur une unique cible (les choses étant ensuite linéaire...).
Comme tu le dis, pour chaque fréquence fi, le radar récupère une amplitude Ai et une phase Phasei. Mais ce qu'il nous donne c'est : Ai*cos(Phasei). On obtient donc un tableau de données réelles contenant les Ai*cos(Phasei). Les fréquences varient linéairement. On a en fait Phasei=2*pi*(f0+i*deltaf)*temp s_propagation. Les Ai sont à peu près égaux, la courbe obtenue est donc un cosinus dont la fréquence est proportionnelle au temps de propagation. Le but est donc de retrouver la fréquence du cosinus pour retrouver le temps de propagation.

Les signaux sont complexe car je crois que le radar envoie 2 sinusoïdes en quadratures, et ressort un couple (Ai*cos(phasei),Ai*sin(phasei) ), soit Ai*exp(phasei).

C'est vrai comprendre le fonctionnement du radar est forcément utile à l’interprétation des données, mais initialement je pensais juste demander de l'aide sur les décompositions, pour que les gens ne soient pas obliger de rentrer dans la partie technique.