Pourquoi les valeurs de fréquence sont-elles arrondies dans le signal à l'aide de la FFT?

Puisque vous souhaitez obtenir une résolution de 0,01 Hz, vous devrez échantillonner au moins 100 secondes de données. Vous pourrez résoudre des fréquences jusqu'à environ 22,05 kHz.

Les intervalles de résultats DFT sont séparés par Fs / N en fréquence, où N est la longueur de la FFT. Ainsi, la durée de votre fenêtre DFT limite la résolution en termes d'espacement des centres de fréquence des bin de résultats DFT.

Mais, pour des pics de fréquence bien séparés à faible bruit (S / N élevé), au lieu d'augmenter la durée des données, vous pouvez à la place estimer les emplacements des pics de fréquence à une résolution plus élevée en interpolant le résultat DFT entre le résultat DFT bacs. Vous pouvez essayer l'interpolation parabolique pour une estimation grossière de l'emplacement du pic de fréquence, mais l'interpolation Sinc fenêtrée (essentiellement la reconstruction de Shannon-Whittaker) fournirait une bien meilleure précision et résolution d'estimation de fréquence (étant donné un plancher de bruit suffisamment bas autour du ou des pics de fréquence d'intérêt, par exemple, pas de sinusoïdes à proximité dans votre cas de forme d'onde artificielle).

En effet, si la trame dure T secondes, les fréquences de la DFT sont k / T Hz, où k est un entier. En conséquence, le suréchantillonnage n'améliore pas la précision de la fréquence estimée, tant que ces fréquences sont identifiées comme des maxima de l'amplitude de la DFT. Au contraire, considérer des trames plus longues d'une durée de 100 s induirait un espacement entre les fréquences DFT de 0,01 Hz, ce qui pourrait être assez bon pour produire la fréquence attendue. Il est possible de faire bien mieux, en estimant la fréquence d'un pic comme sa fréquence moyenne par rapport à la densité de puissance.

Figure 1: même après application d'une fenêtre de Tuckey, la DFT du signal fenêtré n'est pas une somme de Dirac: il y a encore des fuites spectrales en bas des pics. Cette puissance doit être prise en compte lors de l'estimation des fréquences.

Un autre problème est que la longueur de la trame n'est pas un multiple de la période du signal, qui peut de toute façon ne pas être périodique. Néanmoins, la DFT est calculée comme si le signal était périodique mais discontinu au bord de la trame. Il induit des fréquences aiguës décrites comme des fuites spectrales . Le fenêtrage est la méthode de référence pour traiter de tels problèmes et atténuer le problème lié à la discontinuité artificielle. En effet, la valeur d'une fenêtre décroît continuellement jusqu'à zéro près des bords du cadre. Il existe une liste de fonctions de fenêtre et de nombreuses fonctions de fenêtre sont disponibles dans scipy.signal . Une fenêtre est appliquée comme suit:

import numpy as np
from scipy import signal
import scipy

sr = 44100 # sample rate
x = np.linspace(0, 1, sr) # one second of signal
tpi = 2 * np.pi
data = np.sin(261.63 * tpi * x) + np.sin(329.63 * tpi * x) + np.sin(392.00 * tpi * x)

#a window...
tuckey_window=signal.tukey(len(data),0.5,True)
data=data*tuckey_window

data -= np.mean(data)
fft = np.fft.rfft(data, norm="ortho")

def abs2(x):
        return x.real**2 + x.imag**2

fftmag=abs2(fft)[:1000]
peaks, _= signal.find_peaks(fftmag, height=np.max(fftmag)*0.1)
print "potential frequencies ", peaks

#compute the mean frequency of the peak with respect to power density
powerpeak=np.zeros(len(peaks))
powerpeaktimefrequency=np.zeros(len(peaks))
for i in range(1000):
    dist=1000
    jnear=0
    for j in range(len(peaks)):
        if dist>np.abs(i-peaks[j]):
             dist=np.abs(i-peaks[j])
             jnear=j
    powerpeak[jnear]+=fftmag[i]
    powerpeaktimefrequency[jnear]+=fftmag[i]*i


powerpeaktimefrequency=np.divide(powerpeaktimefrequency,powerpeak)
print 'corrected frequencies', powerpeaktimefrequency

À ce stade, les fréquences présentant la plus grande magnitude sont encore 262, 330 et 392. L'application d'une fenêtre ne fait que rendre les pics plus visibles: La DFT du signal fenêtré comporte trois crêtes distinctes, chacune comportant un lobe central et des lobes latéraux, en fonction de la DFT de la fenêtre. Les lobes de ces fenêtres sont symétriques: la fréquence centrale peut donc être calculée comme la fréquence moyenne du pic, par rapport à la densité de puissance.

tuckey_window=signal.tukey(len(data),0.5,True)
data=data*tuckey_window

Les fréquences estimées résultantes sont de 261,6359 Hz, 329,637 Hz et 392,0088 Hz: c'est bien mieux que 262, 330 et 392 Hz et il satisfait la précision requise de 0,01 Hz pour un signal d'entrée aussi pur et silencieux.