findpeaks
Trouver les maxima locaux
Syntaxe
Description
pks = findpeaks(y)y. Un pic local est un échantillon de données qui est plus grand que les deux échantillons voisins ou bien égal à Inf. La sortie répertorie les pics dans l’ordre où ils se produisent. Les points d’extrémité du signal différents de Inf sont exclus. Si un pic est plat, la fonction renvoie seulement le point dont l’indice est le plus faible.
[___] = findpeaks(___, spécifie des options avec des arguments nom-valeur en plus du ou des arguments d’entrée des syntaxes précédentes.Name,Value)
findpeaks(___) sans arguments de sortie trace le signal en le superposant aux valeurs de pic.
Exemples
Définissez un vecteur contenant trois pics et tracez-le.
data = [25 8 15 5 6 10 10 3 1 20 7]; plot(data)
Trouvez les maxima locaux. La sortie répertorie les pics dans l’ordre où ils se produisent. Le premier échantillon n’est pas inclus bien qu’il s’agisse du maximum. Pour le pic plat, la fonction renvoie seulement le point dont l’indice est le plus faible.
pks = findpeaks(data)
pks = 1×3
15 10 20
Utilisez findpeaks sans arguments de sortie pour afficher les pics.
findpeaks(data)
Créez un signal constitué d’une somme de courbes en cloche. Définissez l’emplacement, la hauteur et la largeur de chaque courbe.
x = linspace(0,1,1000); Pos = [1 2 3 5 7 8]/10; Hgt = [3 4 4 2 2 3]; Wdt = [2 6 3 3 4 6]/100; for n = 1:length(Pos) Gauss(n,:) = Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2); end PeakSig = sum(Gauss);
Tracez chacune des courbes ainsi que leur somme.
plot(x,Gauss,'--',x,PeakSig)
Utilisez findpeaks avec les paramètres par défaut pour trouver les pics du signal et leur emplacement.
[pks,locs] = findpeaks(PeakSig,x);
Tracez les pics avec findpeaks et ajoutez-leur des étiquettes.
findpeaks(PeakSig,x) text(locs+.02,pks,num2str((1:numel(pks))'))
Triez les pics du plus haut au plus bas.
[psor,lsor] = findpeaks(PeakSig,x,'SortStr','descend'); findpeaks(PeakSig,x) text(lsor+.02,psor,num2str((1:numel(psor))'))
Créez un signal constitué d’une somme de courbes en cloche sur une période entière d’un cosinus. Définissez l’emplacement, la hauteur et la largeur de chaque courbe.
x = linspace(0,1,1000); base = 4*cos(2*pi*x); Pos = [1 2 3 5 7 8]/10; Hgt = [3 7 5 5 4 5]; Wdt = [1 3 3 4 2 3]/100; for n = 1:length(Pos) Gauss(n,:) = Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2); end PeakSig = sum(Gauss)+base;
Tracez chacune des courbes ainsi que leur somme.
plot(x,Gauss,'--',x,PeakSig,x,base)
Utilisez findpeaks pour localiser et tracer les pics dont la proéminence est supérieure ou égale à 4.
findpeaks(PeakSig,x,'MinPeakProminence',4,'Annotate','extents')
Le pic le plus haut et celui le plus bas sont les seuls à remplir cette condition.
Affichez la proéminence et la largeur à la moitié de la proéminence de tous les pics.
[pks,locs,widths,proms] = findpeaks(PeakSig,x); widths
widths = 1×6
0.0154 0.0431 0.0377 0.0625 0.0274 0.0409
proms
proms = 1×6
2.6816 5.5773 3.1448 4.4171 2.9191 3.6363
Les taches solaires sont un phénomène cyclique. On sait que leur nombre atteint un pic environ tous les 11 ans.
Chargez le fichier sunspot.dat qui contient le nombre moyen de taches solaires observées chaque année de 1700 à 1987. Trouvez et tracez les maxima.
load sunspot.dat
year = sunspot(:,1);
avSpots = sunspot(:,2);
findpeaks(avSpots,year)
Affinez votre estimation de la durée du cycle en ignorant les pics qui sont très proches les uns des autres. Trouvez et tracez à nouveau les pics, cette fois en limitant les écarts pic à pic acceptables à des valeurs supérieures à six ans.
findpeaks(avSpots,year,'MinPeakDistance',6)
Utilisez les emplacements des pics renvoyés par findpeaks pour calculer l’intervalle moyen entre les maxima.
[pks,locs] = findpeaks(avSpots,year,'MinPeakDistance',6);
meanCycle = mean(diff(locs))meanCycle = 10.9600
Créez un tableau datetime avec les données relatives à l’année. Supposez que les taches solaires ont été comptées le 20 mars de chaque année, vers l’équinoxe de printemps. Trouvez les années où le nombre de taches solaires atteint un pic. Utilisez la fonction years pour définir l’écart minimal entre les pics comme une donnée de type duration.
ty = datetime(year,3,20); [pk,lk] = findpeaks(avSpots,ty,'MinPeakDistance',years(6)); plot(ty,avSpots,lk,pk,'o')
Calculez le cycle moyen des taches solaires avec la fonctionnalité datetime.
dttmCycle = years(mean(diff(lk)))
dttmCycle = 10.9600
Créez une timetable avec les données. Définissez la variable temporelle en années. Tracez les données. Affichez les cinq dernières entrées de la timetable.
TT = timetable(years(year),avSpots); plot(TT.Time,TT.Variables)
entries = TT(end-4:end,:)
entries=5×1 timetable
Time avSpots
________ _______
1983 yrs 66.6
1984 yrs 45.9
1985 yrs 17.9
1986 yrs 13.4
1987 yrs 29.3
Chargez un signal audio échantillonné à 7 418 Hz. Sélectionnez 200 échantillons.
load mtlb
select = mtlb(1001:1200);Trouvez les pics séparés par un écart d’au moins 5 ms.
Pour appliquer cette contrainte, findpeaks choisit le plus haut pic du signal et ignore tous les pics situés à moins de 5 ms de celui-ci. La fonction répète ensuite la procédure pour le plus haut pic restant et recommence jusqu’à ce qu’il n’y ait plus de pics à traiter.
findpeaks(select,Fs,'MinPeakDistance',0.005)
Trouvez les pics dont l’amplitude est supérieure ou égale à 1 V.
findpeaks(select,Fs,'MinPeakHeight',1)
Trouvez les pics qui sont supérieurs d’au moins 1 V aux échantillons voisins.
findpeaks(select,Fs,'Threshold',1)
Trouvez les pics qui descendent d’au moins 1 V d’un côté ou de l’autre avant que le signal atteigne une valeur supérieure.
findpeaks(select,Fs,'MinPeakProminence',1)
Les capteurs peuvent renvoyer des lectures tronquées si les données dépassent un point de saturation donné. Vous pouvez considérer ces pics comme inutiles et choisir de les ignorer ou au contraire les intégrer dans votre analyse.
Générez un signal constitué du produit de fonctions trigonométriques de fréquences 5 Hz et 3 Hz noyé dans un bruit blanc gaussien avec une variance de 0,1². Le signal est échantillonné pendant une seconde à une fréquence de 100 Hz. Réinitialisez le générateur de nombres aléatoires pour obtenir des résultats reproductibles.
rng default
fs = 1e2;
t = 0:1/fs:1-1/fs;
s = sin(2*pi*5*t).*sin(2*pi*3*t)+randn(size(t))/10;Simulez une mesure saturée en tronquant chaque lecture supérieure à une limitée donnée de 0,32. Tracez le signal saturé.
bnd = 0.32;
s(s>bnd) = bnd;
plot(t,s)
xlabel('Time (s)')
Localisez les pics du signal. findpeaks renvoie uniquement le front montant de chaque pic plat.
[pk,lc] = findpeaks(s,t); hold on plot(lc,pk,'x')
Utilisez l’argument nom-valeur 'Threshold' pour exclure les pics plats. Imposez une différence d’amplitude minimale de entre un pic et ses voisins.
[pkt,lct] = findpeaks(s,t,'Threshold',1e-4); plot(lct,pkt,'o','MarkerSize',12)
Créez un signal constitué d’une somme de courbes en cloche. Définissez l’emplacement, la hauteur et la largeur de chaque courbe.
x = linspace(0,1,1000); Pos = [1 2 3 5 7 8]/10; Hgt = [4 4 2 2 2 3]; Wdt = [3 8 4 3 4 6]/100; for n = 1:length(Pos) Gauss(n,:) = Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2); end PeakSig = sum(Gauss);
Tracez chacune des courbes ainsi que leur somme.
plot(x,Gauss,'--',x,PeakSig)
grid
Mesurez la largeur des pics en utilisant la moitié de la proéminence comme référence.
findpeaks(PeakSig,x,'Annotate','extents')
Mesurez à nouveau les largeurs, cette fois en utilisant la moitié de la hauteur comme référence.
findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight') title('Signal Peak Widths')
Affinez l’estimation de l’emplacement des deux principaux pics d’un signal en utilisant la méthode des moindres carrés non linéaires avec un noyau de fonction sinc.
Générer le signal
La compression d’impulsion radar d’une forme d’onde FM linéaire produit un spectre en sinus cardinal (sinc) où les emplacements de fréquence des pics sont proportionnels à la distance entre le radar et l’objet détecté. Vous pouvez commencer par estimer l’emplacement et l’amplitude des pics avec findpeaks puis affiner vos estimations avec refinepeaks. Cet exemple montre comment déterminer l’amplitude et l’emplacement des pics d’un signal de compression d’impulsion synthétique non bruité. Il utilise refinepeaks pour affiner les estimations.
Générez un signal composé de deux formes d’onde en sinus cardinal (sinc) avec des pics de 1 et 1,5 respectivement situées à 4,76 kHz et 35,8 kHz. Définissez l’espacement des fréquences à 2,5 Hz.
aTarget = [1 1.5]; fTarget = 1e3*[4.76 35.8]; freqkHzFull = (0:0.0025:50)'; waveFull = abs(sinc([1 0.5].*(freqkHzFull-fTarget/1e3)))*aTarget';
Sous-échantillonnez le signal par un facteur de 200 pour que l’espacement des fréquences entre les échantillons soit de 0,5 kHz. Cet exemple affine les estimations d’amplitude et d’emplacement des pics du signal sous-échantillonné et compare les estimations affinées aux valeurs des signaux d’origine.
freq = downsample(freqkHzFull,200); wave = downsample(waveFull,200); plot(freqkHzFull,waveFull,freq,wave,"*") legend(["Full signal" "Selected samples"],Location="northwest") xlabel("Frequency (kHz)") ylabel("Magnitude")
Affiner les pics avec les moindres carrés non linéaires
Utilisez findpeaks pour obtenir une estimation initiale des amplitudes, emplacements et largeurs à mi-hauteur des deux plus hauts pics du signal.
[PV,PL,PW] = findpeaks(wave,NPeaks=2, ... SortStr="descend",WidthReference="halfheight");
Utilisez refinepeaks pour affiner l’estimation des pics avec la méthode NLS (moindres carrés non linéaires). Définissez les points de fréquence du signal et les largeurs de pics. Les valeurs de pics sont sensiblement plus proches des valeurs attendues de 1,5 et 1 tandis que les emplacements de fréquence se rapprochent respectivement de 35,8 kHz et 4,76 kHz.
LW = max(PW,2);
[Ypk,Xpk] = refinepeaks(wave,PL,freq,Method="NLS",LobeWidth=LW)Ypk = 2×1
1.5063
1.0163
Xpk = 2×1
35.8001
4.7628
Tracez les amplitudes des pics affinées sur l’axe des y et les emplacements mis à jour des pics par rapport aux estimations initiales sur l’axe des x. Les deux pics initialement estimés et les deux échantillons voisins sont à 0,5 kHz d’écart. Les pics affinés signalés par des cercles pleins indiquent les emplacements réels des pics par rapport à ceux initialement estimés ainsi que les amplitudes corrigées.
refinepeaks(wave,PL,freq,Method="NLS",LobeWidth=LW) yline(aTarget) % Theoretical peak amplitudes errorBounds = aTarget.*(1+0.03*[-1;1]); yline(errorBounds(:),":") % ±3% error bounds legend("Peak "+[1 2])
Arguments d'entrée
Arguments nom-valeur
Arguments de sortie
En savoir plus
Conseils
Vous pouvez réaliser une estimation initiale des pics du signal avec findpeaks puis affiner leurs amplitudes et leurs emplacements avec refinepeaks.
Supposons qu’un signal ait les amplitudes y et les emplacements x. L’extrait de code suivant montre comment estimer et affiner les pics à partir de y et x.
[yPeaks,xPeaksIdx] = findpeaks(y); [yRPeaks,xRPeaks] = refinepeaks(y,xPeaksIdx,x)
