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 à une valeur Inf. La fonction renvoie les pics dans l’ordre où ils se produisent et exclut les points d’extrémité du signal différents de Inf. 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.
data = [25 8 15 5 6 10 10 3 1 20 7];
Trouvez les maxima locaux. La fonction findpeaks renvoie les pics dans l’ordre où ils se produisent. Elle n’inclut pas les points d’extrémité de data dans les pics détectés, car ils n’y a pas de voisin de gauche ou de droite auxquels les comparer.
pks = findpeaks(data)
pks = 1×3
15 10 20
Utilisez findpeaks sans arguments de sortie pour afficher les pics. Pour le pic plat, la fonction met seulement en évidence le point dont l’indice est le plus faible.
findpeaks(data)
Créez un signal constitué d’une somme de courbes de Gauss. Définissez la position, 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; y = sum(Hgt.*(exp(-((x-Pos)./Wdt).^2)),1);
Utilisez findpeaks avec les paramètres par défaut pour trouver les pics du signal et leur emplacement.
[pks,locs] = findpeaks(y,x);
Tracez les pics avec findpeaks. Ajoutez une étiquette à chaque pic.
findpeaks(y,x) text(locs+.02,pks,num2str((1:numel(pks))'))
Triez les pics du plus haut au plus bas.
[psor,lsor] = findpeaks(y,x,SortStr="descend");
findpeaks(y,x)
text(lsor+.02,psor,num2str((1:numel(psor))'))
Créez un signal constitué d’une somme de courbes de Gauss 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; y = sum(Hgt.*(exp(-((x-Pos)./Wdt).^2)),1) + base;
Utilisez findpeaks pour localiser et tracer les pics dont la proéminence est supérieure ou égale à 4. Seuls le pic le plus haut et celui le plus bas remplissent la condition de proéminence minimale.
[~,locs4,~,proms4] = findpeaks(y,x,MinPeakProminence=4)
locs4 = 1×2
0.1982 0.4995
proms4 = 1×2
5.5773 4.4171
findpeaks(y,x,MinPeakProminence=4,Annotate="extents")
Affichez l’emplacement, la proéminence et la largeur à la moitié de la proéminence de tous les pics. Les deuxième et quatrième pics sont ceux dont la proéminence est supérieure ou égale à 4.
[~,locs,widths,proms] = findpeaks(y,x)
locs = 1×6
0.1001 0.1982 0.2983 0.4995 0.7017 0.8018
widths = 1×6
0.0154 0.0431 0.0377 0.0625 0.0274 0.0409
proms = 1×6
2.6816 5.5773 3.1448 4.4171 2.9191 3.6363
Extrayez les pics des données relatives au nombre moyen de taches solaires et estimez le cycle moyen des taches solaires.
Chargez le fichier sunspot.dat qui contient le nombre moyen de taches solaires observées chaque année de 1700 à 1987. Tracez le nombre de taches solaires par an et superposez-lui les maxima. Zoomez sur la période comprise entre 1855 et 1920.
load sunspot.dat year = sunspot(:,1); avSpots = sunspot(:,2); findpeaks(avSpots,year) xlabel("Year") xlim([1855 1920])
Les taches solaires sont un phénomène cyclique. On sait que leur nombre atteint un pic environ tous les 11 ans.
Estimez la durée du cycle en déterminant l’écart moyen entre les pics. Commencez par ignorer les pics qui ne sont pas des maxima. Limitez les écarts pic à pic acceptables à des périodes supérieures à 6 ans.
findpeaks(avSpots,year,MinPeakDistance=6)
xlabel("Year")
xlim([1855 1920])
Créez un tableau datetime avec les données annuelles. 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") xlabel("Year")
Calculez le cycle moyen des taches solaires. Le cycle moyen des taches solaires estimé est très proche de la durée connue, qui est de 11 ans.
dttmCycle = years(mean(diff(lk)))
dttmCycle = 10.9600
Il est possible de trouver les pics dont les caractéristiques répondent à certains critères, parmi lesquels la distance ou hauteur de pic minimale ou maximale, la proéminence et le seuil d’écart entre un pic et ses voisins. Utilisez findpeaks pour contraindre une ou plusieurs caractéristiques des pics.
Contraindre les pics selon une caractéristique d’un signal audio
Chargez un signal audio échantillonné à 7 418 Hz. Sélectionnez 200 échantillons.
load mtlb
select = mtlb(1001:1200);Trouvez et mettez en évidence tous les pics du signal audio. Trouvez ensuite les pics séparés par un écart d’au moins 5 ms et mettez-les en évidence. Pour appliquer la contrainte d’écart minimal, 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.
figure tiledlayout("vertical") nexttile findpeaks(select,Fs) legend("No Constraints",Location="southwest") nexttile lgn = @(txt) legend(txt,FontName="FixedWidth",Location="southwest"); findpeaks(select,Fs,MinPeakDistance=0.005) lgn("MinPeakDistance")
Trouvez et tracez les pics qui répondent individuellement à chacune des contraintes suivantes :
Trouvez les pics dont l’amplitude est supérieure ou égale à 1 V.
Trouvez les pics qui sont supérieurs d’au moins 1 V aux échantillons voisins.
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.
figure tiledlayout("vertical") nexttile findpeaks(select,Fs,MinPeakHeight=1) lgn("MinPeakHeight") nexttile findpeaks(select,Fs,Threshold=1) lgn("Threshold") nexttile findpeaks(select,Fs,MinPeakProminence=1) lgn("MinPeakProminence")
Contraindre plusieurs caractéristiques des pics d’un signal chirp
Trouvez et mettez en évidence tous les pics d’un signal chirp de 1 000 échantillons dont la largeur est comprise entre 2,5 et 4 ms.
N = 1000; x = sin(2*pi*(1:N)/N + (10*(1:N)/N).^2); figure findpeaks(x,Fs,Annotate="extents", ... MinPeakWidth=2.5e-3,MaxPeakWidth=4e-3) title("Peaks with Multiple Constraints")
Les capteurs peuvent renvoyer des lectures tronquées si les données dépassent un point de saturation donné. Vous pouvez considérer ces valeurs extrêmes 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². Échantillonnez le signal 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 = 100;
t = 0:1/Fs:1;
s = sin(2*pi*5*t).*sin(2*pi*3*t) + randn(size(t))/10;Synthétisez une mesure saturée en tronquant chaque lecture dont la valeur absolue est supérieure à la limite spécifiée de 0,32. Normalisez ensuite le signal pour qu’il soit compris entre 0,1 et 1. Convertissez l’amplitude normalisée en décibels et tracez le résultat.
bnd = 0.32; s(s>bnd) = bnd; s(s<-bnd) = -bnd; sNorm = 0.1+0.9*(s-min(s))/(max(s)-min(s)); sig = mag2db(sNorm); plot(t,sig) xlabel("Time (s)") ylabel("Magnitude (dB)") ylim([-25 5])
Localisez les pics et les creux du signal. findpeaks identifie uniquement le front montant de chaque pic plat.
[pk,lcp] = findpeaks(sig,t); % Peaks [vl,lcv] = findpeaks(-sig,t); % Valleys hold on plot(lcp,pk,"x",lcv,-vl,"*")
Utilisez l’argument nom-valeur Threshold pour exclure les creux et les pics plats. Imposez une différence d’amplitude minimale de 0,1 dB entre un pic ou un creux et ses voisins.
[pkt,lctp] = findpeaks(sig,t,Threshold=0.1); [vlt,lctv] = findpeaks(-sig,t,Threshold=0.1); plot(lctp,pkt,"o",LineWidth=1.2) plot(lctv,-vlt,"d",LineWidth=1.2)
Créez un signal constitué d’une somme de courbes de Gauss. 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 = [7 6 3 2 2 3]'; Wdt = [3 8 4 3 4 6]'/100; y = sum(Hgt.*(exp(-((x-Pos)./Wdt).^2)),1);
Mesurez la largeur des pics en utilisant la moitié de la proéminence et la moitié de la hauteur comme référence.
tiledlayout("flow") nexttile findpeaks(y,x,Annotate="extents") title("Half-Prominence Peak Widths") nexttile findpeaks(y,x,Annotate="extents",WidthReference="halfheight") title("Half-Height Peak Widths")
Sélectionnez les pics les plus élevés séparés par un écart d’au moins 0,5 unité sur l’axe des x. Mesurez la largeur des pics en utilisant la moitié de la proéminence et la moitié de la hauteur comme référence.
figure tiledlayout("flow") nexttile findpeaks(y,x,MinPeakDistance=0.5,Annotate="extents") title("Half-Prominence Peak Widths") nexttile findpeaks(y,x,MinPeakDistance=0.5,Annotate="extents", ... WidthReference="halfheight") title("Half-Height Peak Widths")
Seuls les premier et dernier pics remplissent la condition d’écart minimal. Les largeurs affichées sur le tracé correspondent donc à ces deux pics. Comme l’étendue de chaque pic reste inchangée, la largeur de pic conserve la même valeur indépendamment des conditions spécifiées et du fait que le pic soit sélectionné ou non.
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.
aTg = [1 1.5]; fTg = 1e3*[4.76 35.8]; freqkHzFull = (0:0.0025:50)'; waveFull = abs(sinc([1 0.5].*(freqkHzFull-fTg/1e3)))*aTg';
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"); table(PV,freq(PL), RowNames="Peak estimate "+(1:numel(PV)), ... VariableNames=["Amplitude" "Frequency (kHz)"])
ans=2×2 table
Amplitude Frequency (kHz)
_________ _______________
Peak estimate 1 1.4824 36
Peak estimate 2 0.9374 5
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); table(Ypk,Xpk, RowNames="Refined peak "+(1:numel(Ypk)), ... VariableNames=["Amplitude" "Frequency (kHz)"])
ans=2×2 table
Amplitude Frequency (kHz)
_________ _______________
Refined peak 1 1.5063 35.8
Refined peak 2 1.0163 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(aTg) % Theoretical peak amplitudes errorBounds = aTg.*(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)
