Matlab - Funcion findpeaks de matlab

   
Vista:

Funcion findpeaks de matlab

Publicado por KELVIN DE JESUS (1 intervención) el 10/06/2015 11:28:31
Hola que tal estimados compañeros, acudo a ustedes haber de que forma me pueden colaborar, resulta que me encuentro trabajando analisis de datos en un espectro de datos obtenidos mediante NMR, pasa que necesito identificar todo lo que represente picos, pero al momento de aplicar la funcion findpeaks a la data, este me arroja picos muchas veces donde en realidad no hay picos, quise probar con la misma funcion pero seleccionando un valor minimo, pero resulta que este tiende luego a ignorarme algunos picos que se encuentran por debajo del ruido pero que en realidad me representan informacion importante.

pueden decirme de que forma puedo modificar la funcion findpeaks para que realice esta operacion de manera optima?

este es el codigo de la funcion:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
function [pks,locs] = findpeaks(X,varargin)
%FINDPEAKS Find local peaks in data
%   PKS = FINDPEAKS(X) finds local peaks in the data vector X.
%
%   [PKS,LOCS]= FINDPEAKS(X) also returns the indices LOCS at which the
%   peaks occur.
%
%   [...] = FINDPEAKS(X,'MINPEAKHEIGHT',MPH) finds only those peaks that
%   are greater than MINPEAKHEIGHT MPH. Specifying a minimum peak height
%   may help in reducing the processing time. MPH is a real valued scalar.
%   The default value of MPH is -Inf.
%
%   [...] = FINDPEAKS(X,'MINPEAKDISTANCE',MPD) finds peaks that are at
%   least separated by MINPEAKDISTANCE MPD. MPD is a positive integer
%   valued scalar. This parameter may be specified to ignore smaller peaks
%   that may occur in close proximity to a large local peak. For example,
%   if a large local peak occurs at index N, then all smaller peaks in the
%   range (N-MPD, N+MPD) are ignored. If not specified, MPD is assigned a
%   value of one. 
%
%   [...] = FINDPEAKS(X,'THRESHOLD',TH)finds peaks that are at least
%   greater than their neighbhors by the THRESHOLD TH. TH is real valued
%   scalar greater than or equal to zero. The default value of TH is zero.
%
%   [...] = FINDPEAKS(X,'NPEAKS',NP) specifies the number of peaks to be
%   found. NP is an integer greater than zero. When specified, the
%   procedure terminates once NP peaks are found and the NP peaks are
%   returned. If not specified, all peaks are returned.
%
%   [...] = FINDPEAKS(X,'SORTSTR',STR) specifies the direction of sorting
%   of peaks. STR can take values of 'ascend','descend' or 'none'. If not
%   specified, STR takes the value of 'none' and the peaks are returned in
%   the order of their occurrence.
%
%   See also DSPDATA/FINDPEAKS
 
%   Copyright 2007-2009 The MathWorks, Inc.
%   $Revision: 1.1.6.10 $  $Date: 2009/10/16 06:41:09 $
 
error(nargchk(1,11,nargin,'struct'));
 
validateattributes(X,{'numeric'},{'nonempty','real','vector'},...
    'findpeaks','X');
 
% Check the input data type. Single precision is not supported.
try
    chkinputdatatype(X);
catch ME
    throwAsCaller(ME);
end
 
%#function dspopts.findpeaks
hopts = uddpvparse('dspopts.findpeaks',varargin{:});
 
Ph  = hopts.MinPeakHeight;
Pd  = hopts.MinPeakDistance;
Th  = hopts.Threshold;
Np  = hopts.NPeaks;
Str = hopts.SortStr;
 
pks = [];
locs = [];
 
 
M = numel(X);
 
if (M < 3)
    datamsgid = generatemsgid('emptyDataSet');
    error(datamsgid,'Data set must contain at least 3 samples.');
else
    Indx = find(X > Ph);
    if(isempty(Indx))
        mphmsgid = generatemsgid('largeMinPeakHeight');
        warning(mphmsgid,'Invalid MinPeakHeight. There are no data points greater than MinPeakHeight.');
    else
        % validate value of Pd and set default values for Pd and Np
        [Pd,Np] = setvalues(Pd,Np,M);
        if(Pd >= M)
            pdmsgid = generatemsgid('largeMinPeakDistance');
            error(pdmsgid,'Invalid MinPeakDistance. Set MinPeakDistance as an integer in the range between 1 and %s.',...
                num2str(M));
        else
            [pks,locs] =getpeaks(X,Indx,Pd,Th,Np);
        end
    end
end
 
if isempty(pks)
    npmsgid = generatemsgid('noPeaks');
    warning(npmsgid,'No peaks found.')
elseif(~strcmp(Str,'none'))
    [pks,s]  = sort(pks,Str);
    if(~isempty(locs))
        locs = locs(s);
    end
end
 
%--------------------------------------------------------------------------
function [Pd,Np] = setvalues(Pd,Np,L)
 
if ~isempty(Pd) && (~isnumeric(Pd) || ~isscalar(Pd) ||any(rem(Pd,1)) || (Pd < 1))
    Nmsgid = generatemsgid('invalidMinPeakDistance');
    error(Nmsgid,'MinPeakDistance should be an integer greater than 0.');
end
 
if(isempty(Pd))
    Pd = 1;
end
 
if(isempty(Np))
    Np = L;
end
 
%--------------------------------------------------------------------------
function [pks,locs] =getpeaks(Data,Indx,Pd,Th,Np)
% This function finds peaks in data set Data whose index set Indx is
% disjoint. Some data points were removed from the original set through
% preprocessing
 
m = 0;                  % counter for peaks found
L = numel(Indx);
M = numel(Data);
 
% Pre-allocate for speed
pks  = zeros(1,Np);
locs = zeros(1,Np);
 
endindx = M;      % last point in the index set
 
j = 0;
% First, the "Pd" neighbors, on either side, of the current data point are
% found. Then the current data point is compared with these neighbors to
% determine whether it is a peak.
 
while (j < L) && (m < Np)
    j = j+1;
 
    currentIdx = Indx(j);
 
    % Update the leftmost neighbor index if there is a peak within "Pd"
    % neighbors of leftmost index
    if(m > 0)
        prevPeakBoundR = min([locs(m)+Pd, endindx-1]);
        if currentIdx < prevPeakBoundR
            k = find(Indx(j+1:end) > prevPeakBoundR,1,'first');
            if ~isempty(k)
                j = j+k;
                currentIdx = Indx(j);
            else
                break;
            end
        end
    end
 
    % leftmost neighbor index
    endL = max(1,currentIdx - Pd);
 
    % rightmost neighbor index
    endR = min(currentIdx + Pd,endindx);
 
    % create neighbor data set
    temp = Data(endL:endR);
 
    % set current data point to -Inf in the neighbor data set
    temp(currentIdx-endL+1) = -Inf;
 
    % Compare current data point with all neighbors
    if(all((Data(currentIdx) > temp+Th)) && (currentIdx ~=1)&& (currentIdx~=M))
        m = m+1;
        locs(m) = currentIdx;  % loctions of peaks
        pks(m)  = Data(currentIdx);  % peaks
    end
end
 
% return all peaks found
if m~=0
    locs = locs(locs > 0);
    pks  = pks(1:length(locs));
else
    locs = [];
    pks = [];
end
 
 
% [EOF]
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder