Calculer et interpréter la fréquence instantanée

Je suis nouveau dans le principe du calcul de la fréquence instantanée et j'ai posé beaucoup de questions à ce sujet. Vous les trouverez tous dans une liste à puces à la fin de ce texte. Le texte peut être un peu long, excusez-moi, mais j'ai vraiment essayé de résoudre ce problème par moi-même.

Je m'intéresse donc à la fréquence instantanée d'un signal de valeur réelle . Le calcul est effectué à l'aide d'un signal analytique , où est la transformation de Hilbert de .$f(t)$$x(t)$$z(t) = x(t) + j y(t)$$y(t)$$x(t)$

Pour calculer les fréquences instantanées à partir du signal analytique j'ai suivi l'article:$z(t)$

Le calcul de la fréquence instantanée et de la bande passante instantanée par Arthur E. Barns à partir de 1992. Dans cet article, il présente plusieurs méthodes pour calculer la fréquence instantanée. J'écris, toutes les formules qu'il a proposées (et j'ai utilisées) dans un instant.

Pour "apprendre", j'ai joué avec un signal très simple, et deux signaux un peu plus complexes, dans MATLAB, et je voulais obtenir leurs fréquences instantanées.

Fs = 1000;                                            % sampling-rate = 1kHz
t = 0:1/Fs:10-1/Fs;                                    % 10s 'Timevector'
chirp_signal = chirp(t,0,1,2);                         % 10s long chirp-signal, signal 1
added_sinusoid = chirp_signal + sin(2*pi*t*10);        % chirp + sin(10Hz), signal 2
modulated_sinusoid = chirp_signal .* sin(2*pi*t*10);   % chirp * sin(10Hz), signal 3

Les tracés dans le domaine temporel de ces trois signaux se présentent comme suit:

Les tracés de toutes les fréquences instantanées que j'ai obtenues après avoir appliqué toutes les méthodes du papier sont les suivants:

Fréquences instantanées du signal gazouillis pur: Fréquences instantanées du signal gazouillis avec sinusoïde ajoutée: Fréquences instantanées du signal gazouillis modulé: Veuillez noter que, dans les trois images, l'axe y des tracés 3 et 4 sont agrandis, de sorte que les amplitudes de celles-ci les signaux sont très petits!

La première possibilité pour passer du signal analytique à la fréquence instantanée est: où est la phase instantanée. Je pense que c'est la méthode la plus couramment utilisée aujourd'hui, du moins sur la page Web de MATLAB, elle est calculée de cette façon. Le code ressemble à ceci:

$$\begin{equation} f_{2}(t) = \frac{1}{2\pi} \frac{d}{dt}\theta(t) \end{equation}$$ $\theta(t)$

function [instantaneous_frequency] = f2(analytic_signal,Fs)
    factor =  Fs/(2*pi);
    instantaneous_frequency = factor * diff(unwrap(angle(analytic_signal)));
    % Insert leading 0 in return-vector to maintain size
    instantaneous_frequency = [0 instantaneous_frequency];
end

Dans l'article, Barns suggère maintenant (ou plutôt dit compile) quatre autres façons de calculer les fréquences instantanées à partir du signal analytique. Il mentionne également la formule supérieure, mais est d'avis qu'elle n'est pas pratique en raison des ambiguïtés de la phase. Je suppose qu'il ne connaissait pas la unwrap()méthode ou, pour être plus précis, les mathématiques derrière. (J'ai moi-même appris cette méthode aujourd'hui seulement, en regardant d'autres codes sources sur les fréquences instantanées)

Dans son article, la formule a l'étiquette Numéro (2), par conséquent, j'ai donné au f (t) l'index 2. Tous les autres index correspondent de la même manière à leurs nombres dans le papier.

Du fait des ambiguïtés en phase, il suggère plutôt:

$$\begin{equation} f_{3}(t) = \frac{1}{2\pi}\cdot\frac{\overbrace{x(t)y'(t)}^\text{a}-\overbrace{x'(t)y(t)}^{\text{b}}}{\underbrace{x(t)^2}_{\text{c}}+\underbrace{y(t)^2}_{\text{d}}} \end{equation}$$ J'ai introduit les symboles "a", "b", "c" et "d", pour faciliter un peu la programmation:

function [instantaneous_frequency] = f3(analytic_signal,Fs,T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    diff_x = diff(x);
    diff_y = diff(y);
    factor = Fs/(2*pi);
    a = x(2:end).*diff_y;
    b = y(2:end).*diff_x;
    c = x(2:end).^2;
    d = y(2:end).^2;
    instantaneous_frequency = factor * ((a-b)./(c+d));
    % Insert leading 0 in return-vector to maintain size
    instantaneous_frequency = [0 instantaneous_frequency];
end

Puis Barner donne trois autres formules qu'il nomme "approximations instantanées de fréquence":

$$\begin{equation} f_{9}(t) = \frac{1}{2\pi T}\cdot\text{arctan}\left[\frac{\overbrace{x(t)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t)}^{\text{b}}}{\underbrace{x(t)x(t+T)}_{\text{c}}+\underbrace{y(t)y(t+T)}_{\text{d}}}\right] \end{equation}$$

function[instantaneous_frequency] = f9(analytic_signal, Fs, T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = Fs/(2*pi*T);
    a = x(1:end-T).*y(1+T:end);
    b = x(1+T:end).*y(1:end-T);
    c = x(1:end-T).*x(1+T:end);
    d = y(1:end-T).*y(1+T:end);
    instantaneous_frequency = factor.*atan((a-b)./(c+d));
    % Append 0 to return-vector to maintain size
    instantaneous_frequency = [instantaneous_frequency zeros(1,T)];
end

$$\begin{equation} f_{11}(t) = \frac{1}{4\pi T}\cdot\text{arctan}\left[\frac{\overbrace{x(t-T)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t-T)}^{\text{b}}}{\underbrace{x(t-T)x(t+T)}_{\text{c}}+\underbrace{y(t-T)y(t+T)}_{\text{d}}}\right] \end{equation}$$

function [instantaneous_frequency] = f11(analytic_signal, Fs, T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = Fs/(4*pi*T);
    a = x(1:end-2*T).*y(1+2*T:end);
    b = x(1+2*T:end).*y(1:end-2*T);
    c = x(1:end-2*T).*x(1+2*T:end);
    d = y(1:end-2*T).*y(1+2*T:end);
    instantaneous_frequency = factor.*atan((a-b)./(c+d));
    % Append and insert 0s to maintain size
    instantaneous_frequency = [zeros(1,T) instantaneous_frequency zeros(1,T)];
end

$$\begin{equation} f_{14}(t)=\frac{2}{\pi T}\left[\frac{\overbrace{x(t)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t)}^{\text{b}}}{\underbrace{(x(t)+x(t+T))^2}_{\text{c}}+\underbrace{(y(t)+y(t+T))^2}_{\text{d}}}\right] \end{equation}$$

function [instantaneous_frequency] = formula14(analytic_signal, Fs, T);
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = 2*Fs/(pi*T);
    a = x(1:end-T).*y(1+T:end);
    b = x(1+T:end).*y(1:end-T);
    c = (x(1:end-T)+x(1+T:end)).^2;
    d = (y(1:end-T)+y(1+T:end)).^2;
    instantaneous_frequency = factor * ((a-b)./(c+d));
    % Append and insert 0s to maintain size
    instantaneous_frequency = [instantaneous_frequency zeros(1,T)];
end

Dans les 3 approximations, les formules T ont été définies sur Fs (T = Fs = 1000 = 1s), comme suggéré dans le document.

Maintenant, ma question est:

• Les formules f2 et f3 renvoient le même résultat pour le signal chirp pur. Je pense que c'est bien, car ils calculent la même chose. Les trois méthodes d'approximation ne renvoient pas la même chose, pas même quelque chose qui s'en rapproche! Pourquoi est-ce le cas? (J'espère que ce n'est pas seulement un bug de programmation ...)
• Bien qu'ils reviennent de la même manière, en particulier à la fin de l'intrigue, ils commencent à beaucoup se tortiller . Quelle en est l'explication? J'ai d'abord pensé à quelque chose comme le repliement, mais ma fréquence d'échantillonnage est assez élevée, par rapport à la fréquence des signaux, donc je pense que cela peut être exclu.

• Au moins f2 et f3 semblent fonctionner correctement sur un signal de chirp pur, mais toutes les méthodes, y compris f2 et f3 semblent échouer horriblement, quand il s'agit de plus d'une fréquence dans le signal. En réalité, avoir plus d'une fréquence dans un signal est plutôt toujours le cas. Alors, comment obtenir la fréquence instantanée (plus ou moins) correcte?

  • En fait, je ne sais même pas à quoi m'attendre, lorsque plus d'une fréquence est présente dans le signal. Le calcul renvoie un nombre pour un moment donné, alors que doit-il faire lorsque, comme ici, plus de fréquences sont présentes? Renvoyer la moyenne de toutes les fréquences ou quelque chose comme ça?

• Et ma question probablement la plus importante est: comment cela est-il géré dans un logiciel réel et élaboré? Disons que je veux connaître la fréquence instantanée du signal modulé à 1,75 s, et j'ai choisi la méthode f2, que je puisse être «chanceux» et obtenir un nombre proche de 6 [Hz] qui est probablement la bonne réponse, ou je choisir mes résultats quelques échantillons à côté et soudain, j'obtiens un résultat filaire, trop élevé, car j'ai malheureusement choisi une valeur dans le pic. Comment cela peut-il être géré? En le post-traitant avec un filtre médian moyen ou mieux? Je pense que même cela pourrait devenir vraiment difficile, surtout dans les régions où de nombreux pics sont côte à côte.

Et une dernière question, pas si importante, pourquoi la plupart des articles que je trouve sur les fréquences instantanées proviennent du domaine de la géographie, en particulier dans le calcul des événements sismographiques comme les tremblements de terre. Le document de Barne prend également cela comme exemple. La fréquence instantanée n'est-elle pas intéressante dans de nombreux domaines?

Ça y est jusqu'à présent, je suis très reconnaissant pour chaque réponse, surtout quand quelqu'un me donne des conseils sur la façon de l'implémenter dans un vrai projet logiciel ;)

Cordialement, Patrick

Pas vraiment une réponse mais peut-être utile: personnellement, j'ai trouvé que le concept de fréquence instantanée n'est utile que pour des signaux à bande suffisamment étroite.

Prenons l'exemple simple de deux ondes sinusoïdales stables, disons 100 Hz et 934 Hz. Dans ce cas, vous pouvez certainement définir et calculer la fréquence instantanée (comme vous le souhaitez), mais quel devrait être le résultat? Quelles informations ou propriétés possibles la fréquence instantanée peut-elle avoir pour dire quelque chose de significatif sur le signal? Appliquer le concept de fréquence instantanée à des signaux qui ont plusieurs fréquences en même temps n'a tout simplement pas beaucoup de sens.

C'est pourquoi vous obtenez des résultats décents pour les balayages mais des courbes bizarres pour le sinus Sweep +. C'est aussi la raison pour laquelle vous voyez les ondulations dans la partie haute du balayage. La bande passante du signal devient trop élevée pour lui attribuer un numéro de fréquence unique et donc les résultats sautent.

Au moins f2 et f3 semblent fonctionner correctement sur un signal de chirp pur, mais toutes les méthodes, y compris f2 et f3 semblent échouer horriblement, quand il s'agit de plus d'une fréquence dans le signal. En réalité, avoir plus d'une fréquence dans un signal est plutôt toujours le cas. Alors, comment obtenir la fréquence instantanée (plus ou moins) correcte?

comme le suggère Hilmar, la méthode de la transformée de Hilbert (ou "signal analytique") ne fonctionne pas sur la large bande car il existe plusieurs composantes de fréquence. vous pouvez faire cette méthode uniquement pour un seul composant sinusoïdal.

donc, avec l'approche du signal analytique, ce que vous voulez faire, c'est utiliser cette identité:

$$\arctan u - \arctan v = \arctan \left( \frac{u-v}{1+uv} \right)$$

$|u-v|$$f_9$

mais il ne doit y avoir qu'une seule sinusoïde variant dans le temps dans le calcul de la transformée de Hilbert pour le faire correctement. et vous feriez mieux d'aligner le composant "en phase" avec la sortie de la transformée de Hilbert (qui est retardée avec un filtre FIR causal). sinon vous aurez de la merde.

Wow, quelle énorme question. Je vais d'abord répondre à la question pas si importante:

Et une dernière question, pas si importante, pourquoi la plupart des articles que je trouve sur les fréquences instantanées proviennent du domaine de la géographie, en particulier dans le calcul des événements sismographiques comme les tremblements de terre. Le document de Barne prend également cela comme exemple. La fréquence instantanée n'est-elle pas intéressante dans de nombreux domaines?

La raison en est que le système sismographique "vibroseis" est utilisé dans l'industrie pétrolière pour effectuer des levés sismiques. Les camions que j'ai reliés vibrent d'environ 5 Hz à environ 90 Hz et peuvent être conçus pour émettre des signaux de gazouillis. Il y a beaucoup d'argent dans l'industrie pétrolière et le traitement des retours de ces signaux peut être très, très lucratif. Par conséquent, de nombreuses personnes ont passé de nombreuses heures à analyser de tels signaux, notamment à étudier des techniques de fréquence instantanée.

---

$^{\rm TM}$

Consultez cet article.

De meilleures approches ont tendance à utiliser des «moyennes pondérées par phase», telles qu'elles sont mises en œuvre ici . Ou ici pour un lien direct vers le matlab .

Désolé de fournir une réponse un an après le fait, mais je suis tombé sur ce post en cherchant des articles sur ce même sujet. Vos questions reflètent les désaccords et interprétations répandus de la "fréquence instantanée" qui ont tourmenté le domaine depuis sa création. De nombreuses personnes vous diront, comme certaines des réponses ici, que l'IF ne s'applique qu'aux signaux "à bande étroite" ou "mono-composante". En fait, ce n'est pas vrai: parfois le FI obtenu par la transformée de Hilbert se comporte parfaitement bien pour des signaux large bande et / ou "multi-composants". Une quantité qui a été proposée pour éviter bon nombre de ces difficultés est la "fréquence instantanée moyenne pondérée (WAIF)", qui peut être mesurée à l'aide d'un spectrogramme.

Voir Loughlin dans J. Acoust. Soc. Am., Janvier 1999. Picinbono (IEEE Trans. Sig. Proc., Mars 1997) et Vakman (IEEE Trans. Sig. Proc., Avril 1996) ont publié d'autres bons articles sur les FI et les idées fausses courantes.