ssa

1
function chss2(pw, Path, Name)
2
%     dt=1/30;     % èíòåðâàë âðåìåííîé äèñêðåòèçàöèè
3
    fmin=0.15; % íèæíÿÿ ãðàíèöà - ÷àñòîòíûé äèàïàçîí äûõàòåëüíîé âîëíû
4
    dt=1/30;     % èíòåðâàë âðåìåííîé äèñêðåòèçàöèè
5
    Nmed=1/(dt*fmin); % àïåðòóðà ôèëüòðà
6
    PathHR=Path+replace(Name,"_nc","")+"_rPPG_output.csv";
7
    hr = LoadHR(PathHR);
8
    % PathHR=Path+replace(Name,"_nc","")+"_Mobi_RR-intervals.rr";
9
    % hr = LoadRR(PathHR);
10

11
    %% Cåãìåíòû pw
12
    N   = length(pw); % êîëè÷åñòâî îòñ÷åòîâ pw
13
    win = 1024;
14
    res = N-win*floor(N/win);
15
    nPart = 20; % êîëè÷åñòâî äîëåé res
16
    res = floor(res/nPart); overlap = (win-res)/win;
17
    S = 1; Imin = 1; Imax = win;
18
    while Imax<=N
19
        ns(S) = S; % íîìåð òåêóùåãî ñåãìåíòà pw
20
        Imin  = Imin+res;
21
        Imax  = Imax+res;
22
        S     = S+1;
23
    end
24
    S = S-1; % êîë-âî ïåðåêðûâàþùèõñÿ ñåãìåíòîâ pw â ïðåäåëàõ N
25
    NSF = win+res*(S-1); % íîìåð ôèíàëüíîãî îòñ÷åòà ôèíàëüíîãî ñåãìåíòà <= N
26
    for j=1:S
27
        for i=1:win
28
            k = (j-1)*res;
29
            spw(i,j) = pw(k+i); % òåêóùèé ñåãìåíò pw äëèííîþ win 
30
        end
31
    end
32
    %% Set general parameters
33
    cad = 30;      % 30 êàäðîâ/ñåê
34
    dt  = 1.0/cad; % èíòåðâàë äèñêðåòèçàöèè âðåìåíè, ñåê
35
    tim(1) = 0.0;
36
    for i=2:N
37
        tim(i) = tim(i-1)+dt; % âðåìÿ â ñåêóíäàõ
38
    end
39
    % tim = tim';
40
    ns  = (1:S)'; % íîìåðà ñåãìåíòîâ pw
41
    % fmp=zeros(S,1);
42
    for j=1:S % öèêë ïî ñåãìåíòàì pw
43
        %    L(j) = floor(cad/fmp(j)); % êîë-âî îòñ÷åòîâ îñíîâíîãî òîíà pw
44
        L(j) = floor(cad/1.5); % êîë-âî îòñ÷åòîâ îñíîâíîãî òîíà pw
45
    end
46
    L = L';
47
    K = 5; % êîë-âî ïåðèîäîâ äëÿ ïàðàìåòðà âëîæåíèÿ
48
    M = K*max(L); % ïàðàìåòð âëîæåíèÿ â òðàåêòîðíîå ïðîñòðàíñòâî 
49
    %% SSA- àíàëèç ñåãìåíòîâ pw
50
    nET = 4;   % êîë-âî ñèíãóëÿðíûõ òðîåê äëÿ ñåãìåíòîâ pw
51
    for j=1:S  % öèêë ïî ñåãìåíòàì
52
        %% SSA time series
53
        %    M = K*L(j); % ïàðàìåòð âëîæåíèÿ â òðàåêòîðíîå ïðîñòðàíñòâî
54
        [C,LBD,RC] = SSA(win,M,spw(:,j),nET);
55
        %% Estimation of the spw(:,j) reconstructed with the pair of RC   
56
        sET12(:,j) = sum(RC(:,1:2),2);   
57
        %    sET34(:,j) = sum(RC(:,3:4),2);   
58
        %% Compare reconstruction and original time series
59
        if j==100 % íîìåð ñåãìåíòà pw äëÿ âèçóàëèçàöèè
60
            figure('name','Covariance matrix'); clf;
61
            imagesc(C); axis square; set(gca,'clim',[-1 1]); colorbar;
62
            figure('name','Eigenvalues'); clf; plot(LBD,'o-');
63
            figure('name','Original time series and reconstruction'); clf;
64
            plot(tim(1:win),spw(:,j),'b-',tim(1:win),sET12(:,j),'r-');
65
            legend('Original','sET12'); xlabel("t,s",'interp','none'); ylabel("sET",'interp','none');
66
        end
67
    end
68
    %% Îöåíêà ÀÊÔ ñèíãóëÿðíûõ òðîåê äëÿ ñåãìåíòîâ pw
69
    lag  = floor(win/10); % íàèáîëüøèé ëàã ÀÊÔ <= win/10
70
    lagS = 2*lag;
71
    for j=1:S
72
        Acf_sET12(:,j) = AcfMed(lagS,win,sET12(:,j)); % íîðìèðîâàííàÿ ÀÊÔ j-ãî ñåãìåíòà
73
        %    Acf_sET12(:,j) = autocorr(sET12(:,j),'NumLags',lag); % íîðìèðîâàííàÿ ÀÊÔ j-ãî ñåãìåíòà
74
    end
75
    %% Âèçóàëèçàöèÿ ÀÊÔ ñèíãóëÿðíûõ òðîåê äëÿ ñåãìåíòîâ pw
76
    lgl = 1:lag; % ñåòêà 3D-ãðàôèêà ÀÊÔ
77
    Time=0:dt:lag*dt-dt;
78
    figure('name','ÀÊÔ ñèíãóëÿðíûõ òðîåê sET12 ñåãìåíòîâ pw'); clf;
79
    % mesh(ns,lgl,Acf_sET12(1:lag,:),'FaceAlpha',0.5,'FaceColor','flat'); colorbar;
80
    % xlabel("ns",'interp','none'); ylabel("lag",'interp','none');
81
    mesh(ns,Time,Acf_sET12(1:lag,:),'FaceAlpha',0.5,'FaceColor','flat'); colorbar;
82
    xlabel("ns",'interp','none'); ylabel("lag,s",'interp','none');
83
    zlabel("Acf",'interp','none'); grid on;
84
    %% Îãèáàþùàÿ ïî êðèòåðèþ ëîêàëüíûõ ìàêñèìóìîâ abs(acf_sET12)
85
    for j=1:S % öèêë ïî ñåãìåíòàì ÀÊÔ
86
        absTS = abs(Acf_sET12(:,j));
87
        AT1   = absTS(1);
88
        AT2   = absTS(2);
89
        maxTS = zeros(lag,1); maxTS(1) = AT1;
90
        maxN  = zeros(lag,1); maxN(1)  = 1;
91
        Nmax = 1;
92
        for m=3:lag
93
            AT3 = absTS(m);
94
            if (AT1<=AT2)&&(AT2>=AT3)
95
                Nmax = Nmax+1; % íîìåð î÷åðåäíîãî óçëà èíòåðïîëÿöèè (ñ÷åò÷èê ìàêñèìóìîâ)
96
                maxN(Nmax) = m-1; % íîìåð î÷åðåäíîãî ìàêñèìóìà äëÿ ðÿäà absTS
97
                maxTS(Nmax) = AT2; % îòñ÷åò î÷åðåäíîãî óçëà èíòåðïîëÿöèè
98
            end
99
            AT1 = AT2;
100
            AT2 = AT3;
101
        end
102
        Nmax = Nmax+1; % êîëè÷åñòâî óçëîâ èíòåðïîëÿöèè
103
        maxN(Nmax)  = lag; % íîìåð îòñ÷åòà absTS ôèíàëüíîãî óçëà èíòåðïîëÿöèè
104
        maxTS(Nmax) = absTS(lag); % îòñ÷åò absTS ôèíàëüíîãî óçëà èíòåðïîëÿöèè
105
        NumMax = maxN(1:Nmax); % íîìåðà ìàêñèìóìîâ ÂÐ absTS
106
        % Èíòåðïîëÿöèÿ îãèáàþùåé ÀÊÔ
107
        % 'pchip','cubic','v5cubic','makima','spline'
108
        EnvAcf_sET12(:,j) = interp1(NumMax,maxTS(1:Nmax),lgl,'pchip');
109
        AcfNrm_sET12(:,j) = Acf_sET12(1:lag,j)./EnvAcf_sET12(:,j); % íîðìèðîâàííûå ÀÊÔ
110
    end
111
    figure('name','Íîðìèðîâàííûå ÀÊÔ ñèíãóëÿðíûõ òðîåê sET12 ñåãìåíòîâ pw'); clf;
112
    % mesh(ns,lgl,AcfNrm_sET12(1:lag,:),'FaceAlpha',0.5,'FaceColor','flat'); colorbar;
113
    % xlabel("ns",'interp','none'); ylabel("lag",'interp','none');
114
    mesh(ns,Time,AcfNrm_sET12(1:lag,:),'FaceAlpha',0.5,'FaceColor','flat'); colorbar;
115
    xlabel("ns",'interp','none'); ylabel("lag,s",'interp','none');
116
    zlabel("Acf_Nrm",'interp','none'); grid on;
117
    %% Ìãíîâåííàÿ ÷àñòîòà íîðìèðîâàííîé ÀÊÔ ñèíãóëÿðíûõ òðîåê sET12 äëÿ ñåãìåíòîâ pw
118
    pi2 = 2.0*pi;
119
    for j=1:S % öèêë ïî ñåãìåíòàì ÀÊÔ
120
        PhaAcfNrm = abs(acos(AcfNrm_sET12(:,j))); % ìãíîâåííàÿ ôàçà
121
        pAcf = pchip(lgl,PhaAcfNrm); % êîýôôèöèåíòû èíòåðïîëÿíòà pchip
122
        for m=2:lag
123
            FrcAcfNrm(m) = abs(pAcf.coefs(m-1,3))/pi2/dt; % ìãíîâåííàÿ ÷àñòîòà íîðìèðîâ-îé ÀÊÔ, Ãö        
124
        end
125
        FrcAcfNrm(1) = FrcAcfNrm(2); 
126
        %    FrcAcfNrm = abs(diff(PhaAcfNrm))/pi2/dt; % ìãíîâåííàÿ ÷àñòîòà íîðìèðîâ-îé ÀÊÔ, Ãö
127
        insFrc_AcfNrm(j) = median(FrcAcfNrm); % ñðåäíÿÿ ìãíîâåííàÿ ÷àñòîòòà j-ãî ñåãìåíòà pw 
128
    end
129
    
130
    disp("Àïåðòóðà ôèëüòðà äëÿ insFrc_AcfNrm: " + Nmed);  
131
%     insFrc_AcfNrm=medfilt1(insFrc_AcfNrm, 5);
132
      
133
    smo_insFrc_AcfNrm = smoothdata(insFrc_AcfNrm, 'rloess', 0.25*S); % smo_insFrc_AcfNrm = smooth(insFrc_AcfNrm,0.25*S,'rlowess');
134
    figure('name','×àñòîòû íîðìèð-îé ÀÊÔ ñèíãóëÿð-õ òðîåê ñåãìåíòîâ pw','Position', [0 0 1400 800]); clf;
135
%             insFrc_AcfNrm=medfilt1(insFrc_AcfNrm,Nmed);
136

137
    p1 = plot(ns,insFrc_AcfNrm,'b','LineWidth',0.8); hold on;
138
    plot(ns,smo_insFrc_AcfNrm,'r','LineWidth',0.8); grid on;
139
    xlabel("ns",'interp','none'); ylabel("insFrc_AcfNrm,Hz",'interp','none');
140
    title("×àñòîòû íîðìèð-îé ÀÊÔ ñèíãóëÿð-õ òðîåê ñåãìåíòîâ pw");
141
%     legend(p1,'sET12');
142
    if length(hr)>100
143
        ns_hr = (length(ns)/length(hr) : length(ns)/length(hr) : length(ns))';
144
        % yyaxis right; 
145
        plot(ns_hr,hr./60,'black'); ylabel("HR[bpm]",'interp','none');
146
        % legend(p1,'insFrc_AcfNrm','rloess','HR[bpm]');
147
        hr_med=medfilt1(hr,Nmed*5);
148
        hr_diff_med=hr-hr_med;
149
        plot(ns_hr,hr_med./60,'cyan--');
150
        plot(ns_hr,hr_diff_med./60,'magenta'); %ylabel("HR[bpm]",'interp','none');
151
        legend('insFrc AcfNrm','rloess','HR','HR[medfilt]','HR[HR-medfilt1]')
152
    end
153
    
154
    %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
155
    
156
    
157
    
158
%         [outs, lower_prct, upper_prct] =  RaznFilter(insFrc_AcfNrm, [1 99]);
159
    subplot(1,3,[1 2]);
160
    plot(ns,insFrc_AcfNrm,'b','LineWidth',0.8); hold on;
161
    plot(ns,smo_insFrc_AcfNrm,'r','LineWidth',0.8); grid on;
162
    
163
    % Ýòî ïîðîã äëÿ ïåðâîé ðàçíîñòè ìîæåò áûòü êàê ÷èñëîì, òàê è
164
    % ïðîöåíòèëÿìè êð÷
165
    threshold = [1 99];
166
    threshold = 0.05;
167
    insFrc_AcfNrm_first_diff_filter = first_diff_filter(insFrc_AcfNrm, threshold);
168
    plot(ns, insFrc_AcfNrm_first_diff_filter  ,'magenta-x' ); grid on;
169
    
170
%     plot(ns,[outs outs(length(outs))],'red-o' ); grid on;
171
%     line('XData', [0 length(outs)], 'YData', [upper_prct upper_prct], 'Color','black','LineStyle','--');
172
%     line('XData', [0 length(outs)], 'YData', [lower_prct lower_prct], 'Color','black','LineStyle','--');
173
    if length(hr)>100
174
        ns_hr = (length(ns)/length(hr) : length(ns)/length(hr) : length(ns))';
175
        plot(ns_hr,hr./60,'yellow'); ylabel("HR[bpm]",'interp','none');
176
        
177
        hr_fd = first_diff_filter(hr, 4.5);
178
        plot(ns_hr,hr_fd./60,'black'); ylabel("HR+FD[bpm]",'interp','none');
179
        
180
    end  
181

182
    xlabel("ns",'interp','none'); ylabel("insFrc_AcfNrm,Hz",'interp','none');
183
    title("×àñòîòû íîðìèð-îé ÀÊÔ ñèíãóëÿð-õ òðîåê ñåãìåíòîâ pw");
184
    
185
    
186
    subplot(1,3,3);
187
    diff_signal = diff(insFrc_AcfNrm);
188
    lower_prct=0; upper_prct=0;
189
    if length(threshold) == 1
190
        lower_prct=-threshold; upper_prct=threshold;
191
    end
192
    if length(threshold) == 2 
193
        prct = prctile(diff_signal, threshold);
194
        lower_prct=prct(1); upper_prct=prct(2);
195
    end
196
    
197
    plot(diff_signal,'red-o'); hold on; grid on;
198
    line('XData', [0 length(diff_signal)], 'YData', [lower_prct lower_prct], 'Color','black','LineStyle','--');
199
    line('XData', [0 length(diff_signal)], 'YData', [upper_prct upper_prct],'Color','black','LineStyle','--');
200
    
201
    
202
%     legend(p1,'sET12');
203
%     if length(hr)>100
204
%         ns_hr = (length(ns)/length(hr) : length(ns)/length(hr) : length(ns))';
205
%         % yyaxis right; 
206
%         plot(ns_hr,hr./60,'black'); ylabel("HR[bpm]",'interp','none');
207
%         % legend(p1,'insFrc_AcfNrm','rloess','HR[bpm]');
208
%         hr_med=medfilt1(hr,Nmed*5);
209
%         hr_diff_med=hr-hr_med;
210
%         plot(ns_hr,hr_med./60,'cyan--');
211
%         
212
%         legend('insFrc AcfNrm','rloess','HR','HR[medfilt]'); 
213
% 
214
%         
215
%         
216
%     subplot(1,3,3); 
217
%         
218
%         plot(outs,'red-o');
219
%         plot((insFrc_AcfNrm-medfilt1(insFrc_AcfNrm, 5)),'black--'); hold on;
220
% %         plot(rmoutliers_emulated(insFrc_AcfNrm, [10 90])./60,'red'); 
221
% %         plot(filloutliers(insFrc_AcfNrm,"next")./60,'red'); 
222
%         
223
% %         plot(outs, 'blue');
224
%         line('XData', [0 200], 'YData', [lower_prct lower_prct], 'Color','red','LineStyle','--');
225
%         line('XData', [0 200], 'YData', [upper_prct upper_prct],'Color','red','LineStyle','--');
226
% 
227
% 
228
%         
229
% %         ðàçíîñòü ïî èñõîäíûì äàííûì êîíòàêòíîãî ñïîñîáà îïðåäåëåíèÿ
230
% %         plot(ns_hr,hr_diff_med./60,'magenta'); %ylabel("HR[bpm]",'interp','none');
231
%         legend('HR[HR-medfilt1]');  ylabel("HR[bpm]",'interp','none'); xlabel("ns",'interp','none'); grid on;
232
%     end
233

234
    %% Îöåíêè ÑÏÌ ñèíãóëÿðíûõ òðîåê äëÿ ñåãìåíîâ pw
235
    smopto = 3; % ïàðàìåòð ñãëàæèâàíèÿ ïåðèîäîãðàììû Òîìñîíà
236
    for j=1:S
237
        %    disp(spw(:,j));
238
        pto_sET12(:,j) = periodogram(spw(:,j),blackmanharris(win),win); % Áëýêìàíà-Õàððèñà
239
        pto_sET12(:,j) = pmtm(sET12(:,j),smopto,win); % ïåðèîäîãðàììà Òîìñîíà
240
    end
241
    %% Âèçóàëèçàöèÿ ÑÏÌ ñèíãóëÿðíûõ òðîåê ñåãìåíîâ pw
242
    fmi  = 40.0/60.0;   % ÷àñòîòà ñðåçà äëÿ 40 óä/ìèí (0.6667 Ãö)
243
    fma  = 240.0/60.0;  % ÷àñòîòà ñðåçà äëÿ 240 óä/ìèí (4.0 Ãö)
244
    Nf   = 1+win/2;     % êîë-âî îòñ÷åòîâ ÷àñòîòû
245
    df   = cad/(win-1); % èíòåðâàë äèñêðåòèçàöèè ÷àñòîòû, Ãö
246
    Fmin = fmi-10*df; Fmax = fma+10*df; % ÷àñòîòà â Ãö
247
    f(1) = 0.0;
248
    for i=2:Nf
249
        f(i) = f(i-1)+df; % ÷àñòîòà â ãåðöàõ
250
        if abs(f(i)-Fmin)<=df
251
            iGmin = i;
252
        end
253
        if abs(f(i)-Fmax)<=df
254
            iGmax = i;
255
        end
256
    end
257
    for i=1:iGmax
258
        fG(i) = f(i); % ñåòêà ÷àñòîò 3D-ãðàôèêà
259
    end
260
    f = f';
261
    figure('name','Ïåðèîäîãðàììà Òîìñîíà sET12 ñåãìåíòîâ pw'); clf;
262
    mesh(ns,fG(iGmin:iGmax),pto_sET12(iGmin:iGmax,:),'FaceAlpha',0.5,'FaceColor','flat');
263
    colorbar; grid on;
264
    xlabel("ns",'interp','none'); ylabel("f,Hz",'interp','none'); zlabel("P(f)",'interp','none');
265
    %% Îöåíêè ñðåäíèõ ÷àñòîò îñíîâíîãî òîíà ñèíãóëÿðíûõ òðîåê ñåãìåíòîâ pw
266
    for j=1:S
267
        [B,I] = sort(pto_sET12(:,j),'descend');
268
        pto_fMAX12(j) = f(I(1)); % I(1) - èíäåêñ ÷àñòîòû(Ãö) ìàêñèìóìà pto_sET12(:,j)
269
    end
270
    pto_fMAX12 = pto_fMAX12';
271
    smo_pto_fMAX12 = smoothdata(pto_fMAX12,'rloess',0.3*S); 
272
    % smo_pto_fMAX12 = smooth(pto_fMAX12,0.3*S,'rloess');
273
    figure('name','×àñòîòû îñíîâíîãî òîíà sET ñåãìåíòîâ pw','Position', [800 0 800 600]); clf;
274
%             pto_fMAX12=medfilt1(pto_fMAX12,Nmed);
275
    p=plot(ns,pto_fMAX12,'b'); hold on;
276
    plot(ns,smo_pto_fMAX12,'r','LineWidth',0.8); grid on;
277
    xlabel("ns",'interp','none'); ylabel("fMAX,Hz",'interp','none');
278
    title("×àñòîòû îñíîâíîãî òîíà sET ñåãìåíòîâ pw");
279
    if length(hr)>100
280
        ns_hr = (length(ns)/length(hr) : length(ns)/length(hr) : length(ns))';
281
        % yyaxis right;
282
        plot(ns_hr,hr./60,'black'); ylabel("HR[bpm]",'interp','none');
283
        % legend(p,'pto_fMAX12','rloess','HR[bpm]');
284
        hr_med=medfilt1(hr,Nmed*5);
285
        hr_diff_med=hr-hr_med;
286
        plot(ns_hr,hr_med./60,'cyan--');
287
        plot(ns_hr,hr_diff_med./60,'magenta'); %ylabel("HR[bpm]",'interp','none');
288
        legend('pto sET12','smoothdata','HR[bpm]','medfilt1','hr-medfilt1')
289
    end
290
    saveas(p,Path+Name+'_×ÑÑ_sET.png');
291

292
    %% Àãðåãèðîâàíèå ñåãìåíòîâ î÷èùåííîé ïóëüñîâîé âîëíû cpw
293
    [NumS,cpw_avr,cpw_med,cpw_iqr] = wav(NSF,S,win,res,sET12);
294
    % figure('name','Pulse wave'); clf;
295
    % plotwave(1,NSF,tim,cpw_avr,cpw_med,cpw_iqr);
296
    %% Íàêîïëåííàÿ ìãíîâåííàÿ ôàçà cpw
297
    % cpw = cpw_avr; % îöåíêà î÷èùåííîé ïóëüñîâîé âîëíû
298
    cpw = cpw_med; % îöåíêà î÷èùåííîé ïóëüñîâîé âîëíû
299
    cutoff = pi; pi2 = 2.0*pi;
300
    H_cpw = hilbert(cpw);
301
    insE_cpw = abs(H_cpw); % ìãíîâåííàÿ îãèáàþùàÿ
302
    unwPha = unwrap(angle(H_cpw),cutoff); % íàêîïëåííàÿ ìãíîâåííàÿ ôàçà
303
    % Íåïðåðûâíàÿ-(ñ) è ðàçðûâíàÿ-(d) êîìïîíåíòû íàêîïëåííîé ìãíîâåííîé ôàçû
304
    unwPc_cpw(1) = unwPha(1); unwPd_cpw(1) = 0.0;
305
    for i=2:NSF
306
        dif = unwPha(i)-unwPha(i-1);
307
        unwPc_cpw(i) = unwPc_cpw(i-1); % íåïðåðûâíàÿ
308
        unwPd_cpw(i) = unwPd_cpw(i-1); % ðàçðûâíàÿ 
309
        if dif>=0.0
310
            unwPc_cpw(i) = unwPc_cpw(i)+dif;
311
        else
312
            unwPd_cpw(i) = unwPd_cpw(i)+dif+pi2;
313
        end
314
    end
315
    unwPc_cpw = unwPc_cpw'; unwPd_cpw = unwPd_cpw';
316
    figure('name','Unwrape phase pulse wave'); clf;
317
    sp1 = subplot(2,1,1); plot(tim(1:NSF),unwPc_cpw); grid on;
318
    xlabel("t,s",'interp','none'); ylabel("Phase cont",'interp','none');
319
    title(sp1,'Íåïðåðûâíàÿ íàêîïëåííàÿ ôàçà pw');
320
    sp2 = subplot(2,1,2); plot(tim(1:NSF),unwPd_cpw); grid on; 
321
    xlabel("t,s",'interp','none'); ylabel("Phase disc",'interp','none');
322
    title(sp2,'Ðàçðûâíàÿ íàêîïëåííàÿ ôàçà pw');
323
    %% Ìãíîâåííàÿ ÷àñòîòà è ýíåðãèÿ î÷èùåííîé ïóëüñîâîé âîëíû
324
    % Îöåíêà ïåðâîé ïðîèçâîäíîé íåïðåðûâíîé êîìïîíåíòû íàêîïëåííîé ìãíîâåííîé ôàçû
325
    % ñ ïîìîùüþ ñîõðàíÿþùåé ôîðìó êóñî÷íî êóáè÷åñêîé èíòðåðïîëÿöèè ïîëèíîìàìè Ýðìèòà
326
    t = 1:NSF;
327
    p = pchip(t,unwPc_cpw);
328
    insF_cpw(1) = 0.0;
329
    for i=2:NSF
330
        insF_cpw(i) = p.coefs(i-1,3)/pi2/dt; % ìãíîâåííàÿ ÷àñòîòà íåïðåðûâíîé êîìïîíåíòû cpw, Ãö
331
    end
332
    insF_cpw(1) = insF_cpw(2);
333
    insF_cpw = insF_cpw';
334
    smo_insF_cpw = smoothdata(insF_cpw, 'rloess', 0.03*NSF); 
335
    % smo_insF_cpw = smooth(insF_cpw,0.03*NSF,'rloess');
336
    res_insF_cpw = insF_cpw-smo_insF_cpw; % êâàäðàò îñòàòêà 
337
    dev_insF_cpw = smoothdata(res_insF_cpw.^2, 'rloess', 0.03*NSF); 
338
    % dev_insF_cpw = smooth(res_insF_cpw.^2,0.03*NSF,'rloess');
339
    std_insF_cpw = abs(sqrt(dev_insF_cpw));
340
    figure('name','Frequencie and energy pulse wave'); clf;
341
    sp1 = subplot(2,1,1); plot(tim(1:NSF),insF_cpw); hold on;
342
    plot(tim(1:NSF),smo_insF_cpw,'Color','r','LineWidth',0.8); hold on;
343
    % plot(tim(1:NSF),std_insF_cpw);
344
    xlabel("t,s",'interp','none'); ylabel("insF,Hz",'interp','none');
345
    grid on; title(sp1,'Ìãíîâåííàÿ ÷àñòîòà cpw'); ylim([1.0 3.0]);
346
    sp2 = subplot(2,1,2); plot(tim(1:NSF),insE_cpw.^2); 
347
    xlabel("t,s",'interp','none'); ylabel("insE^2",'interp','none');
348
    grid on; title(sp2,'Ìãíîâåííàÿ ýíåðãèÿ cpw');
349
    save(Path+Name+"_nc"+".mat");
350
end