ssa

1
function chss(pw)
2
    len   = length(pw); % êîëè÷åñòâî îòñ÷åòîâ pw
3
    cad = 30;      % 30 êàäðîâ/ñåê
4
%     cad=1427/60;
5
    dt  = 1.0/cad; % èíòåðâàë äèñêðåòèçàöèè âðåìåíè, ñåê
6
    tim=zeros(len,1);
7
    tim(1) = 0.0;
8
    for i=2:len
9
       tim(i) = tim(i-1)+dt; % âðåìÿ â ñåêóíäàõ
10
    end
11

12
    nPart = 20; % êîëè÷åñòâî äîëåé res
13
    win = 1024;
14
    fMi  = 40.0/60.0;   % ÷àñòîòà ñðåçà äëÿ 40 óä/ìèí (0.6667 Ãö)
15
    fMa  = 240.0/60.0;  % ÷àñòîòà ñðåçà äëÿ 240 óä/ìèí (4.0 Ãö)
16
%     Nf   = 1+win/2;     % êîë-âî îòñ÷åòîâ ÷àñòîòû
17
%     df   = cad/(win-1); % èíòåðâàë äèñêðåòèçàöèè ÷àñòîòû, Ãö 
18

19

20

21
    %%%%%%%%%%%%%%%%%%%%%%%%
22

23

24
    % Ìàòðèöà ïåðåêðûâàþùèõñÿ ñåãìåíòîâ pw
25
    col = 1; Imin = 1; Imax = win;
26
    res = len-win*floor(len/win); res = floor(res/nPart);
27
    LenArray = round((len-Imax)/res,-1);
28
    ns = zeros(LenArray, 1);
29
    tseg = zeros(LenArray, 1);
30
    while Imax<=len
31
       ns(col) = col; % íîìåð òåêóùåãî ñåãìåíòà pw
32
       Imin    = Imin+res;
33
       Imid  = floor(Imin+res/2);
34
       Imax    = Imax+res;
35
       tseg(col) = tim(Imid);
36
       col     = col+1;
37
    end
38
    col = col-1; % êîë-âî ïåðåêðûâàþùèõñÿ ñåãìåíòîâ â ïðåäåëàõ len
39
%     NSF = win+res*(col-1); % íîìåð ôèíàëüíîãî îòñ÷åòà ôèíàëüíîãî ñåãìåíòà <=len
40
    
41
    spw=zeros(win, col);
42
    for j=1:col
43
       for i=1:win
44
          k        = (j-1)*res;
45
          spw(i,j) = pw(k+i); % òåêóùèé ñåãìåíò pw äëèííîþ win 
46
       end
47
    end
48

49
    % Îöåíêè ÑÏÌ ïåðåêðûâàþùèõñÿ ñåãìåíîâ pw
50
    f(1) = 0.0; df = cad/(win-1); % èíòåðâàë äèñêðåòèçàöèè ÷àñòîòû, Ãö
51
    row  = 1+win/2;
52
    Fmin = fMi-10*df; Fmax = fMa+10*df; % ÷àñòîòà â Ãö
53
    for i=2:row
54
       f(i) = f(i-1)+df; % ÷àñòîòà â ãåðöàõ
55
       if abs(f(i)-Fmin)<=df
56
          iGmin = i;
57
       end
58
       if abs(f(i)-Fmax)<=df
59
          iGmax = i;
60
       end
61
    end
62
    for i=1:iGmax
63
       fG(i) = f(i); % ñåòêà ÷àñòîò 3D-ãðàôèêà
64
    end
65
    f = f'; smopto = 3; % ïàðàìåòð ñãëàæèâàíèÿ ïåðèîäîãðàììû Òîìñîíà
66
    figure();
67
    zpg_spw=zeros(iGmax, col);
68
    zto_spw=zeros(iGmax, col);
69
    for j=1:col
70
       pg_spw(:,j) = periodogram(spw(:,j),blackmanharris(win),win);
71
       stem(f(iGmin:iGmax),pg_spw(iGmin:iGmax,j),'LineStyle','-','Marker','none'); hold on;
72
       pto_spw(:,j) = pmtm(spw(:,j),smopto,win);   % ïåðèîäîãðàììà Òîìñîíà
73
       for i=1:iGmax
74
          zpg_spw(i,j) = pg_spw(i,j);
75
          zto_spw(i,j) = pto_spw(i,j);
76
       end
77
    end
78
    grid on; title('ÑÏÌ spw');
79
    xlabel("f,Hz",'interp','none'); ylabel("Psd",'interp','none');
80
    
81
    
82
ppto=load('Files/pto.txt')';
83
figure(); 
84
subplot(1,2, 1); imagesc(ppto);
85
subplot(1,2, 2); imagesc(pto_spw);
86
% pto_spw=ppto;
87

88

89
    % Îöåíêè ñðåäíèõ ÷àñòîò îñíîâíîãî òîíà äëÿ ñåãìåíòîâ pw
90
    for j=1:col
91
        [~,I] = sort(pg_spw(:,j),'descend');
92
        pg_fMAX(j) = f(I(1)); % I(1) - èíäåêñ ÷àñòîòû(Ãö) ìàêñèìóìà pg_spw(:,j)   
93
        [~,I] = sort(pto_spw(:,j),'descend');
94
        pto_fMAX(j) = f(I(1));
95
     end
96
     pto_fMAX = pto_fMAX';
97
     smo_pto_fMAX = smooth(pto_fMAX,0.3*col,'rloess'); 
98
%       EUT_P1LC3_Fmaxcr = smo_pto_fMAX; 
99
%       save('EUT_P1LC3_Fmaxcr','EUT_P1LC3_Fmaxcr');
100

101
     figure();
102
    % plot(ns,pg_fMAX,'b','Marker','o'); hold on;
103
     plot(pto_fMAX,'g'); hold on; % tseg
104
     plot(smo_pto_fMAX,'k','LineWidth',0.8);
105
     grid off; title('×àñòîòû îñíîâíîãî òîíà ïåðèîäîãðàììû Òîìñîíà spw ');
106
     xlabel("t, s",'interp','none'); ylabel("fMAX,Ãö",'interp','none');
107
end