2017-08-03 2 views
0

J'ai un instrument qui produit des données grossièrement sinusoïdales, mais dont la fréquence varie légèrement dans le temps. J'utilise MATLAB pour prototyper du code pour caractériser la dépendance temporelle, mais je rencontre des problèmes.Détermination de la fréquence dépendant du temps à l'aide d'une fenêtre FFT

Je générant une approximation idéalisée de mes données, I (t) = sin (pi 2 f (t) t), avec f (t) variable, mais actuellement testé comme linéaire ou quadratique. Je mettre en place ensuite une fenêtre de Hamming de glissement (de largeur w) pour générer un ensemble de transformées de Fourier F [I (t), t '] correspondant aux points de données dans I (t), et chaque F [I (t), t '] est compatible avec un Gaussien pour déterminer plus précisément l'emplacement du pic.

Mon code actuel de Matlab est:

fs = 1000; %Sample frequency (Hz) 
tlim = [0,1]; 
t = (tlim(1)/fs:1/fs:tlim(2)-1/fs)'; %Sample domain (t) 
N = numel(t); 

f = @(t) 100-30*(t-0.5).^2; %Frequency function (Hz) 
I = sin(2*pi*f(t).*t); %Sample function 

w = 201; %window width 
ww=floor(w/2); %window half-width 

for i=0:2:N-w 

    %Take the FFT of a portion of I, convolved with a Hamming window 
    II = 1/(fs*N)*abs(fft(I((1:w)+i).*hamming(w))).^2; 
    II = II(1:floor(numel(II)/2)); 
    p = (0:fs/w:(fs/2-fs/w))'; 

    %Find approximate FFT maximum 
    [~,maxIx] = max(II); 
    maxLoc = p(maxIx); 

    %Fit the resulting FFT with a Gaussian function 
    gauss = @(c,x) c(1)*exp(-(x-c(2)).^2/(2*c(3)^2)); 
    op = optimset('Display','off'); 
    mdl = lsqcurvefit(gauss,[max(II),maxLoc,10],p,II,[],[],op);  

    %Generate diagnostic plots 
    subplot(3,1,1);plot(p,II,p,gauss(mdl,p)) 
    line(f(t(i+ww))*[1,1],ylim,'color','r'); 

    subplot(3,1,2);plot(t,I); 
    line(t(1+i)*[1,1],ylim,'color','r');line(t(w+i)*[1,1],ylim,'color','r') 

    subplot(3,1,3);plot(t(i+ww),f(t(i+ww)),'b.',t(i+ww),mdl(2),'r.'); 
    hold on 
    xlim([0,max(t)]) 
    drawnow 
end 
hold off 

Mon processus de pensée est que l'emplacement de pointe dans chaque F [I (t), t '] devrait être une bonne approximation de la fréquence au centre de la fenêtre qui a été utilisée pour le produire. Cependant, cela ne semble pas être le cas, expérimentalement. J'ai eu un certain succès en utilisant l'analyse de Fourier discrète pour des problèmes d'ingénierie dans le passé, mais je n'ai fait que des cours sur les transformées de Fourier continues - donc il peut y avoir quelque chose d'évident qui me manque. Aussi, c'est ma première question sur StackExchange, donc les critiques constructives sont les bienvenues.

+2

Si vous êtes juste essayer de suivre une sinusoïde variant lentement, vous pourriez être mieux d'utiliser un (https [PLL (Phase Locked Loop)]: // fr .wikipedia.org/wiki/Phase-locked_loop # Implementing_a_digital_phase-locked_loop_in_software). –

+0

Sur une note de côté, vous pouvez regarder dans cette fonction, ['spectrogram'] (https://in.mathworks.com/help/signal/ref/spectrogram.html). Il calcule essentiellement le DFT pour les fenêtres coulissantes (également appelé "Transformée de Fourier à court terme"), au lieu d'essayer manuellement de mettre en œuvre le même. – crazyGamer

+1

J'ai réussi à faire fonctionner mon code, mais après quelques recherches, j'ai l'impression qu'une PLL est probablement plus applicable à ce que j'essaie de faire - je vais m'y pencher. –

Répondre

0

Il s'avère donc que mon problème était une mauvaise compréhension des mathématiques de la fonction sinus. J'avais supposé que la fréquence de l'onde était égale à tout ce qui était multiplié par la variable de temps (par exemple le f dans sin (ft)). Cependant, il s'avère que la fréquence est en fait définie par la dérivée de l'argument entier de la fonction sinus - le taux de changement de la phase.

Pour constant f les deux définitions sont égaux, puisque d (ft)/dt = f. Mais pour, disons, f (t) = sin (t):

d (f (t) t)/dt = d (sin (t) t)/dt = t cos (t) + sin (t)

La fréquence varie en fonction d'une fonction très différente de f (t). Modification de la définition de fonction fixe mon problème suivant:

f = @(t) 100-30*(t-0.5).^2; %Frequency function (Hz) 
G = cumsum(f(t))/fs; %Phase function (Hz) 
I = sin(2*pi*G); %Sampling function